雷达学报 Journal of Radars ISSN 2095 -283X,CN 10-1030/TN      《雷达学报》网络首发论文  题目：  基于DEM辅助后向投影模型的 InSAR 高程反演方法 作者：  胡晓宁，汪丙南，向茂生，王钟斌 收稿日期：  2020 -11 -28  网络首发日期：  2021 -03 -26 引用格式：  胡晓宁，汪丙南，向茂生，王钟斌．基于 DEM辅助后向投影模型的 InSAR高程反演方法[J/OL] ．雷达学报.  https://kns.cnki.net/kcms/detail/10.1030.TN.20210324.1406.002.html        网络首发：在编辑部工作流程中，稿件从录用到出版要经历录用定稿、排版定稿、整期汇编定稿等阶段。录用定稿指内容已经确定，且通过同行评议、主编终审同意刊用的稿件。排版定稿指录用定稿按照期刊特定版式（包括网络呈现版式）排版后的稿件，可暂不确定出版年、卷、期和页码。整期汇编定稿指出版年、卷、期、页码均已确定的印刷或数字出版的整期汇编稿件。录用定稿网络首发稿件内容必须符合《出版管理条例》和《期刊出版管理规定》的有关规定；学术研究成果具有创新性、科学性和先进性，符合编辑部对刊文的录用要求，不存在学术不端行为及其他侵权行为；稿件内容应基本符合国家有关书刊编辑、出版的技术标准，正确使用和统一规范语言文字、符号、数字、外文字母、法定计量单位及地图标注等。为确保录用定稿网络首发的严肃性，录用定稿一经发布，不得修改论文题目、作者、机构名称和学术内容，只可基于编辑规范进行少量文字的修改。 出版确认：纸质期刊编辑部通过与《中国学术期刊（光盘版）》电子杂志社有限公司签约，在《中国学术期刊（网络版）》出版传播平台上创办与纸质期刊内容一致的网络版，以单篇或整期出版形式，在印刷出版之前刊发论文的录用定稿、排版定稿、整期汇编定稿。因为《中国学术期刊（网络版）》是国家新闻出版广电总局批准的网络连续型出版物（ISSN  2096-4188，CN 11-6037/Z ），所以签约期刊的网络版上网络首发论文视为正式出版。  基于DEM辅助后向投影模型的InSAR高程反演方法 胡晓宁①②③    汪丙南*①②③    向茂生①②③    王钟斌①②③ ①（中国科学院空天信息创新研究院 北京 100190） ②（中国科学院微波成像技术国家级重点实验室 北京 100190） ③（中国科学院大学 北京 100049） 摘 要：利用干涉合成孔径雷达(InSAR)技术获取数字高程模型(DEM)时，在地形起伏剧烈区域，干涉条纹十分密集，增加了相位解缠的难度，影响相位展开和高程反演的精度。为了解决该问题，本文提出了一种基于DEM辅助后向投影模型的InSAR高程反演方法。该方法可以在统一的后向投影成像空间中实现成像和InSAR高程反演，并且引入外源DEM作为辅助信息，去除大部分地形相位，有效的降低了干涉条纹的密度，减少了干涉相位的缠绕。此外，该方法在多数情况下可以避免图像配准和相位解缠过程，简化了传统InSAR的处理流程，并且可以实现高精度的高程反演。通过仿真实验和X波段机载双天线InSAR数据的处理验证了该方法的有效性。 关键词：高程反演； InSAR；辅助DEM；BP算法；SAR 中图分类号：TN958 DOI: 10.12000/JR20144 引用格式：胡晓宁，汪丙南，向茂生, 等. 基于DEM辅助后向投影模型的InSAR高程反演方法[J]. 雷达学报, 待出版. DOI: 10.12000/JR20144.   Reference format: HU Xiaoning，WANG Bingnan，XIANG Maosheng, et al. InSAR Elevation Inversion Method Based on Backprojection Model with External DEM [J]. Journal of Radars, in press. DOI: 10.12000/JR20144. InSAR Elevation Inversion Method Based on Backprojection Model with External DEM HU Xiaoning①②③  WANG Bingnan*①②③  XIANG Maosheng①②③  WANG Zhongbin①②③ ①(Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100190, China) ②(National Key Laboratory of Microwave Imaging Technology, Beijing 100190, China) ③(University of Chinese Academy of Sciences, Beijing 100049, China) Abstract: When interferometric synthetic aperture radar (InSAR) is used to obtain the digital elevation model (DEM), highly sloped terrain will make interferometric fringes dense and increase the  difficulty  of  phase  unwrapping,  which  will  affect  the  accuracy  of  phase  unwrapping  and elevation inversion. In order to solve this problem, an InSAR elevation inversion method based on backprojection (BP) model with external DEM is proposed. This model achieves  imaging and InSAR DEM inversion in a uniform backprojection geographic space and introduces external DEM as auxiliary information, which can remove most of the phase of terrain and reduce the density of interferometric  fringes and phase  wrapping.  Additionally,  the  proposed  method  can  avoid  the procedure  of  image  registration  and  phase  unwrapping  in  most  case,  which  simplifies  the traditional InSAR processing and achieves higher processing accuracy. Simulation experiment and X-band InSAR data processing have been carried out to verify the effectiveness of this method.                                                              收稿日期：2020-11-28*通信作者：汪丙南 wbn@mail.ie.ac.cn  *Corresponding author: WANG Bingnan, wbn@mail.ie.ac.cn 基金项目：国家自然科学基金资助项目(62073306)，中国科学院青年创新促进会项目资助，国家重点研发计划资助( 2017YFC0822400) Foundation Item: The National Natural Science Foundation of China (62073306), Youth Innovation Promotion Association CAS, National Key R&D Program of China (2017YFC0822400) 责任主编：靳国旺    Corresponding Editor: JIN Guowang 网络首发时间：2021-03-26 11:06:14网络首发地址：https://kns.cnki.net/kcms/detail/10.1030.TN.20210324.1406.002.htmlKey word: Elevation inversion; InSAR; External DEM; Backprojection (BP) algorithm; SAR 1．引言 数字高程模型(Digital Elevation Model, DEM)在地形测绘、水文、农业等诸多领域均有着广泛的应用[1]-[3]。利用干涉合成孔径雷达(Interferometric Synthetic Aperture Radar, InSAR)技术反演DEM具有全天时、全天候、高效率、高精度等优势，是目前常用的DEM获取手段之一[4]。在地形起伏剧烈的区域，InSAR干涉条纹变得十分密集，特别是星载平台。在如横断山脉等地形突变区域，甚至会出现干涉相位模糊现象。密集的干涉条纹不仅增加了相位展开的难度，并且会在解缠后的干涉相位中引入更多的数据处理误差，进而影响高程反演的精度。 通过引入外源DEM可以降低干涉条纹的密度，已有一些文献对该方法进行了研究[5]- [8]。考虑到在中国已有覆盖全国的地形数据，以及已有的各国研究团队通过各种方法获取的一些全球数字高程数据，如利用InSAR 技术获取的SRTM  DEM[9],  利用光学成像的立体像对技术获取的ASTER  GDEM[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分别代表主、副天线相位中心的位置，下视角为�，基线长度�，基线与水平方向之间的夹角为�。由于目标�的真实高度与辅助DEM之间存在高度偏差，主副天线将�分别投影到Q1，Q2，而成像网格中对应的离散采样点为C。天线相位中心A1的高程为�。  点C, Q1  和Q2  的高程值分别为  �c,  �Q1  和�Q2。投影曲面上的采样点C与目标点P的真实位置之间的高程偏差为∆ℎ， �为A1Q1̅̅̅̅̅̅̅与垂直方向的夹角。 首先分析图像对偏移量Q1Q2̅̅̅̅̅̅̅在XY平面上的投影Q1Q2̅̅̅̅̅̅̅Proj。根据图1中的几何关系，A1P̅̅̅̅̅=A1Q1̅̅̅̅̅̅̅,A2P̅̅̅̅̅=A2Q2̅̅̅̅̅̅̅，可得到下式：  ( ) ( ) aa æö = - - - + - + - ç÷ èø 1222 22 1 2 1 2 Proj QQ QQ AP cos AP sin H z B H B z (1)  j - D = -c1APsin h H z   (2) 根据余弦定理可得：  ( ) aj = + - × × +222 1 1AP AP 2 AP cos BB   (3) 由式(1)-(3)，计算Q1Q2̅̅̅̅̅̅̅Proj对∆ℎ的泰勒展开形式，得到图像对偏移量的表达式如式(4)，其中OQ2̅̅̅̅̅̅Proj是OQ2̅̅̅̅̅̅在XY平面上的投影。  ajj+» - DProjPro122 jsin( )QQQ cos OBh   (4) 由式(4)可以看出，图像对偏移量Q1Q2̅̅̅̅̅̅̅Proj与辅助DEM的高程误差有关，针对目前可以获取的外源DEM精度，其高程误差引起的图像对偏移量少于一个分辨单元的大小。以第3.2 节中的机载InSAR的参数为例，根据式(4)可以计算当投影曲面高程误差为10m时，图像对偏移量大概为1e-3m，即只有当∆ℎ达到百米量级时，图像对偏移量才可能超过一个分辨单元。因此可以认为经离散采样后点目标P位于BP成像网格的同一采样点处。也就是说，基于BP算法的InSAR模型在成像过程中可以实现自配准，无需后续图像配准操作，主辅图像直接共轭相乘即可得到干涉相位，简化了干涉处理的流程。 基于上述结论，下面推导基于DEM辅助后向投影模型的InSAR高程反演公式。在BP成像过程中，主图像和辅图像理论上需要补偿的多普勒相位项分别为−4�A1P̅̅̅̅̅/�和−4�A2P̅̅̅̅̅/�。但在实际的成像过程中，补偿的多普勒相位项分别为4�A1C̅̅̅̅̅/� 和4�A2C̅̅̅̅̅/�，因此由于BP成像网格化导致主辅图像对之间存在下式所示的相位差。  ( ) ( )( ) ( ) ( ) ( )pp jll pp ll j j j jé ù é ùD = - - -ê ú ê úë û ë ûé ù é ù= - -ê ú ê úë û ë û+ = D - D D - D111 1 2 21 2 1 1 2 1 1 1 2 1 1 2P Q Q C44AP-AC AP-AC44AP-AP - AQ -AQ AQ -AQ - AC) ( (-AC)  (5) 在传统InSAR的干涉相位模型中，P和Q1两点间的干涉相位差与其高程差�ℎPQ1有关，根据高程相位的表达式[20]，该干涉相位差如式(6)所示。  ( ) qa pjjl q-D - D = - D11 P Q PQ1cos4APsinBh   (6) 根据图1中的几何关系，Q1和C点之间的干涉相位差可以表示为∆�CQ1引起的平地相位和�ℎCQ1引起的高程相位之和的形式，如式(7)所示。其中，�ℎCQ1和∆�CQ1分别是Q1和C两点之间的高程差与到雷达平台的斜距差。  ( ) qa p p q ajjl q l q- -D - D = - D - D1 1 1 Q C CQ CQ1 1 1cos4 4 cos( )APsin AQ tanB Bhr   (7) 在图1中，�ℎ=�ℎPQ1+�ℎCQ1，所以干涉相位可以写为式(8)的形式。  ( ) qa p p q ajl q l qjj--D = - D - D= D + D1CQ1 1 1cos4 4 cos( )APsin AQ tanhrB Bhr   (8) 当成像投影面为平面时，即为�ℎCQ1=0的特例，式(8)仍然成立。 分析式(8)，后向投影InSAR高程反演模型得到的干涉相位由两部分组成——∆ℎ的高度差引起的干涉相位差∆�ℎ和CQ1̅̅̅̅̅对应的斜距变化引起的干涉相位差∆��。在实际情况中，∆��远小于∆�ℎ，可以忽略不计（以第3.2 节中的实验参数为例，根据式(8)计算高程误差为10m时对应的∆��与∆�ℎ分别为-0.0166rad和-1.7rad）。因此最终得到后向投影干涉模型下高程反演公式如式(9)。  ( )1cos4AP sinBhqa pjl q-D » - × × D×   (9) 由式(9)可以看出，本文基于后向投影算法的InSAR高程反演模型获得的干涉相位去除了成像参考面对应的地形相位，仅包含辅助DEM的高程误差对应的高程相位，因此无需去平地即可获得条纹密度大大降低的干涉相位图，并且与基于频域成像算法的传统InSAR高程反演方法中得到的去平地后的干涉相位相比，其条纹更加稀疏。  图  2  高程模糊度几何示意图 Fig. 2 Geometric description of the height of ambiguity 下面分析本文方法中相位缠绕的发生条件。图  2为高程模糊度几何示意图，P1、P2为干涉相位发生2π缠绕的两个点目标，h1、he1分别为P1的真实高程值和辅助DEM对应点处的高程值，h2、he2分别为P2的真实高程值和辅助DEM对应点处的高程值，Δh1和Δh2则分别为P1和P2辅助DEM对应点处的高程误差。根据式(9)计算基于后向投影模型得到的干涉相位变化2π，即出现相位缠绕现象时对应的辅助DEM高程误差变化值，用Δh2π表示，如式(10)所示。Δh2π对应图2中的∆ℎ2−∆ℎ1。根据高程模糊度的定义以及图2中的几何关系，推导本文基于辅助DEM的后向投影InSAR高程反演模型的高度模糊数计算公式，如式(11)所示。与传统InSAR方法中的高度模糊数相比，本文方法的高度模糊数更大，即干涉条纹更加稀疏，且更不易发生干涉相位缠绕。  ( )21AP sin2 coshBplqqaD×=×-   (10)  ( )2 2 12 1 e2 e1e2 e11AP sin2 cos( ) ( )()h h hh h h hhhBplqqa=-= D - D + -= + -××-   (11) 以第3.2 节中的X波段机载SAR系统参数为例，根据式(10)计算得到Δh2π大于100m，即辅助DEM的高程误差变化大于100m时，干涉相位才会缠绕，而对于目前可以获取的外源DEM而言，其精度远高于该值，因此即使在地形起伏较大的区域也不会出现相位缠绕的现象，从而避免了相位解缠处理。在实际情况中，考虑到基线误差会引入随距离向变化的相位误差，当辅助DEM高程误差与基线误差导致的干涉相位小于2π时，基于后向投影算法的InSAR模型得到的干涉相位不会出现相位缠绕，避免相位解缠过程。 基于DEM辅助后向投影模型的InSAR 高程反演方法在简化干涉处理流程的同时，减少了解缠过程中的相位损失，提高了DEM反演的精度。综上所述，本文所提的高程反演方法具有精确的成像信号模型、运动误差补偿模型以及很高的相位保持能力。 2.3 基于DEM辅助后向投影模型的InSAR高程反演流程 本文给出了一种基于DEM辅助后向投影模型的InSAR 高程反演方法，处理流程如图  3所示，并与传统基于频域算法的InSAR 处理流程进行对比。传统InSAR 采用频域成像算法得到聚焦的SAR图像后，需要进行图像配准、去平地、干涉相位滤波、相位解缠等处理，然后根据高程反演模型得到DEM重建结果。而基于DEM辅助后向投影模型的InSAR 高程反演方法简化了数据处理流程，无需图像配准以及相位解缠操作，减少了数据处理过程中的相位损失，为获取高精度DEM提供了保障。并且基于后向投影的InSAR模型位于地理空间坐标系下，可以避免后续复杂的图像拼接和三维定位等处理。具体步骤如下：   (1)  BP成像。获取实验场景的外源DEM数据，并插值作为成像的参考面，利用时域的后向投影算法获取聚焦的SAR图像。 (2)  图像干涉。主副天线获取的复图像直接共轭相乘得到后向投影InSAR模型下的干涉相位。 (3)  相位滤波。为了降低相位噪声对高程反演精度的影响，对干涉相位进行滤波处理。 (4)  高程反演。根据式(9)，实现干涉相位到高程的转换。由于式(9)计算的是BP算法成像过程中参考面相对真实地形的高程偏差，需要将其加回到成像过程中采用的辅助DEM，即可得到最终反演的DEM结果。 回波数据去平地频域成像图像配准相位滤波相位解缠DEM重建DEM回波数据后向投影成像相位滤波DEM重建DEM辅助DEM (a)  本文方法处理流程     (b) 传统InSAR处理流程 (a) The chain of proposed method   (b) The chain of traditional InSAR 图  3  干涉处理流程 Fig. 3 InSAR processing chain  3．实验及分析 3.1 仿真实验 利用仿真实验验证基于DEM辅助后向投影模型的InSAR 高程反演方法的有效性及优越性。仿真实验的参数如表  1所示。实验模拟地形起伏较大的区域，其高程设置如图  4(a)所示。根据InSAR回波信号模型，仿真生成主副天线的回波信号，并加入-30dB的高斯白噪声。实验过程中分别采用了本文所提的基于后向投影的InSAR方法与传统InSAR方法对回波信号进行成像和干涉处理。 表  1  仿真参数 Tab. 1 Simulation parameters 成像仿真参数  参数值  载频  9.6GHz 距离向带宽 100MHz 距离向采样率 120MHz 脉冲宽度  3.7μs 脉冲重复频率 300Hz 平台平均速度 113.5m/s 平台平均高度 3286.5m 中心下视角 0.8727rad 基线长度  2.189m 基线倾角  0  首先利用后向投影InSAR高程反演方法，引入存在高程误差的辅助DEM作为BP算法成像的参考面并生成聚焦的SAR图像。干涉相位由主辅图像对直接共轭相乘即可得到，无需图像配准，如图  4 (b)所示。根据2.2 节的分析，基于后向投影算法的InSAR高程反演模型获得的干涉相位已去除了成像过程中大部分的地形相位，因此干涉条纹密度大大降低。同时由于辅助DEM的高程误差较小，残余的地形相位数值较小不足以引起干涉相位的缠绕，如图  4 (b)所示，避免了相位解缠操作。 然后根据传统InSAR高程反演方法，采用距离多普勒（RD）算法成像，对得到的复图像对进行配准、去平地以及相位滤波处理，获取的干涉相位如图  4(c)所示。由于仿真区域地形起伏较大，因此干涉条纹十分密集，并且在高程剧烈变化的区域出现了条纹模糊的问题，如图  4 (c)中红色矩形框区域所示，采用常规的干涉处理将无法得到该区域正确的相位解缠结果。 针对图  4 (c)中的相位模糊问题，本文进一步在传统InSAR方法中引入DEM辅助相位解缠处理。将根据SAR图像生成的原始干涉相位与基于辅助DEM模拟的干涉相位进行差分处理，得到残余干涉相位，如图  4  (d)所示。该干涉相位去除了辅助DEM对应的地形相位，大大降低了干涉条纹的密度。将残余干涉相位与模拟的干涉相位求和即可得到原始干涉相位图的解缠结果。 由图  4 (b)与图  4 (d)可以看出，辅助DEM的引入避免了条纹密集区域的相位模糊问题，并且可以得到正确的相位解缠结果。但是在传统InSAR方法中引入辅助DEM需要模拟DEM对应的干涉相位，并需要对模拟相位与原始干涉相位进行配准，增加了干涉处理流程的复杂度。而本文所提后向投影InSAR 高程反演方法在解决相位模糊问题的同时简化了干涉处理流程，因此具有更大的优势。   (a)  仿真场景的高程图  (a) The height of simulation scene  (b) 本文所提基于DEM辅助后向投影模型的InSAR方法获取的干涉相位图 (b) Interferogram obtained by InSAR method based on backprojection model with external DEM  (c) 传统InSAR方法去平地后的干涉相位    (c) Interferogram removing flat terrain obtained by traditional InSAR  (d) 传统InSAR方法引入辅助DEM后的残余干涉相位 (d) Residual interferogram obtained by traditional InSAR based on external DEM 图  4  仿真实验结果图 Fig. 4 Simulation results  在生成InSAR回波信号时，仿真的实验场景中设置了25个强散射点作为地面控制点，用于统计高程反演的误差。本文所提后向投影InSAR方法的实验结果如图  5(a)所示，25个控制点处的高程反演误差的均值为-0.0326m，标准差为0.2510m。 进一步对比传统InSAR高程反演方法，采用RD算法成像并在相位解缠过程中引入DEM辅助。统计控制点处的高程反演误差，如图  5(b)所示，均值为0.0349m，标准差为0.3484m。由于BP算法的成像模型更加精确，相位保持精度更高，因此后向投影InSAR高程反演精度略优于基于RD成像的传统InSAR方法，即本文所提的简化的基于辅助DEM的后向投影InSAR高程反演方法可以实现高精度DEM的反演。   (a) 本文所提基于DEM辅助后向投影模型的InSAR方法的高程反演误差 (a) The elevation inversion errors of InSAR method based on backprojection model with external DEM  (b) DEM辅助下传统InSAR方法的高程反演误差 (b) The elevation inversion errors of traditional InSAR method with external DEM 图  5  两种方法的高程反演误差 Fig. 5 The elevation inversion errors of two methods  仿真实验验证了本文方法在反演高精度DEM的同时，无需图像配准操作，大大降低干涉条纹的密度，甚至避免了相位解缠步骤，简化了干涉处理的流程。在地形起伏剧烈的区域，基于后向投影成像模型的InSAR高程反演方法通过辅助DEM的引入，可以有效地避免相位模糊问题。  3.2  机载InSAR实测数据处理 一、平地区域实验结果 在本节中，采用第2节中提出的基于后向投影的InSAR高程反演方法对实测机载InSAR数据进行处理。实验所用的机载InSAR 数据是在中国内蒙古获取的X波段双天线数据，高程测量标称精度1m，具体参数如表  2所示。 表  2 X波段机载InSAR系统参数 Tab. 2 X-band airborne InSAR system parameters 机载InSAR参数  参数值 载频  9.6GHz 距离向带宽  300MHz 距离向采样率  500MHz 脉冲宽度  15μs  脉冲重复频率  1000Hz 平台平均速度  108m/s 平台平均高度  4874m 中心下视角  0.7854rad 基线长度  1.05m 基线倾角  -0.2358rad  实验中选取的成像区域大小为5km×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高程反演模型可以实现自配准，避免图像配准过程。    (a) BP 算法得到的主辅图像中强散射点位置 (a) The position of intense scatterer in master and slave images obtained by BP algorithm (b) RD算法得到的主辅图像中强散射点位置 (b) The position of intense scatterer in master and slave images obtained by RD algorithm  (c) 16倍插值后BP成像目标沿距离向的分布 (c) Distributions of BP imaging target along the range direction after 16 times interpolation  (d) 16倍插值后RD成像目标沿距离向的分布 (d) Distributions of RD imaging target along the range direction after 16 times interpolation 图  6  强散射目标点成像结果 Fig. 6 The focused images of intense scatterer  BP算法得到的主图像幅度图如图  7(a)所示。实验过程中在该成像区域布放的地面控制点位置信息如图 7(a)中红色三角形标注。由于主辅图像已实现配准，因此直接共轭相乘即可得到干涉相位，如图  7(b)所示。根据第二节中的理论分析，基于DEM辅助后向投影模型的InSAR 高程反演方法由于在成像过程中引入了外源DEM，可直接得到去地形相位后的干涉相位，该干涉相位由辅助DEM成像参考面与真实地形间高程偏差引起，因此残余干涉相位数值较小，如图  7(b)所示，变化范围小于2π/5，避免了相位解缠处理步骤。相位滤波后直接进行高程反演，辅助DEM的高程偏差如图  7(c)所示。将计算的偏差与辅助DEM的高程值相加，最终获得成像区域的高程反演结果，如图  7(d)所示。         (a)  幅度图 (a) SAR amplitude image (b) 干涉相位图 (b) Interferometric phase (c) 辅助DEM偏差 (c) Inversion results of external DEM offset (d) 高程反演结果 (d) DEM inversion result 图  7  机载InSAR数据处理结果图 Fig. 7 Airborne InSAR data processing result  利用实验中布放的部分地面控制点（未参与干涉定标）作为检查点，对传统方法与本文方法高程反演误差进行对比验证，如表  3所示。对比结果表明，传统InSAR成像和干涉处理结果能够达到1m标称高程精度要求（高程精度0.85m），然而利用本文提出的DEM辅助干涉处理算法可以达到更高的处理精度（高程精度0.65m）。由于后向投影成像和干涉处理几何，是一种完全时域的无误差成像算法，在成像信号模型、运动误差补偿模型上理论最为精确，相位保持精度更高；同时在外源DEM辅助下，避免配准和相位解缠过程，减少处理过程中的相位损失。因此，本文方法处理精度优于传统InSAR 成像和干涉处理算法，可以获得高精度的DEM反演结果，且具有简化的干涉处理流程。 表  3  地面检查点处高程反演误差 Tab. 3 The elevation inversion errors of ground detection points 地面检查点  本文方法  (m)  传统InSAR方法(m) 1  -0.5765  -0.8017 2  -0.1310  -0.0824 3  0.7075  0.8840 标准差  0.6519  0.8459  二、丘陵区域实验结果 本文另外选取了丘陵区域的实验数据，进一步验证所提的基于DEM辅助后向投影模型的InSAR 高程反演方法。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高程反演模型的有效性。    (a) 基于后向投影InSAR方法的干涉相位图 (a) Interferometric phase based on backprojection InSAR (b) 传统InSAR方法的干涉相位图 (b) Interferometric phase based on traditional InSAR  (c)  幅度图 (c) SAR amplitude image  (d) 高程反演结果 (d) DEM inversion result 图  8  机载InSAR山区数据处理结果图 Fig. 8 Airborne InSAR data processing results of mountainous area 4．结论 本文提出了一种基于DEM辅助后向投影模型的InSAR高程反演方法。该方法利用BP算法运动补偿精度高和保相性能好的优势，同时在BP成像过程中引入外源DEM辅助，降低干涉条纹的密度，在多数情况下还实现了干涉处理流程的简化——无需图像配准与相位解缠操作。最后本文设计了仿真实验，并对X波段双天线机载InSAR数据进行了处理，通过所提方法中简化的干涉处理流程实现了高精度DEM的反演，验证了该方法的有效性。然而，在基于BP算法的干涉模型中还存在一些问题，例如实际的InSAR系统中天线相位中心存在测量误差，在高分辨成像时需要考虑对方位向配准产生的影响，以及本文模型的高程反演误差理论分析等，这将是今后的研究方向。 参 考 文 献 [1]  Wechsler S P. Uncertainties associated with digital elevation models for hydrologic applications: a review[J]. Hydrology and Earth System Sciences, 2007, 11(4): 1481-1500. doi: 10.5194/hess-11-1481-2007. [2]  Yang  K,  Smith L  C,  Chu V  W, et  al. A  caution  on  the  use  of  surface  digital  elevation  models  to  simulate supraglacial  hydrology  of  the  greenland  ice  sheet[J].  IEEE  Journal  of  Selected  Topics  in  Applied  Earth Observations and Remote Sensing, 2015, 8(11): 5212-5224. doi: 10.1109/JSTARS.2015.2483483. [3]  Shukla  G,  Garg R  D,  Srivastava H  S, et  al. Performance  analysis  of  different  predictive  models  for  crop classification across an aridic to ustic area of indian states[J]. Geocarto International, 2016, 33(3): 240-259. doi: 10.1080/10106049.2016.1240721. [4]  Rufino G, Moccia A, and Esposito S. DEM generation by means of ERS tandem data[J]. IEEE Transactions on Geoscience and Remote Sensing, 1998, 36 (6): 1905-1912. doi: 10.1109/36.729362. [5]  靳国旺,  徐青,  张红敏.  合成孔径雷达干涉测量[M]. 国防工业出版社, 2014. JIN  Guowang,  XU  Qing, ZHANG  Hongmin.  Synthetic  aperture  radar  interferometry[M].  National  Defense Industry Press, 2014. [6]  Seymour M, and Cumming I. InSAR terrain height estimation using low-quality sparse DEMs[C]. Processing of 3rd ERS Symposium, Florence, Italy, 17-21 March 1997, pp. 421-426.  [7]  LIAO Mingsheng,  WANG Teng,  LU Lijun, et  al.  Reconstruction  of  DEMs from ERS-1/2  Tandem data  in mountainous area facilitated by SRTM data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2007, 45 (7): 2325-2335. doi:10.1109/TGRS.2007.896546. [8]  刘辉,靳国旺,张红敏,等.DEM 辅助的山区InSAR相位解缠[J].测绘科学技术学报, 2017, 34(02): 215-220. LIU  Hui,  JIN  Guowang,  ZHANG  Hongmin,  et  al.  Phase unwrapping  assisted  by  DEM  of  InSAR  for mountainous  terrain[J].  Journal  of  Geomatics  Science  and  Technology,  2017,  34(02):  215-220.  doi: 10.3969/j.issn.1673-6338.2017.02.019. [9]  Farr T G, Rosen P A, Caro E, et al. The shuttle radar topography mission[J]. Reviews of Geophysics, 2007, 45 (2): RG2004. doi:10.1029/2005RG000183. [10] Fujisada  H,  Bailey G  B,  Kelly G  G, et  al.  ASTER  DEM  performance[J]. IEEE  Transactions  on Geoscience and Remote Sensing, 2005, 43 (12): 2707-2714. doi: 10.1109/TGRS.2005.847924. [11] Crippen R, Buckley S, Agram P, et al. NASADEM global elevation model: methods and progress[J]. ISPRS - International  Archives  of  the  Photogrammetry,  Remote  Sensing  and  Spatial  Information  Sciences, XLI-B4:125-128. doi:10.5194/isprsarchives-XLI-B4-125-2016. [12] LIN  Jianhe,  LV  Xiaolei,  and  LI  Rui.  Automatic  registered  back-projection  approach  based  on  object orientation for airborne repeat-track interferometric SAR[J]. IET Radar, Sonar and Navigation, 2018, 12 (9): 1066-1076. doi:10.1049/iet-rsn.2018.5053. [13] Duersch M I. Backprojection for synthetic aperture radar[D]. [Ph.D. dissertation], Brigham Young University, 2013. [14] 韦顺军,  师君,  张晓玲,  等.  基于曲面投影的毫米波InSAR  数据成像方法[J]. 雷达学报,  2015,  4(1): 49-59. Wei Shun-jun, Shi Jun, Zhang Xiao-ling, et al.. Millimeter-wave interferometric synthetic aperture radar data imaging based on terrain surface projection[J]. Journal of Radars, 2015, 4(1): 49-59. doi: 10.12000/JR14137. [15] Zebker  H. User-friendly  InSAR  products – do  we  need  range-Doppler?[C]. 12th  European  Conference  on Synthetic Aperture Radar, Aachen, Germany, 4-7 June 2018, pp. 1-4. [16] 师君,  马龙,  韦顺军,  等.  基于导航数据的Ka波段InSAR 成像处理与分析[J].  雷达学报,  2014, 3(01):19-27. SHI Jun, MA Long, WEI Shunjun, et al. Ka-band InSAR imaging and analysis based on IMU data[J]. Journal of Radars, 2014, 3(01):19-27. doi: 10.3724/SP.J.1300.2014.13142. [17] 潘舟浩,  李道京,  刘波,  等.  基于BP算法和时变基线的机载InSAR数据处理方法研究[J]. 电子与信息学报, 2014, 36(7):1585-1591. PAN Zhouhao, LI Daojing, LIU Bo, et al.. Processing of the airborne InSAR data based on the BP algorithm and  the  time-varying  baseline[J]. Journal  of  Electronics & Information Technology, 2014,  36(7):1585-1591. doi: 10.3724/SP.J.1146.2013.00715. [18] 韦顺军,  张晓玲,  师君,  等.  同空间后向投影InSAR 成像及干涉相位提取方法[J].  宇航学报,  2015, 36(03):336-343. WEI Shunjun, ZAHGN Xiaoling, SHI Jun, et al.. InSAR imaging and interferogram extraction-based the same space  back-projection[J].  Journal  of  Astronautics,  2015,  36(03):336-343.  doi: 10.3873/j.issn.1000-1328.2015.03.013. [19] Soumekh M. Synthetic aperture radar signal processing with MATLAB algorithms. New York: Wiley. 1999. [20] 王超,  张红,  刘智.  星载合成孔径雷达干涉测量[M]. 科学出版社, 2002. WANG Chao, ZHANG Hong, LIU Zhi. Spaceborne synthetic aperture radar interferometry[M]. Science Press, 2002.  作 者 简 介 胡晓宁(1994-)，女，山东临沂人，2016年获得北京理工大学学士学位，现于中国科学院空天信息创新研究院攻读博士学位。研究方向为InSAR信号处理算法。 E-mail: huxiaoning16@mails.ucas.ac.cn    汪丙南(1984-)，男，安徽安庆人。2011年在中国科学院电子学研究所获博 士学位，现为中国科学院空天信息创新研究院副研究员，硕士生导师。主要研究方向为干涉合成孔径雷达信号仿真和处理方法。 E-mail: wbn@mail.ie.ac.cn   向茂生(1964-)，男，湖南人。1999年在中国科学院遥感应用研究所获博士学位，现为中国科学院空天信息创新研究院研究员，博士生导师。主要研究双天线干涉、多基线干涉、极化干涉、阵列天线干涉等理论与方法，以及干涉  SAR 面向高精度测绘、复杂地物定位与识别、组合导航等的应用技术。 E-mail: xms@mail.ie.ac.cn   王钟斌（1994-），男，甘肃定西人，2017年获得重庆大学学士学位，现于中国科学院空天信息创新研究院攻读博士学位。研究方向为干涉合成孔径雷达技术。 E-mail: wangzhongbin17@mails.ucas.ac.cn