-
嫦娥二号(CE-2) 卫星是我国探月工程任务实施的第二颗环月轨道运行卫星。它搭载了一台两线阵CCD立体相机，以线阵推扫成像方式沿卫星飞行方向成像，可获取前视与后视两个线阵影像。与CE-1号相比，CE-2号影像地面分辨率由120 m提高到7 m，月表地貌细节更加丰富，数据量亦大大增加。
由于月面地物类型单一，与地球遥感影像相比，月球影像存在反照率低、对比度弱、纹理信息匮乏等问题。此外，由于成像时刻太阳高度角的原因，月球影像上有较多撞击坑内壁存在太阳直射过曝光或处于阴影区等而导致纹理丢失。在现有的研究中，目前比较流行的SIFT (scale invariant feature transform)、SURF(speeded up robust features)等特征算法应用于CE-1号影像匹配可以获得一定数量的稀疏匹配点[1, 2]；邹崇尧等[3]采用基于RPC近似核线约束的影像匹配方法对CE-1号影像进行了试验；李春来等[4]利用基于SIFT点特征与线特征相结合的匹配算法获取了密集匹配点，并制作了CE-1号500 m分辨率全月球DEM (digital elevation model)产品。随着影像分辨率的提高，生成三维地形数据要求更加密集的同名点。文献[5, 6]对CE-2号CCD影像提取DEM方法都进行了研究，采用特征匹配方法获取一定数量同名点后，再利用ENVI软件提取DEM数据，该方法虽然能够提取出月表DEM信息，但从结果上看存在依靠特征匹配无法获取稠密的同名点, DEM精度受限、自动化程度不高, 无法应用于大范围数据等实际生产等问题。
针对CE-2号影像的特点，本文提出了一种适用于CE-2号影像的提取月表三维地形信息的方法，具体流程如下。
1) 基于RPC(rational polynomial coefficients)模型生成近似核线影像，在此基础上构建金字塔影像，以金字塔上层影像匹配结果预测下层匹配点的视差初值，实现了原始影像上的逐像素匹配；
2) 根据匹配点坐标和RPC模型计算其月面坐标，再采用基于反距离插值方法插值生成三维地形；
3) 对于图像阴影或过曝光的地形漏洞区域，采用基于径向基函数的方法进行插值填补。
本文提出的方法实现了逐像素匹配，能够获取稠密的同名点以及自动提取出月表三维地形，可用于大范围数据的任务生产。
-
CCD立体相机图像数据首先经过了数据预处理[7]，并在此基础上开展后期处理。月球影像存在明显的地域性特征特点：在高纬地区图像清晰、层次丰富，月表地物界限明显；而在低纬度区域及月海区域，影像的亮度和对比度较低，图像质量相对较差。此外，嫦娥二号图像数据在星上多经过了8倍压缩，对图像质量有一定的影响。因此，在进行影像匹配前，根据月球影像数据的特点进行预处理工作，比如增强图像的对比度、减小图像的噪声对匹配的影响等。本文采用不同的滤波算法进行了试验，结果表明，采用Wallis滤波进行影像预处理可以得到较好的结果。Wallis滤波[8]是一种局部的影像变换，它使影像反差小的区域反差增大，影像反差大的区域反差减小，使得影像中灰度的微小变化信息得到增强。该方法在计算影像的局部灰度均值和方差时引入了平滑算子，亦能压制噪声，改善影像的质量。
-
目前对卫星影像核线关系的各种描述中，基于投影轨迹法的核线定义建立在成像的几何约束条件之上，其理论最为严密[9, 10]。经过左影像上p点的成像光线的每一个物点，投影到右影像上的投影点轨迹将形成一条曲线，将这条曲线定义为像点p的核线。类似于框幅式成像的核线几何关系，点p在右像上的同名点p′必定位于这条“核曲线”上。已有研究表明，该核曲线具有小范围内的近似直线性和局部共轭性[11]。
在实际的单轨近似核线影像生成时，通过RPC模型建立地面点与前后视影像上的像素坐标的关系，如式(1)、式(2)，前视影像与后视影像RPC文件中分别有80个多项式系数。
$$ ({\rm{Lat, Lon}}) = f(S, L, {\rm{Height}}) $$ (1) $$ \left( {S, L} \right) = f\prime ({\rm{Lat, Lon, Height}}) $$ (2) 式中，f和f′代表原始影像RFM坐标正解和反解函数; (S, L)为标准化后的像点坐标; (Lat, Lon, Height)为标准化后的地面点大地坐标。
-
金字塔影像匹配策略是一项应用较广泛的匹配策略，其优点如下[3]。
1) 采用了低通滤波使得高层影像在保留了原始影像大的整体结构特征时，去除了非平稳信息和噪声，减少了匹配时的不确定性；
2) 金字塔上层影像匹配的结果作为下层影像匹配的初始值和约束，减少了搜索范围时，又提高了匹配的可靠性。构建金字塔影像的方法有移动平均法、高斯滤波法等，本文在近似核线影像基础上采用高斯滤波法构建了3层金字塔。
影像匹配时，采用归一化互相关系数(NCC)作为匹配测度：
$$ \begin{array}{l} \rho \left( {c, r} \right) = \\ \;\frac{{\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^n {({g_{i, j}} \times {{g'}_{i + r, j + c}})}- \frac{1}{{mn}}\left( {\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^n {{g_{i, j}}} } } \right)\left( {\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^n {{{g'}_{i + r, j + c}}} } } \right)} }}{{\sqrt {\left[{\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^n {g_{i, j}^2-\frac{1}{{mn}}{{\left( {\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^n {{g_{i, j}}} } } \right)}^2}} } } \right]\left[{\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^n {{{g'}_{i + r, j + c}}} }-\frac{1}{{mn}}{{\left( {\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^n {{{g'}_{i + r, j + c}}} } } \right)}^2}} \right]} }} \end{array} $$ (3) 当ρ=ρmax(c, r)且ρmax(c, r)>ρthreshold时，接受该像素作为待匹配点的同名点，对于一维相关时, r=0。ρthreshold从金字塔顶层到下层可逐步减小。
在对每层金字塔影像匹配时，以前视影像为基准，在后视影像上搜索待匹配点，搜索范围d可定义为：
$$ d = \left\{ \begin{array}{l} {\rm{width}}/2, {\rm{ }}i = 1\\ 2({u_{i-1}}-1) \sim 2({u_{i-1}} + 1), {\rm{ }}i = 2, 3, 4 \end{array} \right. $$ (4) 式中，ui－1为该像元对应上一层影像的视差值，i表示金字塔影像的层数。对于金字塔顶层影像，为了加快处理速度，在前视影像按照M×N的间隔定义了一个规则格网，只对格网点在后视影像上搜索待匹配点。而对于其他层影像，则进行逐像元匹配。
待所有点匹配结束后，需要对匹配结果进行视差一致性检查，剔除视差变化较大的错误匹配点。视差一致性检查的思路是：假设该点视差u(x, y)与其相邻点视差的差值为Δu=|u(x, y)－u(x±1, y±1)|，若共有m个点的Δu < Δuthreshold则接受该点，否则剔除该点。视差值检查后，对没有匹配成功的像素采用双性内插方法得出其视差值。最后，对于第4层(即原始影像层)，原始图像在图像阴影和过曝光的区域的点，如果处在该区域，说明匹配结果是不可靠的，删除这些点的视差值。
-
影像匹配计算得到视差文件后，可以根据视差文件前视影像上每一个像元，得到对应后视影像上像素的行列号坐标。利用RPC模型计算每个三维点坐标，得到在物方的一些三维散点，进而插值生成DEM规则格网，其中每个格网点的高程值通过反距离加权法(inverse distance weighted，IDW)计算得到。反距离加权插值算法是基于相近相似的原理，每个采样点都对插值点具有一定的影响，即权重。权重随着采样点和插值点之间距离的增加而减弱，距离插值点越近的采样点的权重越大。在任一插值点处的值是各采样点权重之和，表示为[12]：
$$ {z_p} = \sum\limits_{i = 1}^n {(d_i^{-u} \times {z_i})} /\sum\limits_{i = 1}^n {d_i^{-u}} $$ (5) 式中，zp为插值点的高程值; di为采样点到插值点的距离; u为权指数(常取值1或2)。以待插值点为中心，5×DEM输出分辨率为半径内的点作为采样点参与计算。
-
在地形漏洞填补算法的选择上，经过试验对比，本文选择了基于径向基函数(radial basis function，RBF)的曲面拟合算法来进行漏洞的填补。由于撞击坑边缘的地势变化往往比较缓慢而均匀，这种算法适合于解决撞击坑边缘的插值问题。RBF是一个随着距某一位置的距离变化而变化的函数，可以表示为[13]：
$$ F\left( x \right) = \sum\limits_{i = 1}^N {{w_i}\varphi } \left( {\left\| {x-{x_i}} \right\|} \right) $$ (6) 式中，φ(‖x－xi‖)是径向基函数；wi是权值。
常用的径向基函数有Multiquadrics、Inverse Multiquadrics、Gauss函数等不同形式，本文采用的是Multiquadrics函数，其数学表达形式：
$$ \varphi \left( r \right) = ({r^2} + {c^2})1/2, c > 0, 且r \in R $$ (7) 实际处理时，以漏洞的外接矩形中心为圆心，以漏洞的外接矩形的短边像素数加上一定外扩像素数作为半径划圆，将圆内有效格网点值参与插值计算。为加快数据处理速度，将该圆内的点抽稀取样10×10的格网，格网点的有效值参与计算。曲面插值计算完后，在原有的阴影边界处再进行了羽化接边处理，保证DEM插值补漏洞的连续性。
-
月球的形貌特征主要分为高地和月海，本文在高地和月海分别选取了一轨数据(轨道号分别为035 308和059 512)，从数据处理效率、影像匹配结果和地形构建结果准确性等方面对本文方法进行评价。本文在HP Z800图形工作站(Inter(R) Xeon (R) CPU X5670 @ 2.93 GHz (2处理器))上进行试验，将算法进行了多核并行处理以提高处理效率，整个数据处理时间见表 1。从表 1中可以看出，匹配处理和地形构建处理均耗时较短，可用于大范围的数据生产任务中。
表 1 本文方法数据处理耗时情况
Table 1. Computational Time of Image Processing
轨道号 像素个数/
(像数×像数)影像匹配
/s地形构建
/s合计
/s035 308 6 144×39 760 140 90 230 059 512 6 144×33 980 135 85 220 本文对最终匹配结果的灰度相关性进行了统计。一般而言，在图像匹配处理中，生成用于相对定向的连接点时，对点的精度要求较高，灰度相关系数往往大于0.7。而在生成三维地形数据时的同名点时，灰度相关系数大于0.5以上的点可以接受为同名点，在纹理匮乏的地区如沙漠、屋顶等系数还可适当降低。从表 2中可以看出，高地地区纹理细节丰富，图像质量较好，灰度相关系数大于0.5占总像素的98.6%，匹配成功率较高；月海地区纹理匮乏，图像亮度对比度低，灰度相关系数大于0.5占总像素的77.3%，低于高地地区图像，但由于月海地区地形起伏小、特征地物少，该匹配点数已足以表达该地区地形特征。此外，采用的匹配策略是匹配点精度由高到低、从大到小逐渐缩小同名点搜索区域，灰度相关系数低于0.5的也有较大部分亦是可靠的，实际数据生产时是可适当调整参数。
表 2 影像匹配结果灰度相关性统计/%
Table 2. The Correlation of Grey Scale BetweenImage Matching Points/%
轨道号 占总像素数的百分比 >0.5 >0.6 >0.7 >0.8 >0.9 035 308 98.6 96.3 89.3 71.8 32.3 059 512 77.3 64.6 47.9 29.5 10.4 图 1是位于月球南半球正面的第谷撞击坑(中心经纬度为43°23′S, 11°10′W)的影像和本文处理后得到的20 m分辨率DEM数据的色彩渲染图。从图 1中可以看出，月表的地形起伏、地形地貌特征都得到了较好的表现，真实反映了实际的月表地形。取经过撞击坑内外的一条高程剖面线，将该数据与CE-1号500 mDEM(数据来源：http:// moon.bao.ac.cn/)、美国30 m LOLA (数据来源：http://imbrium.mit.edu/DATA/)进行比较，对比结果如图 2(a)。从图 2中可以看出，本文方法得到的数据与CE-1号500 m DEM数据的地形起伏趋势吻合，与美国30 m分辨率的LOLA数据相比非常吻合。但其在地形细节上较LOLA数据更加丰富，如图 2(b)、2(c)、2(d)中本文制作的DEM数据在撞击坑的边缘表现更加丰富，表达出了LOLA数据丢失的较小的地形起伏信息。
图 1 第谷撞击坑原始影像与本文DEM色彩渲染图
Figure 1. The Original Image and DEM Shaded Relief Map Produced by Our Method in Tycho Crater
图 2 本文制作的DEM数据与CE-1、LOLA数据的高程对比图
Figure 2. The Elevation Comparisons Between our DEM, CE-1 DEM and LOLA DEM
此外，本文还对图像阴影区域或者过曝光的区域进行了插值填补，由于该区域图像信息丢失，通过影像匹配无法获得同名点，因此, 只能通过周围像素的地形信息插值填补。图 3(a)是一个撞击坑内壁处在阴影区域的影像，本文经过插值填补后得到了该撞击坑的DEM数据，从图 3(b)中可以看出，碗型撞击坑的形态特征完整，处于阴影位置处的撞击坑内壁地形无明显异常。本文将其与相同位置的LOLA数据进行了比较，取经过地形漏洞划一条高程剖面线，从图 3(c)上可以看出，本文地形插值算法得到的高程值与LOLA数据相差很小，说明本文的地形漏洞插值填补算法是有效的。但任何插值算法都存在一定局限性，本文采用的插值算法在填补面积较大阴影区域造成的地形漏洞时难免会产生误差，需根据实际情况加以调整。
图 3 本文地形漏洞填补后DEM数据与LOLA数据的高程对比图
Figure 3. The Elevation Comparisons Between the DEM After Terrain Hole Filling and LOLA
-
摘要: 嫦娥二号卫星搭载了一台两线阵CCD立体相机，以线阵推扫成像方式获取了覆盖全月球的7 m分辨率影像。与嫦娥一号影像相比，影像分辨率大大提高，月表地貌细节更加丰富，制作高分辨率月表DEM更加复杂而困难，现有的基于嫦娥影像制作方法难以满足应用需求。针对嫦娥二号影像特点，在基于RPC模型生成近似核线影像基础上，进行金字塔影像匹配，实现了原始影像的逐像素匹配；计算匹配点的月面坐标后，采用基于反距离插值方法生成月表三维地形；对于图像阴影等地形漏洞，采用基于径向基函数进行插值填补。试验结果表明，该方法生成的月表三维地形DEM数据，能准确真实地表达高分辨率的月表地形地貌细节，并较好地应用于嫦娥二号全月球三维地形DEM数据产品的生产任务中。Abstract: The two-line CCD stereo camera carried by the Chang E-2 (CE-2) satellite has successfully captured the global lunar imagery at 7 m resolution by the linear push-broom imaging manner. Compared with the image of Chang E-1(CE-1), the spatial resolution of CE-2 is greatly improved, which enables the more detailed presentation of the lunar surface features and topography, and consequently makes the production of high-resolution lunar DEM more complex and difficult. The existing methods based on the Chang E image can hardly meet the application requirements. In this paper, based on the previous studies on the characteristic of CE-2 image, we proposed methods of image matching and the Digital Elevation Model (DEM) extraction. The working flow started with the generation of the approximate epipolar images based on the RPC model and the construction of the pyramid images. By matching with the pyramid hierarchical, a pixel-wise matching of the original image was realized. We then calculated the three-dimension coordinates on the lunar surface of the matching points, and generated the DEM using the inverse distance interpolation method. Finally, the radial basis function interpolation was employed to fill the terrain hole caused by the image shadows. The results show that the lunar DEM generated by our method can accurately depict the most detailed and complete lunar surface topography, and perform better in the production of the global Lunar DEM products of CE-2 mission.
-
Key words:
- Chang E-2 /
- image matching /
- DEM extraction /
- Moon
-
图 1 第谷撞击坑原始影像与本文DEM色彩渲染图
Figure 1. The Original Image and DEM Shaded Relief Map Produced by Our Method in Tycho Crater
图 2 本文制作的DEM数据与CE-1、LOLA数据的高程对比图
Figure 2. The Elevation Comparisons Between our DEM, CE-1 DEM and LOLA DEM
图 3 本文地形漏洞填补后DEM数据与LOLA数据的高程对比图
Figure 3. The Elevation Comparisons Between the DEM After Terrain Hole Filling and LOLA
表 1 本文方法数据处理耗时情况
Table 1. Computational Time of Image Processing
轨道号 像素个数/
(像数×像数)影像匹配
/s地形构建
/s合计
/s035 308 6 144×39 760 140 90 230 059 512 6 144×33 980 135 85 220 
下载: 导出CSV
表 2 影像匹配结果灰度相关性统计/%
Table 2. The Correlation of Grey Scale BetweenImage Matching Points/%
轨道号 占总像素数的百分比 >0.5 >0.6 >0.7 >0.8 >0.9 035 308 98.6 96.3 89.3 71.8 32.3 059 512 77.3 64.6 47.9 29.5 10.4
下载: 导出CSV
-
[1] 贾博, 姜挺, 郭丽, 等. SIFT和相关系数在嫦娥一号月球影像匹配中的应用[J].测绘科学技术学报, 2009, 26(5):363-366 http://www.cnki.com.cn/Article/CJFDTOTAL-JFJC200905015.htm Jia Bo, Jiang Ting, Guo Li, et al.Application of SIFT and Correlation Coefficient in Chang E-1 Lunar Image Matching[J].Journal of Geomatics Science and Technology, 2009, 26(5):363-366 http://www.cnki.com.cn/Article/CJFDTOTAL-JFJC200905015.htm [2] 何钰, 蓝朝桢, 邢帅.一种基于嫦娥一号CCD影像的影像匹配方法[J].测绘科学技术学报, 2012, 29(4):281-284 http://www.cnki.com.cn/Article/CJFDTOTAL-JFJC201204012.htm He Yu, Lan Chaozhen, Xing Shuai.A Method of Image Matching Based on CE-1 CCD Image[J].Journal of Geomatics Science and Technology, 2012, 29(4):281-284 http://www.cnki.com.cn/Article/CJFDTOTAL-JFJC201204012.htm [3] 邹崇尧, 何培培, 赵双明.基于RPC近似核线约束的CE-1影像匹配方法[J].地理空间信息, 2013, 11(3):1-3 http://www.cnki.com.cn/Article/CJFDTOTAL-DXKJ201303000.htm Zou Chongyao, He Beibei, Zhao Shuanming. Research on CE-1 Image Matching Method Based on RPC Approximate Epipolar Constraint[J]. Geospatial Information, 2013, 11(3):1-3 http://www.cnki.com.cn/Article/CJFDTOTAL-DXKJ201303000.htm [4] 李春来.嫦娥一号三线阵CCD数据摄影测量处理及全月球数字地形图[J].测绘学报, 2013, 42(6):853-860 http://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201306012.htm Li Chunlai. Photogrammetric Processing and Lunar Global Topographic Map from the Chang'E-13 Line-array CCD Data[J].Acta Geodaetica et Cartographica Sinica, 2013, 42(6):853-860 http://www.cnki.com.cn/Article/CJFDTOTAL-CHXB201306012.htm [5] 张亭. 嫦娥二号线阵CCD影像的月表DEM提取技术研究[D]. 成都: 成都理工大学, 2012 http://cdmd.cnki.com.cn/Article/CDMD-10616-1012500545.htm Zhang Ting. Research on the Lunar DEM Extraction Technology of Based Chang'E-2 Line Array CCD Image[D]. Chengdu:Chengdu University of Technology, 2012 http://cdmd.cnki.com.cn/Article/CDMD-10616-1012500545.htm [6] 尹永会. 嫦娥二号CCD影像DEM提取方法及月表地形分析初步研究[D]. 成都: 成都理工大学, 2013 http://cdmd.cnki.com.cn/Article/CDMD-10616-1013288925.htm Yin Yonghui. Research on the Lunar DEM Extraction Method and Terrain Analysis Technology Based on CE-2 Image[D]. Chengdu:Chengdu University of Technology, 2013 http://cdmd.cnki.com.cn/Article/CDMD-10616-1013288925.htm [7] 刘建军, 任鑫, 谭旭, 等.嫦娥二号CCD立体相机数据预处理与数据质量评价[J].武汉大学学报·信息科学版, 2013, 38(2):186-190 http://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201302014.htm Liu Jianjun, Ren Xin, Tan Xu, et al. Lunar Image Data Preprocessing and Quality Evaluation of CCD Stereo Camera on Chang'E-2[J]. Geomatics and Information Science of Wuhan University, 2013, 38(2):186-190 http://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201302014.htm [8] 张力, 张祖勋, 张剑清. Wallis滤波在影像匹配中的应用[J].武汉测绘科技大学学报, 1999, 24(1):24-27 http://www.cnki.com.cn/Article/CJFDTOTAL-WHCH199901005.htm Zhang Li, Zhang Zuxun, Zhang Jianqing.The Image Matching Based on Wallis Filtering[J].Journal of Wuhan Technical University of Surveying and Mapping, 1999, 24(1):24-27 http://www.cnki.com.cn/Article/CJFDTOTAL-WHCH199901005.htm [9] Al-Rousan N, Cheng P, Petrie G, et al. Automated DEM Extraction and Orthoimage Generation from SPOT Level 1B Imagery[J]. Photogrammetric Engineering and Remote Sensing, 1997, 63(8):965-974 https://www.researchgate.net/publication/290516781_Automated_DEM_extraction_and_orthoimage_generation_from_SPOT_Level_1B_imagery [10] Kim T. A Study on the Epipolarity of Linear Pushbroom Images[J]. Photogrammetric Engineering and Remote Sensing, 2000, 66(8):961-966 https://www.researchgate.net/publication/283431921_A_Study_on_the_Epipolarity_of_Linear_Pushbroom_Images [11] 胡芬, 王密, 李德仁, 等.基于投影基准面的线阵推扫式卫星立体影像对近似核线影像生成方法[J].测绘学报, 2009, 38(5):428-436 http://www.cnki.com.cn/Article/CJFDTOTAL-CHXB200905012.htm Hu Fen, Wang Mi, Li Deren, et al. Generation of Approximate Epipolar Images from Linar Pushbroom Satellite StereoImagery Based on Projection Reference Plane[J]. Acta Geodaetica et Cartographica Sinica, 2009, 38(5):428-436 http://www.cnki.com.cn/Article/CJFDTOTAL-CHXB200905012.htm [12] 卢华兴. DEM误差模型研究[D]. 南京: 南京师范大学, 2008 http://d.wanfangdata.com.cn/Thesis/Y1330990 Lu Huaxing. Research on DEM Error Model[D]. Nanjing:Nanjing Normal University, 2008 http://d.wanfangdata.com.cn/Thesis/Y1330990 [13] Powell M J D. The Theory of Radial Basis Function Approximation in 1990[M]. Cambridge:University of Cambridge Press, 1990 -
点击查看大图
图(3) / 表(2)
计量
- 文章访问数: 1421
- HTML全文浏览量: 139
- PDF下载量: 510
- 被引次数: 0
