1 引 言

地形高程信息的提取及数字高程模型(digital elevation model，DEM)的生成是合成孔径雷达(synthetic aperture radar，SAR)最重要的应用之一。雷达摄影测量(StereoSAR，又称Radargrammetry)是雷达遥感提取DEM的一种重要方法，它利用SAR影像的幅度信息和立体像对的成像几何关系来提取地面目标的三维信息^[1]。StereoSAR技术对立体像对获取的时间、基线和地表覆盖类型等条件的要求并不苛刻，数据获取相对容易，提取的DEM可靠性较高，尤其是随着新一代高分辨率星载SAR系统的投入运行，StereoSAR技术提取DEM的精度逐渐提高，具备了广阔的应用前景。

经过20余年的发展，雷达摄影测量的基础理论已趋于成熟^{[2, 3]}。各国学者相继开展了多项利用StereoSAR技术提取DEM的应用研究。文献^[4]研究开发了机载立体SAR解析测图系统SARMS，并将其应用于1∶5万以下中小比例尺地形图的生产和修测。文献^[5]利用RADARSAT立体像对在加拿大西部进行了DEM提取试验，高程精度约为25 m。文献^[6]在华南高山区利用RADARSAT数据进行了DEM提取试验，在使用9个控制点的情况下得到的DEM高程精度约为35 m。文献^[7]开展了利用地面控制点对星载SAR与光学影像联立构建广义立体像对进行复合式定位的研究，取得了10 m左右的定位结果精度。随着SAR影像的分辨率越来越高，StereoSAR提取DEM的精度也逐渐提高。文献^[8]利用聚束模式TerraSAR-X 立体像对和20多个控制点获取到的DEM高程精度可达米级。文献^[9]采用聚束模式COSMO-SkyMed数据得到的DEM高程精度优于5 m。本文采用了条带模式TerraSAR-X立体像对，在无地面控制点的情况下，也取得了可观的结果。

SAR立体像对的匹配是StereoSAR提取DEM处理过程中的关键问题。目前SAR影像匹配主要有两类方法：基于灰度和基于特征。基于特征的匹配算法是通过计算特征间的相似度来确定影像间的同名特征，从而拟合建立影像间的几何映射关系，其匹配过程通常具有稀疏性，不适用于逐像素的匹配，因此在StereoSAR数据处理中应用较少。基于灰度的匹配方式虽然在影像信息贫乏的区域成功率不高，但若能采用适当的匹配策略以及误差控制方法，也能取得较好的结果。鉴于此，本文采用了基于灰度的匹配方法。

SAR特有的侧视斜距成像方式及由此导致的透视收缩、叠掩和阴影等现象，造成SAR影像在地形起伏区域普遍存在着几何畸变^[10]。由于成像视角不同，SAR立体像对中的几何畸变会呈现显著的差异，从而使得立体像对之间的匹配非常困难。因此，如何利用尽可能多的信息来提高匹配精度是一项值得考虑的事情。文献^[11]提出的利用像方和物方信息的多影像匹配方法采用金字塔分层匹配策略，在每一层金字塔上匹配时，通过物方信息剔除部分误匹配点，用于下一层金字塔影像的匹配，提高了匹配成功率。本研究同样采用了金字塔匹配策略，利用已有粗分辨率的DEM作为辅助信息，控制金字塔分层匹配中误匹配点的传递，目的是提高SAR立体像对匹配的有效性和精度，从而提高结果DEM的高程精度和可靠性。

StereoSAR是利用从不同的入射角分别观测获取的两幅具有一定重叠度的SAR影像，通过同名点的影像坐标提取对应地物的高程信息。SAR立体像对的获取可以采用同侧或者是异侧立体观测两种方式。由于异侧观测获取的雷达立体像对的色调和几何畸变差异非常大，难以实现准确可靠的匹配，因此试验中一般采用同侧观测获取的立体像对^{[12 ,13]}。

StereoSAR提取高程信息的主要依据是视差效应。雷达影像上存在由高差引起的投影差，不同的视角将造成同一地物有不同的投影差，从而在立体像对上形成左右视差。任何一个地物目标的视差大小与它的高程值有关，当成像几何不变时，同一水平位置上，高程高的点的视差比高程低的点的视差大^[14]。图 1简化表示了同侧立体观测获取的地距影像对中高程与视差之间的关系。

图 1 立体视差与高程的关系示意图 Fig. 1 Sketch map of the relationship between stereo parallax and elevation

图选项

假设O₁和O₂是构成立体像对的两次观测成像时SAR传感器所处的位置，AB为地面目标，且与O₁O₂共面。雷达量测天线至目标的斜距，对于星载SAR系统来说，斜距相对于目标尺寸足够大，因此可以近似认为其投影线与视线垂直。地面目标AB的顶端B在构成立体像对的两幅SAR影像上的投影分别为B₁和B₂。假设目标的底端A在地面上，其对应视差为0，则在立体像对上，B点对应的视差为p₂-p₁。

由三角几何原理可知，B点上的视差与高程具有如下关系^[14]

点击浏览公式

该高程是B点相对于A点的高程。式中,分母项对于同侧成像采用减法，对于异侧成像采用加法。θ₁和θ₂为卫星侧视角;m为影像的比例尺系数，一般取为影像列方向采样间隔。

需要注意的是，图 1中传感器在O₁和O₂处的飞行方向相互平行，主辅影像在方位向上已对齐，即无上下视差，高程仅与距离向上的左右视差有关。上下视差是指同名点在方位向上的坐标差。由于卫星轨道较飞机的飞行轨迹要稳定得多，星载SAR影像在方位向上变形较小，通过线性变换矫正即可将立体SAR影像对在方位向上配准，从而消除上下视差。本文的重点在于研究SAR立体像对的匹配，因此对于如何消除立体像对的上下视差不作详细阐述，采用的试验数据是已经消除了上下视差的SAR立体像对。手工量测选取若干同名点对进行检查，上下视差均控制在两个像元以内，通过合理设置匹配窗口的滑动范围，可以消除上下视差对匹配结果的影响。此外，当采用金字塔影像匹配策略时，在粗分辨率的金字塔分层上同名点几乎不存在上下视差。

由于斑点噪声以及不同视角造成的几何辐射畸变差异，搜索同名点是立体SAR数据处理流程中最为困难也最耗时的一步。为了提高匹配的效率，同时保证匹配结果精度，本文采用了由粗到精的金字塔分层匹配策略，先从顶层金字塔搜索同名点，将其匹配结果传递给下一层，逐渐精化匹配过程^[16]。金字塔的建立采用了经典的移动平均法，下一层金字塔上相邻4个像元通过灰度平均生成上一层金字塔上的一个对应像元。金字塔的层数根据估计的最大左右视差来决定，即需要保证最大左右视差在顶层金字塔上能大于一个像元。搜索时，从顶层金字塔出发，通过匹配得到主影像上每个像元的左右视差。在下一层金字塔匹配时，主影像上的每一个像元将会有一个初始的左右视差，取值为上一层金字塔上对应像元的左右视差的两倍。根据初始左右视差可以定位到该像元在辅影像上的同名点的初始位置，搜索时就只需要以该位置为中心，上下左右搜索，从而大大减小搜索范围。

SAR影像的匹配有多种匹配准则，如归一化相关系数法、Tupin提出的基于比值的方法、互信息法等^{[16, 17, 18]}。这些匹配准则在理论上来说很难分出优劣，都是基于灰度的统计来进行匹配。相关系数法是摄影测量学中影像匹配的经典方法，可以克服线性畸变的影响，是一种比较稳健成熟的方法。互信息法的计算量较相关系数法大很多，Tupin提出的基于比值的方法尽管在某些试验中取得了更好的效果，但应用在本研究所采用的立体像对中，匹配结果的精度并没有得到明显提高。在本文的试验中，相关系数法取得了最佳的结果，因此本试验采用了基于灰度的归一化相关系数匹配算法，其匹配测度计算如下^[19]

点击浏览公式

虽然金字塔分层匹配策略中低层金字塔影像能够从高层金字塔匹配结果中获得同名点视差初值，然而，顶层金字塔却没有良好的同名点视差初值。顶层金字塔的匹配结果是下层金字塔匹配的基础，若粗差较多，将对后续的匹配造成不利影响，因此保证顶层金字塔的匹配正确率是非常重要的。此外，高层金字塔匹配的结果向低层金字塔传递的过程中也不可避免地传递了误差，因为高层金字塔的匹配结果不正确，将造成下一层金字塔选取初始同名点的错误，从而影响搜索范围。若不对这些误差加以控制，将对最终DEM的精度造成极大的影响。因此，有必要对每一层金字塔匹配结果进行误差控制。在很多区域，由于低空间分辨率的DEM数据已经存在，能够用来对金字塔匹配结果进行误差控制，可作为StereoSAR技术的有效辅助数据。

2000年2月，美国、德国及意大利合作开展的“航天飞机雷达地形测量任务”(SRTM)计划，利用航天飞机搭载双天线SAR系统获取的数据进行雷达干涉测量处理，获取了全球大部分陆地的DEM。3弧秒分辨率(水平分辨率约90 m)的SRTM DEM已被公布在互联网上供全球的科学家免费使用^[20]。SRTM DEM的标称高程精度达到了20 m，是目前能够免费获取的质量最佳的全球覆盖DEM数据，因此本研究中选用其作为外部粗分辨率DEM来辅助SAR立体像对的匹配。采用SRTM DEM辅助SAR立体像对匹配的基本流程如图 2所示。

图 2 StereoSAR匹配流程图 Fig. 2 Flow chart of StereoSAR image matching

图选项

首先对已经消除方位向视差的立体像对生成金字塔影像，将外部DEM通过雷达编码转化到雷达影像坐标系下，并重采样到每层金字塔影像的分辨率下，根据式(1)中视差与高程的关系将高程转换成每个像元上的初始视差。考虑到传统的基于距离-多普勒定位模型的雷达编码过程需要耗费大量的计算时间，本文采用了有理函数模型来替代距离-多普勒模型用于雷达编码，以有效提高雷达编码的处理效率^[20]。

在立体像对匹配的过程中外部DEM的作用主要体现在以下两方面：

首先，在顶层金字塔同名点搜索的过程中，外部DEM可以提供匹配初值。距离向视差的初值可以根据外部DEM高程转化得到的视差来确定，搜索范围以初值为中心，左右各搜索若干像元。方位向可以只搜索一行，也可以上下各多搜索一行，根据方位向视差消除的效果而定。

其次，在从金字塔当前层向下一层传递的过程中，可以通过已有DEM转到雷达坐标系下的高程点阵图和当前金字塔上的高程点阵图比较从而判断粗差点，避免粗差点向下一层金字塔传递。试验中主要有两种像元被认为是粗差点，一种是相关系数低于一定阈值的像元，另一种是由匹配视差计算出高程值后，该高程与已有DEM高程之差大于一定阈值的像元。高程差的阈值选取根据切比雪夫法则来设定^[20]。根据切比雪夫法则，误差分布在μ±4δ间隔内的概率至少为94%，μ和δ分别为样本均值和标准差。因此，试验中设SRTM DEM高程误差的4倍为阈值，仅选择高程差小于阈值的那些点来传递到下一层金字塔。而粗差点的视差将被从SRTM DEM计算出来的视差替代，从而避免了粗差的累计传递。

本试验使用的是德国TerraSAR-X卫星升轨获取的覆盖马来西亚瓜拉江沙地区的一对同侧立体像对，分别成像于2008-09-18和2008-09-13，侧视角分别为42.7°和21.3°。两景均为条带模式地距影像，分辨率约为3 m。该区域海拔高程变化范围约为0～1600 m，沿河流两岸地势较为平坦，而在河流西北和东南方向均为覆盖有茂密植被的高山区。立体像对幅度影像如图 3所示，可以看出，由于成像视角相差了一倍，两幅影像中地形起伏巨大的山区在几何畸变上存在显著的差异，这就使得立体像对的精确匹配十分困难。因此，本试验中采用90 m分辨率的SRTM DEM作为外部DEM来辅助这对TerraSAR-X立体像对的匹配。

图 3 试验所用TerraSAR-X立体像对( DLR/Infoterra) Fig. 3 TerraSAR-X stereo data pair as the experimental data( DLR/Infoterra)

图选项

试验中，根据目视判读所得最大视差计算出该立体像对合适的金字塔层数为6层。影像金字塔由上到下，自顶层起分辨率依次倍增，即每个上层影像像元对应于4个下层影像像元。考虑到各层影像在分辨率和信噪比方面具有的明显差异，各层影像匹配采用的参数配置有所差别，具体参数配置如表 1所示。由于高层影像是由低层影像经低通滤波构建得到，因此高层影像信噪比更高，故而设置了更高的相关系数阈值。

搜索是影像匹配中一个非常耗时的过程，所以搜索区域应该尽可能地小，以降低搜索时间。搜索区域大，匹配时间就长，而且搜索区域过大，会增加错误匹配的概率。由于匹配之前对主从影像消除了方位向视差，因此在低分辨率的高层金字塔影像上只在距离向上搜索，而在较高分辨率的低层金字塔影像上，在方位向上下各自多搜索一个像元。距离向上的搜索范围应根据最大视差来设置，顶层金字塔影像匹配时距离向搜索范围与最大视差相当。

为了对常规的和采用外部DEM辅助的两种立体SAR影像匹配方法进行对比，本文分别采用这两种方法对该立体像对进行了处理。在这两种方法中，金字塔各层级对应的高程点阵图上均只保留相关系数不小于所设定阈值的像元。采用外部DEM辅助的方法中对于相关系数低于阈值的像元则采用由外部DEM估计得到的偏移量予以替换。

图 4(a)是采用常规匹配方法得到的主影像坐标系下的高程点阵图。从中可以看出，由于粗差点的影响，许多地形细节变得模糊不清，且存在较多的错误点。尤其是在图像边缘地区和地形起伏较大的地区，高程误差往往高达数百米，致使图中左下角存在一个明显的空洞。这是因为左下角的山顶只在参考影像上出现，而在待匹配影像上不存在，因此该山顶不应出现在高程点阵图上，然而由于缺少先验信息，在匹配时仍然会进行搜索，找出相关系数最大的点，认为是同名点。实际上，搜索出来的相关系数最大点并不是真实的同名点。与SRTM的高程点阵图比较发现，该山区处的高程错误高达900 m。

图 4 雷达坐标系下的高程点阵图(单位为m) Fig. 4 Elevation maps in the SAR coordinate system (unit: m)

图选项

图 4(b)为经过外部DEM辅助匹配后生成的影像坐标系下的高程点阵图。可以看出，该高程点阵图的目视效果有了明显改善，地形起伏较大的地区纹理也非常清晰。该高程点阵图中左下角的山顶并没有匹配出来，这得益于外部DEM为匹配提供了初值，由于发现该初始同名点超出了辅影像的范围，因此不进行搜索，从而避免了错误的产生。

由于左下角山顶处为明显的错误值，而由高程点阵图生成DEM的过程中需要采用SAR定位模型对高程点阵图进行地理编码，若这些错误的高程点参与到地理编码的过程中，将不仅造成结果DEM的高程错误，而且会造成平面位置的显著偏移，从而降低结果DEM的质量。所以在试验中设定边缘部分(图 4中方框外部分)不参与地理编码。图 5为地理编码之后的DEM。

图 5 地理编码后的StereoSAR结果DEM Fig. 5 Geocoded result DEMs of StereoSAR

图选项

为了评价图 5所示的两幅结果DEM的高程精度，本文采用Infoterra公司提供的由商业软件Pixel Factory处理生成的一幅DEM作为参考数据。该DEM通过升轨和降轨获取的各一对立体像对数据分别生成高程点阵图后经融合处理得到，且经过了异常值过滤等人工编辑处理，其采样间隔为10 m，经检验其标称平面定位精度约为1 m，高程测量精度约为5 m。将两幅DEM分别与Infoterra DEM作差值运算，得到高程差异分布统计结果如表 2所示。

统计结果表明，常规匹配方法得到的结果与Infoterra DEM之差的标准差为25.7 m；而采用外部DEM辅助得到的结果与Infoterra DEM之差，标准差减少到了17.4 m。采用外部DEM辅助之后，粗差点也明显减少，误差在25 m以内的点的比例由79%提高到86%。图 6展示了试验得到的DEM与Infoterra DEM之差在不同地形条件区域统计的直方图，可以看出，在平坦地区，两种方法得到的结果并无太大差别，与Infoterra DEM相比，标准差均能达到5 m。而在山区，采用外部DEM进行辅助的优势则较为明显，不仅大于100 m的粗差点明显减少，且高差均值也更接近于0。由此可以看出，外部DEM对雷达摄影测量提取山区地形高程的辅助作用是十分明显而有效的，其得到的结果可靠性更高。

图 6 分区域统计的高差直方图 Fig. 6 Histograms of elevation differences within two sub-areas of different terrain conditions

图选项

此外，笔者还采用了Infoterra公司提供的26个GPS地面测量点对StereoSAR生成的DEM的高程精度进行验证，这些点记录的坐标均为WGS-84参考系下的经纬度坐标。由于这些GPS点对应的影像坐标未知，只能量测结果DEM上与GPS点地理坐标相同的点位上的高程值，可以评价高程精度，但无法评价平面定位精度。由于DEM地理编码的误差，平面位置可能存在偏移，从而影响比较结果，但总体上也能初步评价本试验提取的DEM的精度。

统计结果表明，采用外部DEM辅助得到的结果相对于GPS点的高程差异均值为-1 m，差异值统计标准差为16 m。最大误差为60 m，出现在高程约为1400 m的山峰上，该点在雷达影像中位于叠掩区，其匹配结果可靠性很低，出现较大的高程误差可以理解。若将该点排除在外，DEM的精度能达到10 m。

采用雷达摄影测量方法提取DEM的精度虽然还不能与光学航空摄影测量提取DEM的精度相匹敌，但由于星载SAR系统能够全天时全天候观测获取数据，使得雷达数据资源丰富、获取便利，因此，雷达摄影测量方法也不失为提取DEM的重要途径之一，尤其是在东南亚、我国西南山区等难以获取高质量光学数据的多云多雨地区可发挥不可替代的重要作用。

在StereoSAR数据处理过程中，外部DEM的应用可为立体像对的匹配提供良好的初值，并能够抑制由于误匹配造成的粗差，从而大幅提高所生成DEM的精度。然而，由于地理编码的误差和匹配误差，目前的试验结果尚不能称作高精度。尽管如此，将外部DEM作为辅助数据在提高StereoSAR提取DEM精度上具有重要意义。进一步改善StereoSAR高程信息提取精度可以采取的措施包括：① 改进立体像对的匹配方法，提高匹配精度；② 采用控制点对定位误差进行修正，提高绝对定位精度。这将是我们今后研究中的重要内容。

致谢：本文试验采用的TerraSAR-X立体像对、GPS地面测量数据和Pixel Factory生成DEM由Infoterra公司提供，在此表示衷心感谢。