基于DEM辅助后向投影模型的InSAR高程反演方法

胡晓宁 汪丙南* 向茂生 王钟斌

①(中国科学院空天信息创新研究院 北京 100190)

②(中国科学院微波成像技术国家级重点实验室 北京 100190)

③(中国科学院大学 北京 100049)

1 引言

数字高程模型(Digital Elevation Model,DEM)在地形测绘、水文、农业等诸多领域均有着广泛的应用[1–3]。利用干涉合成孔径雷达(Interferometric Synthetic Aperture Radar,InSAR)技术反演DEM具有全天时、全天候、高效率、高精度等优势，是目前常用的DEM获取手段之一[4]。在地形起伏剧烈的区域，InSAR干涉条纹变得十分密集，特别是星载平台。在如横断山脉等地形突变区域，甚至会出现干涉相位模糊现象。密集的干涉条纹不仅增加了相位展开的难度，并且会在解缠后的干涉相位中引入更多的数据处理误差，进而影响高程反演的精度。

通过引入外源DEM可以降低干涉条纹的密度，已有一些文献对该方法进行了研究[5–8]。考虑到中国已有覆盖全国的地形数据，以及已有的各国研究团队通过各种方法获取的一些全球数字高程数据，如利用InSAR技术获取的SRTM(Shuttle Radar Topography Mission)DEM[9],利用光学成像的立体像对技术获取的ASTER GDEM(The Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model)[10]，以及通过对SRTM进行再处理，并将其与精细化后的ASTER GDEM高程相融合而生成的NASADEM[11]。上述数字高程数据可通过互联网下载获取，为在InSAR处理中引入外源DEM提供了条件。

后向投影(BackProjection,BP)算法是一种完全时域的无误差成像算法，具有聚焦效果好、相位保持精度高、运动误差补偿精度高且算法简单、适用于各种工作模式以及任意运动轨迹的雷达等优点，在合成孔径雷达(Synthetic Aperture Radar,SAR)成像中得到了广泛的应用[12]。文献[13]对基于BP算法的InSAR方法开展了相关研究，理论分析了基于BP算法和基于频域成像算法的InSAR方法的异同，但没有进行实验验证。文献[14]研究了毫米波的后向投影InSAR成像方法，并进行了实验验证，但在实测InSAR数据处理实验中仅获取了干涉相位，未进行后续的高程反演研究。后向投影成像算法也被用于获取干涉处理所需的SAR图像，以简化时序干涉处理的流程[15]。总体而言，针对基于BP算法的干涉处理目前研究较少，尤其是针对理论模型和实测InSAR数据的实验还缺乏深入的研究[16–18]。

针对上述问题，为了降低干涉相位的条纹密度，本文将外源DEM与BP算法相结合，提出了一种基于DEM辅助后向投影模型的InSAR高程反演方法。该方法利用BP算法获得聚焦的SAR图像，并在成像过程中引入辅助DEM，然后建立基于BP成像算法的InSAR高程反演模型，实现高精度的DEM反演。外源DEM辅助以及BP算法的引入，不仅可以降低干涉条纹的密度，减少干涉相位的缠绕，多数情况下甚至可以避免相位缠绕，同时图像对无需干涉配准，简化了干涉处理的流程。仿真实验与机载InSAR系统获取的实测数据处理结果证明了算法的有效性。

2 后向投影InSAR高程反演模型

传统InSAR处理中多采用频域算法成像，获取的干涉相位中包含平地相位和高程相位。平地效应会大大增加干涉条纹的密集程度，去平地操作虽然可以降低干涉条纹的密度，但在地形起伏较大的区域干涉条纹仍然十分密集，甚至出现相位模糊的问题，进而影响高程反演。在本文方法中，外源DEM的辅助和BP算法的引入大大降低了干涉条纹的密度，并且多数情况下避免了相位缠绕。本节将从成像处理和干涉处理两个环节对该方法进行详细介绍，并给出具体的流程。

2.1 后向投影成像算法

BP成像算法的基本原理是对距离压缩后的数据精确补偿投影面上的采样点到天线相位中心的回波延时相位，并逐点逐脉冲地进行相干累积，从而得到聚焦图像[19]。

主要包括以下步骤：

(1)根据成像区域确定成像空间，对成像空间进行网格划分，并将辅助DEM插值到划分的网格中作为BP成像过程中的投影面。

(2)对原始回波数据进行距离压缩。

(3)计算投影面上的采样点到天线相位中心的距离延时，进行多普勒相位补偿，对合成孔径内每个天线相位中心处相位补偿后的数据进行相干累加，以实现该采样点的方位向聚焦。

(4)遍历成像空间中的每一个采样点，得到聚焦成像的结果。

BP算法不存在雷达和目标的斜距近似误差，成像信号模型十分精确，相位保持精度也很高，并且BP算法根据雷达轨迹信息进行运动误差补偿，具有很高的补偿精度，为SAR图像在后续的干涉处理中获取高精度的DEM提供了基础；另外BP算法的成像过程处于空间直角坐标系下，易于引入外源辅助DEM，并且可以避免复杂的图像拼接和地理编码过程。这也是选择BP算法进行成像处理的原因。

2.2 基于DEM辅助后向投影模型的InSAR高程反演方法

由于BP算法的方位向聚焦步骤中补偿的多普勒相位项包含了网格点到雷达平台的距离历史信息，而这正是传统InSAR中干涉相位的来源，因此基于BP成像算法的干涉几何模型有别于传统InSAR几何模型。考虑到实际情况下辅助DEM存在地形起伏，本文建立了基于BP算法曲面投影下的InSAR高程反演模型，下面进行具体的分析。

成像过程中引入的辅助DEM为粗精度DEM，与目标的真实高程之间存在误差，因此目标在投影曲面上的投影位置将偏离目标的真实位置，并且由于BP成像过程中网格的离散化采样，目标的投影位置也会偏离成像网格的采样点位置，因此目标真实的散射中心与成像过程中投影的网格点分别到雷达天线的斜距存在偏差，进而导致BP成像方位向聚焦过程中多普勒相位补偿存在误差。该误差正是基于时域BP算法获取的SAR图像进行干涉处理时干涉相位的来源，几何模型如图1所示。

图1 基于BP成像算法的干涉模型Fig.1 The interferometric model based on BP imaging algorithm

图1中z轴为高度向，y轴指向雷达平台的飞行方向，x轴由右手坐标系确定。A1和A2分别代表主、副天线相位中心的位置，下视角为φ，基线长度B，基线与水平方向之间的夹角为α。由于目标P的真实高度与辅助DEM之间存在高度偏差，主副天线将P分别投影到Q1,Q2，而成像网格中对应的离散采样点为C。天线相位中心A1的高程为H。点C,Q1和Q2的高程值分别为投影曲面上的采样点C与目标点P的真实位置之间的高程偏差为Δh,θ为与垂直方向的夹角。

根据余弦定理可得

在传统InSAR的干涉相位模型中，P和Q1两点间的干涉相位差与其高程差有关，根据高程相位的表达式[20]，该干涉相位差如式(6)所示。

根据图1中的几何关系，Q1和C点之间的干涉相位差可以表示为引起的平地相位和引起的高程相位之和的形式，如式(7)所示。

分析式(8)，后向投影InSAR高程反演模型得到的干涉相位由两部分组成：Δh的高度差引起的干涉相位差对应的斜距变化引起的干涉相位差Δφr。在实际情况中，Δφr远小于Δφh，可以忽略不计(以3.2节的实验参数为例，根据式(8)计算高程误差为10m时对应的Δφr与Δφh分别为–0.0166rad和–1.7rad)。因此最终得到后向投影干涉模型下高程反演公式为

由式(9)可以看出，本文基于后向投影算法的InSAR高程反演模型获得的干涉相位去除了成像参考面对应的地形相位，仅包含辅助DEM的高程误差对应的高程相位，因此无需去平地即可获得条纹密度大大降低的干涉相位图，并且与基于频域成像算法的传统InSAR高程反演方法中得到的去平地后的干涉相位相比，其条纹更加稀疏。

下面分析本文方法中相位缠绕的发生条件。图2中P1,P2为干涉相位发生2π缠绕的两个点目标，h1,he1分别为P1的真实高程值和辅助DEM对应点处的高程值，h2,he2分别为P2的真实高程值和辅助DEM对应点处的高程值，Δh1和Δh2则分别为P1和P2辅助DEM对应点处的高程误差。根据式(9)计算基于后向投影模型得到的干涉相位变化2π，即出现相位缠绕现象时对应的辅助DEM高程误差变化值，用Δh2π表示，如式(10)所示。Δh2π对应图2中的Δh2−Δh1。根据高程模糊度的定义以及图2中的几何关系，推导本文基于辅助DEM的后向投影InSAR高程反演模型的高度模糊数计算公式，如式(11)所示。与传统InSAR方法中的高度模糊数相比，本文方法的高度模糊数更大，即干涉条纹更加稀疏，且更不易发生干涉相位缠绕。

图2 高程模糊度几何示意图Fig.2 Geometric description of the height of ambiguity

以3.2节的X波段机载SAR系统参数为例，根据式(10)计算得到Δh2π大于100m，即辅助DEM的高程误差变化大于100m时，干涉相位才会缠绕，而对于目前可以获取的外源DEM而言，其精度远高于该值，因此即使在地形起伏较大的区域也不会出现相位缠绕的现象，从而避免了相位解缠处理。在实际情况中，考虑到基线误差会引入随距离向变化的相位误差，当辅助DEM高程误差与基线误差导致的干涉相位小于2π时，基于后向投影算法的In-SAR模型得到的干涉相位不会出现相位缠绕，避免相位解缠过程。

基于DEM辅助后向投影模型的InSAR高程反演方法在简化干涉处理流程的同时，减少了解缠过程中的相位损失，提高了DEM反演的精度。综上所述，本文所提的高程反演方法具有精确的成像信号模型、运动误差补偿模型以及很高的相位保持能力。

2.3 基于DEM辅助后向投影模型的InSAR高程反演流程

本文给出了一种基于DEM辅助后向投影模型的InSAR高程反演方法，处理流程如图3所示，并与传统基于频域算法的InSAR处理流程进行对比。传统InSAR采用频域成像算法得到聚焦的SAR图像后，需要进行图像配准、去平地、干涉相位滤波、相位解缠等处理，然后根据高程反演模型得到DEM重建结果。而基于DEM辅助后向投影模型的InSAR高程反演方法简化了数据处理流程，无需图像配准以及相位解缠操作，减少了数据处理过程中的相位损失，为获取高精度DEM提供了保障。并且基于后向投影的InSAR模型位于地理空间坐标系下，可以避免后续复杂的图像拼接和三维定位等处理。具体步骤如下：

图3 InSAR处理流程Fig.3 InSAR processing chain

(1)BP成像。获取实验场景的外源DEM数据，并插值作为成像的参考面，利用时域的后向投影算法获取聚焦的SAR图像。

(2)图像干涉。主副天线获取的复图像直接共轭相乘得到后向投影InSAR模型下的干涉相位。

(3)相位滤波。为了降低相位噪声对高程反演精度的影响，对干涉相位进行滤波处理。

(4)高程反演。根据式(9)，实现干涉相位到高程的转换。由于式(9)计算的是BP算法成像过程中参考面相对真实地形的高程偏差，需要将其加回到成像过程中采用的辅助DEM，即可得到最终反演的DEM结果。

3 实验及分析

3.1 仿真实验

利用仿真实验验证基于DEM辅助后向投影模型的InSAR高程反演方法的有效性及优越性。仿真实验的参数如表1所示。实验模拟地形起伏较大的区域，其高程设置如图4(a)所示。根据InSAR回波信号模型，仿真生成主副天线的回波信号，并加入–30dB的高斯白噪声。实验过程中分别采用了本文所提的基于后向投影的InSAR方法与传统InSAR方法对回波信号进行成像和干涉处理。

表1 仿真参数Tab.1 Simulation parameters

首先利用后向投影InSAR高程反演方法，引入存在高程误差的辅助DEM作为BP算法成像的参考面并生成聚焦的SAR图像。干涉相位由主辅图像对直接共轭相乘即可得到，无需图像配准，如图4(b)所示。根据2.2节的分析，基于后向投影算法的In-SAR高程反演模型获得的干涉相位已去除了成像过程中大部分的地形相位，因此干涉条纹密度大大降低。同时由于辅助DEM的高程误差较小，残余的地形相位数值较小不足以引起干涉相位的缠绕，如图4(b)所示，避免了相位解缠操作。

图4 仿真实验结果图Fig.4 Simulation results

然后根据传统InSAR高程反演方法，采用距离多普勒(RD)算法成像，对得到的复图像对进行配准、去平地以及相位滤波处理，获取的干涉相位如图4(c)所示。由于仿真区域地形起伏较大，因此干涉条纹十分密集，并且在高程剧烈变化的区域出现了条纹模糊的问题，如图4(c)中红色矩形框区域所示，采用常规的干涉处理将无法得到该区域正确的相位解缠结果。

针对图4(c)中的相位模糊问题，本文在传统InSAR方法中进一步引入DEM辅助相位解缠处理。将根据SAR图像生成的原始干涉相位与基于辅助DEM模拟的干涉相位进行差分处理，得到残余干涉相位，如图4(d)所示。该干涉相位去除了辅助DEM对应的地形相位，大大降低了干涉条纹的密度。将残余干涉相位与模拟的干涉相位求和即可得到原始干涉相位图的解缠结果。

由图4(b)和图4(d)可以看出，辅助DEM的引入避免了条纹密集区域的相位模糊问题，并且可以得到正确的相位解缠结果。但是在传统InSAR方法中引入辅助DEM需要模拟DEM对应的干涉相位，并需要对模拟相位与原始干涉相位进行配准，增加了干涉处理流程的复杂度。而本文所提后向投影In-SAR高程反演方法在解决相位模糊问题的同时简化了干涉处理流程，因此具有更大的优势。

在生成InSAR回波信号时，仿真的实验场景中设置了25个强散射点作为地面控制点，用于统计高程反演的误差。本文所提后向投影InSAR方法的实验结果如图5(a)所示，25个控制点处的高程反演误差的均值为–0.0326m，标准差为0.2510m。

图5 两种方法的高程反演误差Fig.5 The elevation inversion errors of two methods

进一步对比传统InSAR高程反演方法，采用RD算法成像并在相位解缠过程中引入DEM辅助。统计控制点处的高程反演误差，如图5(b)所示，均值为0.0349m，标准差为0.3484m。由于BP算法的成像模型更加精确，相位保持精度更高，因此后向投影InSAR高程反演精度略优于基于RD成像的传统InSAR方法，即本文所提简化的基于辅助DEM后向投影InSAR高程反演方法可以实现高精度DEM的反演。

仿真实验验证了本文方法在反演高精度DEM的同时，无需图像配准操作，大大降低干涉条纹的密度，甚至避免了相位解缠步骤，简化了干涉处理的流程。在地形起伏剧烈的区域，基于后向投影成像模型的InSAR高程反演方法通过辅助DEM的引入，可以有效地避免相位模糊问题。

3.2 机载InSAR实测数据处理

(1)平地区域实验结果

在本节中，采用第2节提出的基于后向投影的InSAR高程反演方法对实测机载InSAR数据进行处理。实验所用的机载InSAR数据是在中国内蒙古获取的X波段双天线数据，高程测量标称精度1m，具体参数如表2所示。

表2 X波段机载InSAR系统参数Tab.2 X-band airborne InSAR system parameters

实验中选取的成像区域大小为5.0km×1.6km，引入30m分辨率的SRTM DEM作为辅助DEM，利用BP算法获取聚焦的SAR图像，并与RD算法得到的SAR图像对比。选取成像场景中的强散射点目标，对该目标邻近区域进行16倍FFT插值，图6(a)和图6(b)分别是插值后地面强散射点处BP算法和RD算法的主辅图像对成像结果，图6(c)和图6(d)分别为目标沿距离向的能量分布。可以看出，BP算法得到的主辅图像已配准，而RD算法得到的主辅图像对中目标位置在距离向偏移了59个单元。插值结果证明了本文所提后向投影InSAR高程反演模型可以实现自配准，避免图像配准过程。

图6 强散射目标点成像结果Fig.6 The focused images of intense scatterer

BP算法得到的主图像幅度图如图7(a)所示。实验过程中在该成像区域布放的地面控制点位置信息如图7(a)中红色三角形标注。由于主辅图像已实现配准，因此直接共轭相乘即可得到干涉相位，如图7(b)所示。根据第2节的理论分析，基于DEM辅助后向投影模型的InSAR高程反演方法由于在成像过程中引入了外源DEM，可直接得到去地形相位后的干涉相位，该干涉相位由辅助DEM成像参考面与真实地形间高程偏差引起，因此残余干涉相位数值较小，如图7(b)所示，变化范围小于2π/5，避免了相位解缠处理步骤。相位滤波后直接进行高程反演，辅助DEM的高程偏差如图7(c)所示。将计算的偏差与辅助DEM的高程值相加，最终获得成像区域的高程反演结果，如图7(d)所示。

图7 机载InSAR数据处理结果图Fig.7 Airborne InSAR data processing results

利用实验中布放的部分地面控制点(未参与干涉定标)作为检查点，对传统方法与本文方法高程反演误差进行对比验证，如表3所示。对比结果表明，传统InSAR成像和干涉处理结果能够达到1m标称高程精度要求(高程精度0.85m)，然而利用本文提出的DEM辅助干涉处理算法可以达到更高的处理精度(高程精度0.65m)。由于后向投影成像和干涉处理几何，是一种完全时域的无误差成像算法，在成像信号模型、运动误差补偿模型上理论最为精确，相位保持精度更高；同时在外源DEM辅助下，避免配准和相位解缠过程，减少处理过程中的相位损失。因此，本文方法处理精度优于传统InSAR成像和干涉处理算法，可以获得高精度的DEM反演结果，且具有简化的干涉处理流程。

表3 地面检查点处高程反演误差Tab.3 The elevation inversion errors of ground detection points

(2)丘陵区域实验结果

本文另外选取了丘陵区域的实验数据，进一步验证所提的基于DEM辅助后向投影模型的In-SAR高程反演方法。BP算法得到的主辅图像对直接共轭相乘得到干涉相位如图8(a)所示。虽然实验区域地形起伏较大，但由于辅助DEM的引入，基于后向投影的InSAR模型在成像过程去除了绝大多数的地形相位，图8(a)中的干涉相位范围小于π，无需相位解缠，简化了后续的干涉处理。而传统InSAR方法得到的干涉相位如图8(b)所示，由于平地效应与地形起伏，图8(b)中存在3个缠绕的干涉条纹，与图8(a)相比后续干涉处理步骤更加复杂。图8(c)是利用BP算法得到的幅度图。根据DEM辅助下基于后向投影的InSAR模型得到的干涉相位及式(9)反演辅助DEM的偏差，并加回到辅助DEM，最终得到该成像区域的DEM如图8(d)所示。由丘陵区域实验数据的处理结果可以看出，本文基于辅助DEM与BP成像算法的InSAR方法可以有效地减少干涉相位的缠绕，避免相位解缠处理，验证了所提基于后向投影的InSAR高程反演模型的有效性。

图8 机载InSAR山区数据处理结果图Fig.8 Airborne InSAR data processing results of mountainous area

4 结论

本文提出了一种基于DEM辅助后向投影模型的InSAR高程反演方法。该方法利用BP算法运动补偿精度高和保相性能好的优势，同时在BP成像过程中引入外源DEM辅助，降低干涉条纹的密度，在多数情况下还实现了干涉处理流程的简化—无需图像配准与相位解缠操作。本文最后设计了仿真实验，并对X波段双天线机载In-SAR数据进行了处理，通过所提方法中简化的干涉处理流程实现了高精度DEM的反演，验证了该方法的有效性。然而，在基于BP算法的干涉模型中还存在一些问题，例如实际的InSAR系统中天线相位中心存在测量误差，在高分辨成像时需要考虑对方位向配准产生的影响，以及本文模型的高程反演误差理论分析等，这将是今后的研究方向。