留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

嫦娥二号影像密集匹配及DEM制作方法

王奋飞 任鑫 刘建军 施俊 李春来

王奋飞, 任鑫, 刘建军, 施俊, 李春来. 嫦娥二号影像密集匹配及DEM制作方法[J]. 武汉大学学报 ● 信息科学版, 2017, 42(7): 918-923. doi: 10.13203/j.whugis20150147
引用本文: 王奋飞, 任鑫, 刘建军, 施俊, 李春来. 嫦娥二号影像密集匹配及DEM制作方法[J]. 武汉大学学报 ● 信息科学版, 2017, 42(7): 918-923. doi: 10.13203/j.whugis20150147
WANG Fenfei, REN Xin, LIU Jianjun, SHI Jun, LI Chunlai. The Methods of CE-2 Image Dense Matching and Lunar DEM Extraction[J]. Geomatics and Information Science of Wuhan University, 2017, 42(7): 918-923. doi: 10.13203/j.whugis20150147
Citation: WANG Fenfei, REN Xin, LIU Jianjun, SHI Jun, LI Chunlai. The Methods of CE-2 Image Dense Matching and Lunar DEM Extraction[J]. Geomatics and Information Science of Wuhan University, 2017, 42(7): 918-923. doi: 10.13203/j.whugis20150147

嫦娥二号影像密集匹配及DEM制作方法

doi: 10.13203/j.whugis20150147
基金项目: 

国家自然科学基金 11273037

详细信息
    作者简介:

    王奋飞, 助理研究员, 主要从事月球遥感图像处理算法研究。wangff@bao.ac.cn

    通讯作者: 任鑫, 副研究员。renx@nao.cas.cn
  • 中图分类号: TP751;P208

The Methods of CE-2 Image Dense Matching and Lunar DEM Extraction

Funds: 

The National Natural Science Foundation of China 11273037

More Information
    Author Bio:

    WANG Fenfei, assistant researcher, specializes in the lunar surveying and mapping methods. E-mail:wangff@bao.ac.cn

    Corresponding author: Ren Xin, associate researcher. E-mail: renx@nao.cas.cn
图(3) / 表(2)
计量
  • 文章访问数:  1421
  • HTML全文浏览量:  139
  • PDF下载量:  510
  • 被引次数: 0
出版历程
  • 收稿日期:  2015-05-12
  • 刊出日期:  2017-07-05

嫦娥二号影像密集匹配及DEM制作方法

doi: 10.13203/j.whugis20150147
    基金项目:

    国家自然科学基金 11273037

    作者简介:

    王奋飞, 助理研究员, 主要从事月球遥感图像处理算法研究。wangff@bao.ac.cn

    通讯作者: 任鑫, 副研究员。renx@nao.cas.cn
  • 中图分类号: TP751;P208

摘要: 嫦娥二号卫星搭载了一台两线阵CCD立体相机,以线阵推扫成像方式获取了覆盖全月球的7 m分辨率影像。与嫦娥一号影像相比,影像分辨率大大提高,月表地貌细节更加丰富,制作高分辨率月表DEM更加复杂而困难,现有的基于嫦娥影像制作方法难以满足应用需求。针对嫦娥二号影像特点,在基于RPC模型生成近似核线影像基础上,进行金字塔影像匹配,实现了原始影像的逐像素匹配;计算匹配点的月面坐标后,采用基于反距离插值方法生成月表三维地形;对于图像阴影等地形漏洞,采用基于径向基函数进行插值填补。试验结果表明,该方法生成的月表三维地形DEM数据,能准确真实地表达高分辨率的月表地形地貌细节,并较好地应用于嫦娥二号全月球三维地形DEM数据产品的生产任务中。

English Abstract

王奋飞, 任鑫, 刘建军, 施俊, 李春来. 嫦娥二号影像密集匹配及DEM制作方法[J]. 武汉大学学报 ● 信息科学版, 2017, 42(7): 918-923. doi: 10.13203/j.whugis20150147
引用本文: 王奋飞, 任鑫, 刘建军, 施俊, 李春来. 嫦娥二号影像密集匹配及DEM制作方法[J]. 武汉大学学报 ● 信息科学版, 2017, 42(7): 918-923. doi: 10.13203/j.whugis20150147
WANG Fenfei, REN Xin, LIU Jianjun, SHI Jun, LI Chunlai. The Methods of CE-2 Image Dense Matching and Lunar DEM Extraction[J]. Geomatics and Information Science of Wuhan University, 2017, 42(7): 918-923. doi: 10.13203/j.whugis20150147
Citation: WANG Fenfei, REN Xin, LIU Jianjun, SHI Jun, LI Chunlai. The Methods of CE-2 Image Dense Matching and Lunar DEM Extraction[J]. Geomatics and Information Science of Wuhan University, 2017, 42(7): 918-923. doi: 10.13203/j.whugis20150147
  • 嫦娥二号(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)

      式中,ff′代表原始影像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)

      式中,φ(‖xxi‖)是径向基函数;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
      合计
      /s
      035 3086 144×39 76014090230
      059 5126 144×33 98013585220

      本文对最终匹配结果的灰度相关性进行了统计。一般而言,在图像匹配处理中,生成用于相对定向的连接点时,对点的精度要求较高,灰度相关系数往往大于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 30898.696.389.371.832.3
      059 51277.364.647.929.510.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

参考文献 (13)

目录

    /

    返回文章
    返回