一种改进的高分辨率卫星影像数字表面模型生成方法

刘世杰，晏飞，3，王卫安，李荣兴

（1.同济大学 空间信息科学及可持续发展应用中心，上海200092；2.同济大学 测绘与地理信息学院，上海200092；3.重庆市地理信息中心，重庆401121）

卫星遥感影像立体匹配是航天摄影测量中的关键技术和核心问题［1］，许多学者提出了各种匹配方法，包括最小二乘匹配［2］、归一化互相关（normalized cross correlation，NCC）［3］等灰度区域匹配、全局／半全局匹配［4］以及相位相关等频域匹配算法［5］.各种匹配方法差异反映了影像立体匹配的复杂性和高难度.在图像获取的过程中，三维空间投影到二维的平面图像上，导致许多信息丢失，因此图像匹配属于逆向工程类问题，这是一个病态问题［6］.除了因遮挡造成匹配失败，还有可能由于重复纹理导致误匹配，以及由于图像存在噪声和几何变形使得匹配算法不稳定等问题.针对这些问题，有学者提出和研究了基于尺度不变特征转换（SIFT）的立体匹配算法［7］，以克服立体影像尺度差异和几何变形造成的匹配困难，但由于特征比较稀疏且分布容易不均，不利于生成较高分辨率的数字表面模型（digital surface model，DSM）；针对DSM生成，有学者研究了基于物方的匹配方法如铅垂线轨迹法（VLL）［8］，该方法以地面为采样计算单元，通过高程值的改变计算像方最大相关系数，从而在成功匹配的同时获取地面DSM，该方法速度快，但仍面临影像畸变和遮挡带来的问题.

本文综合运用SIFT算子和NCC算子，并考虑核线和视差约束，提出了一种改进的基于三角网视差传递约束的立体影像分层匹配方法，并利用上海崇明区域的WorldView-1影像和浙江舟山区域的资源三号卫星影像进行实验，验证了本文方法的有效性，获得了高精度的DSM.

1 影像成像模型

有理函数模型（rational function model，RFM）是目前高分辨率卫星广泛采用的成像模型，卫星影像供应商直接提供有理函数模型参数.有理函数模型是将像点坐标表示为以相应地面点空间坐标为自变量的多项式的比值，如式（1）所示［9］.

式中：（rn，cn）是像平面坐标；（Xn，Yn，Zn）是地面坐标；它们均是经过标准化后的坐标值，以减小计算过程中由于数量级差别过大引入的误差；对于标准有理函数模型，pi（i=1，2，3，4）是一个20项三阶多项式.

使用原始影像提供的有理函数模型参数（RPCs）直接进行定位具有明显的系统误差，需利用地面控制点进行误差补偿，提高定位精度.一般说有两种补偿方案［10］，一种是在物方直接对定位结果进行误差纠正，另一种是先纠正像平面坐标的误差，再进行立体定位.后一种方案理论上严密，对于高分辨率卫星影像，一般仿射变换模型能满足精度要求，且较稳健，故本文采用像方仿射变换误差补偿模型.

2 影像立体匹配和DSM生成

本文提出的影像立体匹配生成DSM的流程如图1所示，包括影像预处理、金字塔逐层约束匹配获取密集匹配点、立体定位和点云内插生成DSM等.

2.1 影像预处理

预处理包括影像去噪、增强和影像金字塔构建.数字影像获取不可避免会存在噪声，因此需要进行去噪处理.对于图像平滑去噪算子而言，不仅需要过滤掉噪声，更重要的是要保持边缘信息.本文采用Saint-Marc提出的自适应平滑算子［11］，属于非线性滤波，通过一个迭代的过程，达到去除噪声并保留边缘的目的.影像增强采用 Wallis滤波［12］，用于增强纹理并提高信噪比，便于特征提取和匹配.在建立金字塔影像过程中，先利用高斯模糊算子对当前层影像进行低通滤波，然后进行降采样，得到上一层影像，最终建立包含4层的影像金字塔.

图1 高分辨率立体影像生成DSM方法流程图Fig.1 Workflow of DSM generation from high resolution satellite stereo images

2.2 影像匹配

2.2.1 SIFT初匹配

影像初始匹配点将作为下层加密匹配的约束条件，对准确性的要求非常高.SIFT特征具有对尺度、旋转、亮度变化保持不变的稳定特性，故采用SIFT算子进行初匹配，结合RANSAC（随机抽样一致）检验算法［13］，检验并剔除可能的误匹配，保证初匹配的准确性.

2.2.2 特征点和格网点匹配

本文采用Förstner特征算子［14］，它能给出特征点的类型且精度较高.其基本思想是：对于角点，对最佳窗口内通过每个像元的边缘直线（垂直于梯度方向）进行加权中心化，得到角点的坐标；对于圆状点，对最佳窗口内通过每个像元的梯度直线进行加权中心化，得到圆心的坐标.

利用特征算子提取的特征点主要分布在影像灰度存在变化具有明显特征的地方，其分布不均匀且数量有限.为了生成DSM，需要更多密集的匹配点.因此，除特征点外，在影像上按一定间隔划分格网，并去除灰度方差低于一定阈值的格网点，然后对保留的格网点逐个在待匹配影像上进行搜索匹配.

特征点和格网点匹配采用NCC算法.在传统的匹配方法中，一般采用矩形窗口计算相关系数，这种方法适用于几何变形较小的影像匹配，而对于几何变形较大的情况，矩形窗口不再适用.因此，本文提出一种基于视差几何纠正的相关系数计算方法.在该方法中，对于待匹配点，根据附近的视差三角形格网点上的视差计算几何变换系数，纠正待匹配窗口并重采样，使左右影像匹配窗口内容一致.另外，由于采用固定模板进行匹配会降低匹配成功率，本文采用一种自动改变匹配窗口大小和形状的匹配策略，它通过把匹配窗口向四个方向扩展，直到相关系数不再增大为止，从而提高匹配的准确性和成功率.2.2.3 特征线匹配

特征线是影像和地面上的重要特征，反映了影像灰度或地面高程的不连续性.通过对特征线的提取和匹配，可以实现对地面高程不连续的精确定位，对精确DSM的构建具有重要意义.对于特征线段的提取，本文采用Akinlar提出的EDLines算法［15］.该方法无需对参数进行调整而可得到精确的结果，且速度快.该检测算子利用亥姆霍兹规则进行线段验证，剔除错误的线段.

对线段的匹配，除几何位置、长度和方向（与核线夹角）外，还考虑线段周围的灰度信息，采用平行于线段的矩形滑动窗口计算灰度相关系数.在匹配搜索过程中，综合利用已匹配点所形成的三角网约束和核线约束来确定候选匹配线段.

2.2.4 匹配搜索约束

（1）三角网约束

对于滩涂、岛礁等特征不明显、纹理重复或缺乏的区域，直接进行影像匹配很容易造成误匹配，因此本文采用三角网约束的策略进行匹配［16］.点的匹配首先在最上层低分辨率的影像上进行，匹配成功的点传递到较高分辨率的下一层重新进行匹配以获得更高的匹配精度，并把匹配上的点建立一个视差表面的狄洛尼三角网，用作加密匹配点的搜索约束.这个过程不断重复，直到最下层原始分辨率影像.在每个影像层上，同名点搜索区域根据上一层的视差确定.从影像最高层到最低层，匹配上的同名点越来越多，获得的减小搜索范围的视差信息就越精确.

（2）核线约束

线阵CCD（电荷耦合器件）推扫式影像遵循一种“动态”的行中心投影成像方式，因此无法像框幅式中心投影影像那样基于成像几何关系给出严格直观的核线定义.目前，在卫星影像近似核线几何关系的各种描述中，基于投影轨迹法的核线定义建立在成像的几何约束条件之上，在理论上最为严密.线阵CCD推扫式影像核线类似于双曲线，但在影像范围内可以近视看作直线［17］.

根据三角网视差约束条件，得到左影像上某个待匹配点在右影像上的估计位置.然后以估计的位置为中心，建立矩形搜索窗口.再由左影像上的待匹配点计算在右影像上的核线.匹配搜索的时候，就只需在这个搜索窗口内的核线段上搜索即可.由于计算的核线会存在一定的误差，因此实际情况下会把核线左右两侧扩宽一定像素，然后沿着这条具有一定宽度的核线段搜索.这样不仅可以大大提高搜索效率，而且有效地降低了匹配的不确定性.

2.2.5 误匹配剔除

（1）双向匹配

双向匹配算法的基本思想是在第一次匹配时，以左影像某个特征点或格网点作为待匹配点，在右影像上一定视差搜索范围内计算所有候选点的相关系数，取相关系数最大值点为有效匹配点，即进行所谓的正匹配.然后将上一次匹配的结果（右影像上的点）作为待匹配点，在左影像上一定视差范围内进行搜索，以同样的方法进行二次匹配，即所谓的逆匹配.如果两次结果的差异小于设定的阈值（如1个像素），则接受该匹配结果，否则舍弃，达到剔除误匹配点的目的.

（2）视差平滑约束

尽管在匹配过程中采用了三角网视差约束和双向匹配的算法，仍然不能保证所有的匹配点都是正确的.在本文的匹配算法中，引入了视差平滑约束条件.对于每一对同名点，使用它们的邻域内的同名点拟合一个视差平面（不跨越边线）.如果该点视差与视差平面的距离超过阈值（3倍中误差），认为该点是误匹配并剔除.实验表明，该方法可有效剔除绝大多数错误匹配.

2.3 DSM构建

成功匹配的同名点通过成像模型立体定位得到对应地面点的三维坐标，但影像匹配定位得到的三维点的排列是不规则的，因此需要通过内插获得规则格网的DSM.综合考虑时间复杂度和插值精度，本文采用反距离加权插值方法构建DSM.

3 实验与分析

3.1 崇明WorldView-1卫星影像生成DSM

（1）实验数据

该实验使用上海市崇明岛东滩区域的WorldView-1影像，分辨率0.5m，为不同时期（2012／01／25 和 2012／04／06）异 轨 立 体，交 会 角28.7°.实验区域大小约为18km×18km，见图2a和图2b，区域内地势平坦，主要为农田和滩涂；该区域位于相邻两景影像内（供应商按一定分幅规格进行分幅，且相邻影像具有一定重叠度），两期共四幅影像.为了纠正影像几何定位误差，使用了17个地面控制点，其地面坐标由GPS（全球定位系统）测得；同时，在重叠区域选取了9个连接点参与平差解算，以保证上下两幅影像生成的DSM高程连续.

图2 崇明WorldView-1影像及控制点与连接点分布（三角点：控制点，圆点：连接点）和生成的DSMFig.2 Chongming WorldView-1image and distribution of control points and tie points （triangle：control point，circle：tie point）and generated DSM

（2）DSM 和精度分析

首先利用地面控制点和影像连接点纠正影像成像模型误差（采用像方仿射模型），平差后成像模型误差在0.5像素以内，最后生成的拼接后的DSM如图2c所示.从生成的DSM上可清晰看出地面纹理，特别是线性特征，如道路和不同地表覆盖的边线等，反映了线特征提取和匹配的作用.利用该区域内均匀分布的100个高程注记点作为检核来评价DSM的精度，计算得DSM的高程中误差约为0.7m.

3.2 舟山资源三号卫星影像生成DSM

（1）实验数据

实验影像数据为中国资源三号卫星2013年1月24日拍摄的全色前、后视立体影像，前、后视夹角为44°，地面分辨率为3.6m，区域大小约为58km×58km，如图3a所示.该区域包含众多岛屿，高程变化在海拔550m范围内.

（2）DSM和精度分析

由于没有该区域的地面控制点，故直接采用影像附带有理函数模型进行地面定位，生成的DSM如图3b所示.为了检验DSM的精度，利用该区的ASTER GDEM（先进星载热辐射反射探测仪全球数字高程模型）作为参考，先将生成的DSM与ASTER GDEM进行配准，然后比较并统计高程差异.考虑到舟山群岛区域地形复杂，按山地（高程大于20m的山地区域）、城区（城镇区域，主要分布在沿海）和平地（高程小于20m的沿海空旷区域）三类地形检验DSM精度.结果显示，山地区域DSM与ASTER GDEM的偏差均方根为15.5m，城区为8.9m，平地为4.6m，两者吻合较好.该差异一部分是因为ASTER GDEM 本身存在误差，另一部分是因为ASTER和资源三号影像拍摄相隔10余年（ASTER GDEM是利用1999发射的ASTER传感器所获取的15m分辨率影像生产得到，资源三号影像为2013年拍摄），由于人类活动而造成了地面变化（如山体开挖、城区建设等）.

4 结语

本文综合运用SIFT算子的可靠性和灰度相关算子的密集匹配特性，并考虑核线和视差约束，提出了一种改进的基于三角网视差传递约束的卫星影像金字塔分层匹配方法，实现特征点、格网点和特征线的匹配，用于从高分辨率卫星立体影像生成DSM.该方法在金字塔分层匹配中引入三角网视差传递约束和核线约束以提高匹配效率和准确性，同时采用自适应窗口改进来提高匹配成功率，对于误匹配点通过双向匹配和视差平滑约束进行剔除，最后将匹配成功的点经立体几何模型定位生成点云，并内插获得DSM.由于采用了自适应窗口和多种约束及误匹配检测算法，该DSM生成方法具有较好的鲁棒性.实验证明，该方法能提高匹配的准确性和成功率，获得高精度的DSM.

图3 舟山资源三号影像和生成的DSMFig.3 Zhoushan ZY-3image and generated DSM

［1］ Gruen A.Development and status of image matching in photogrammetry［J］.The Photogrammetric Record，2012，27（137）：36.

［2］ Zhang L.Automatic digital surface model（DSM）generation from linear array images［D］，Zurich：Swiss Federal Institute of Technology Zurich，2005.

［3］ Tsai D M，Lin C T，Chen J F.The evaluation of normalized cross correlations for defect detection［J］.Pattern Recognition Letters，2003，24：2525.

［4］ Alobeid A，Jacobsen K，Heipke C.Comparison of matching algorithms for DSM generation in urban areas from ikonos imagery［J］.Photogrammetric Engineering &Remote Sensing，2010，76（9）：1041.

［5］ Keller Y，Averbuch A.A projection-based extension to phase correlation image alignment［J］.Signal Processing，2007，87（1）：124.

［6］ Terzopoulos D.Regularization of inverse visual problems involving discontinuities［J］.Pattern Analysis and Machine Intelligence，IEEE Transactions on，1986，4：413.

［7］ Lowe D G.Distinctive image features from scale-invariant keypoints［J］.International Journal of Computer Vision，2004，60（2）：91.

［8］ 郑顺义，张祖勋，张剑清.基于物方影像匹配和概率松弛的断面自动提取［J］.测绘信息工程，2004，29（2）：26.ZHENG Shunyi，ZHANG Zuxun，ZHANG Jianqing.Automatic profile extraction based on object space image matching and probability relaxation algorithm［J］.Journal of Geomatics，2004，29（2）：26.

［9］ Madani M.Real-time sensor-independent positioning by rational functions［C］／／Proceedings of ISPRS Workshop on Direct Versus the HRSC data.Indirect Methods of Sensor Orientation.Barcelona：［s.n.］，1999：64-75.

［10］ 刘世杰.高分辨率卫星遥感影像成像模型与定位技术研究［D］.上海：同济大学，2008 LIU Shijie.Imaging models and geo-positioning technology for high resolution satellite imagery ［D］.Shanghai： Tongji University，2008.

［11］ Saint-Marc P，Chen J S，Medioni G.Adaptive smoothing：a general tool for early vision［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1991，13（6）：514.

［12］ 张力，张祖勋，张剑清.Wallis滤波在影像匹配中的应用［J］.武汉测绘科技大学学报，1999（1）：24.ZHANG Li，ZHANG Zuxun，ZHANG Jianqing.The image matching based on wallis filtering［J］.Journal of Wuhan Technical University of Surveying and Mapping，1999（1）：24.

［13］ Fischler M A，Bolles R C.Random sample consensus：a paradigm for model fitting with applications to image analysis and automated cartography［J］.Communications of the ACM，1981，24（6）：381.

［14］ Forstner W，Gulch E.A fast operator for detection and precise location of distinct points，corners and centres of circular features ［C ］／／Proceedings of ISPRS intercommission conference on fast processing of photogrammetric data，Interlaken：［s.n.］，1987：281-305.

［15］ Akinlar C，Topal C.EDLines：a real-time line segment detector with a false detection control［J］.Pattern Recognition Letters，2011，32（13）：1633.

［16］ Hwangbo J W.Integration of orbital and ground imagery for automation of rover localization ［D］.Ohio：Ohio State University，2010.

［17］ 胡芬，王密，李德仁，等.利用投影轨迹的卫星影像近似核线重排快速算法［J］.武汉大学学报：信息科学版，2009，34（12）：1431.HU Fen，WANG Mi，LI Deren，etal.Approximate epipolar rearrangement algorithm of satellite stereo-imagery by projection trajectory［J］.Geomatics and Information Science of Wuhan University，2009，34（12）：1431.