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

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

胡晓宁, 汪丙南, 向茂生, 等. 基于DEM辅助后向投影模型的InSAR高程反演方法[J]. 雷达学报, 待出版. doi:  10.12000/JR20144
引用本文: 胡晓宁, 汪丙南, 向茂生, 等. 基于DEM辅助后向投影模型的InSAR高程反演方法[J]. 雷达学报, 待出版. doi:  10.12000/JR20144
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
Citation: 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

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

doi: 10.12000/JR20144
基金项目: 国家自然科学基金(62073306),中国科学院青年创新促进会项目资助,国家重点研发计划资助(2017YFC0822400)
详细信息
    作者简介:

    胡晓宁(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

    通讯作者:

    汪丙南 wbn@mail.ie.ac.cn

  • 责任主编:靳国旺 Corresponding Editor: JIN Guowang
  • 中图分类号: TN958

InSAR Elevation Inversion Method Based on Backprojection Model with External DEM

Funds: The National Natural Science Foundation of China (62073306), Youth Innovation Promotion Association CAS, National Key R&D Program of China (2017YFC0822400)
More Information
图(8) / 表 (3)
计量
  • 文章访问数:  36
  • HTML全文浏览量:  3
  • PDF下载量:  6
  • 被引次数: 0
出版历程
  • 收稿日期:  2020-11-28
  • 修回日期:  2021-03-29
  • 网络出版日期:  2021-04-08

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

doi: 10.12000/JR20144
    基金项目:  国家自然科学基金(62073306),中国科学院青年创新促进会项目资助,国家重点研发计划资助(2017YFC0822400)
    作者简介:

    胡晓宁(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

    通讯作者: 汪丙南 wbn@mail.ie.ac.cn
  • 责任主编:靳国旺 Corresponding Editor: JIN Guowang
  • 中图分类号: TN958

摘要: 利用干涉合成孔径雷达(InSAR)技术获取数字高程模型(DEM)时,在地形起伏剧烈区域,干涉条纹十分密集,增加了相位解缠的难度,影响相位展开和高程反演的精度。为了解决该问题,该文提出了一种基于DEM辅助后向投影模型的InSAR高程反演方法。该方法可以在统一的后向投影成像空间中实现成像和InSAR高程反演,并且引入外源DEM作为辅助信息,去除大部分地形相位,有效地降低了干涉条纹的密度,减少了干涉相位的缠绕。此外,该方法在多数情况下可以避免图像配准和相位解缠过程,简化了传统InSAR的处理流程,并且可以实现高精度的高程反演。通过仿真实验和X波段机载双天线InSAR数据的处理验证了该方法的有效性。

注释:
1)  责任主编:靳国旺 Corresponding Editor: JIN Guowang

English Abstract

胡晓宁, 汪丙南, 向茂生, 等. 基于DEM辅助后向投影模型的InSAR高程反演方法[J]. 雷达学报, 待出版. doi:  10.12000/JR20144
引用本文: 胡晓宁, 汪丙南, 向茂生, 等. 基于DEM辅助后向投影模型的InSAR高程反演方法[J]. 雷达学报, 待出版. doi:  10.12000/JR20144
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
Citation: 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
    • 数字高程模型(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系统获取的实测数据处理结果证明了算法的有效性。

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

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

      主要包括以下步骤:

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

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

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

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

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

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

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

      图  1  基于BP成像算法的干涉模型

      Figure 1.  The interferometric model based on BP imaging algorithm

      图1z轴为高度向,$y $轴指向雷达平台的飞行方向,x轴由右手坐标系确定。$ {A}_{1} $$ {A}_{2} $分别代表主、副天线相位中心的位置,下视角为$ \varphi $,基线长度$ B $,基线与水平方向之间的夹角为$ \alpha $。由于目标$ P $的真实高度与辅助DEM之间存在高度偏差,主副天线将$ P $分别投影到$ {Q}_{1} $, $ {Q}_{2} $,而成像网格中对应的离散采样点为C。天线相位中心$ {A}_{1} $的高程为$ H $。点C, $ {Q}_{1} $$ {Q}_{2} $ 的高程值分别为 $ {z}_{c} $, $ {z}_{{Q}_{1}} $$ {z}_{{Q}_{2}} $。投影曲面上的采样点C与目标点$ P $的真实位置之间的高程偏差为$ \Delta h $, $ \theta $$ \overline {{\mathrm{A}}_{1}{\mathrm{Q}}_{1}} $与垂直方向的夹角。

      首先分析图像对偏移量$ \overline {{\mathrm{Q}}_{1}{\mathrm{Q}}_{2}} $xy平面上的投影$ {\overline {{\mathrm{Q}}_{1}{\mathrm{Q}}_{2}}}_{\mathrm{P}\mathrm{r}\mathrm{o}\mathrm{j}} $。根据图1中的几何关系,$ \overline {{\mathrm{A}}_{1}\mathrm{P}}=\overline {{\mathrm{A}}_{1}{\mathrm{Q}}_{1}}, $$ \overline {{\mathrm{A}}_{2}\mathrm{P}}=\overline {{\mathrm{A}}_{2}{\mathrm{Q}}_{2}} $,可得到式(1)和式(2)。

      $$\begin{split} {\overline {{{\rm{Q}}_{\rm{1}}}{{\rm{Q}}_2}} _{{\rm{Proj}}}} =& \sqrt {{{\overline {{{\rm{A}}_{\rm{1}}}{\rm{P}}} }^2} - {{\left({H - {z_{{Q_{\rm{1}}}}}} \right)}^2}} - \bigg(B\cos \alpha \\ & + \sqrt {{{\overline {{{\rm{A}}_{\rm{2}}}{\rm{P}}} }^2} - {{\left({H + B\sin \alpha - {z_{{Q_{\rm{2}}}}}} \right)}^2}} \bigg)\;\; \end{split}$$ (1)
      $$\Delta h = H - {z_{\rm{C}}} - \overline {{{\rm{A}}_{\rm{1}}}{\rm{P}}} \sin \varphi $$ (2)

      根据余弦定理可得

      $$\overline {{{\rm{A}}_{\rm{2}}}{\rm{P}}} = \sqrt {{B^2} + {{\overline {{{\rm{A}}_{\rm{1}}}{\rm{P}}} }^2} - 2B \cdot \overline {{{\rm{A}}_{\rm{1}}}{\rm{P}}} \cdot {\rm{cos}}\left({\alpha + \varphi } \right)} $$ (3)

      由式(1)—式(3),计算$ {\overline {{\mathrm{Q}}_{1}{\mathrm{Q}}_{2}}}_{\mathrm{P}\mathrm{r}\mathrm{o}\mathrm{j}} $$ \Delta h $的泰勒展开形式,得到图像对偏移量的表达式如式(4),其中$ {\overline {\mathrm{O}{\mathrm{Q}}_{2}}}_{\mathrm{P}\mathrm{r}\mathrm{o}\mathrm{j}} $$ \overline {\mathrm{O}{\mathrm{Q}}_{2}} $xy平面上的投影。

      $${\overline {{{\rm{Q}}_{\rm{1}}}{{\rm{Q}}_{\rm{2}}}} _{{\rm{Proj}}}} \approx - \frac{{B\sin (\alpha + \varphi )}}{{{{\overline {{\rm{O}}{{\rm{Q}}_{\rm{2}}}} }_{{\rm{Proj}}}}\cos \varphi }}\Delta h$$ (4)

      由式(4)可以看出,图像对偏移量$ {\overline {{\mathrm{Q}}_{1}{\mathrm{Q}}_{2}}}_{\mathrm{P}\mathrm{r}\mathrm{o}\mathrm{j}} $与辅助DEM的高程误差有关,针对目前可以获取的外源DEM精度,其高程误差引起的图像对偏移量少于一个分辨单元的大小。以3.2节的机载InSAR的参数为例,根据式(4)可以计算当投影曲面高程误差为10 m时,图像对偏移量大概为10–3~3 m,即只有当$ \Delta h $达到百米量级时,图像对偏移量才可能超过一个分辨单元。因此可以认为经离散采样后点目标P位于BP成像网格的同一采样点处。也就是说,基于BP算法的InSAR模型在成像过程中可以实现自配准,无需后续图像配准操作,主辅图像直接共轭相乘即可得到干涉相位,简化了干涉处理的流程。

      基于上述结论,下面推导基于DEM辅助后向投影模型的InSAR高程反演公式。在BP成像过程中,主图像和辅图像理论上需要补偿的多普勒相位项分别为$ -4{\pi }\overline {{\mathrm{A}}_{1}\mathrm{P}}/\lambda $$ -4{\pi }\overline {{\mathrm{A}}_{2}\mathrm{P}}/\lambda $。但在实际的成像过程中,补偿的多普勒相位项分别为$ 4{\pi }\overline {{\mathrm{A}}_{1}\mathrm{C}}/\lambda $$ 4{\pi }\overline {{\mathrm{A}}_{2}\mathrm{C}}/\lambda $,因此由于BP成像网格化导致主辅图像对之间存在式(5)所示的相位差。

      $$\begin{split} \Delta \varphi =& \left[ { - \frac{{4\pi }}{\lambda }\left({\overline {{{\rm{A}}_{\rm{1}}}{\rm{P}}} {\rm{ - }}\overline {{{\rm{A}}_{\rm{1}}}{\rm{C}}} } \right)} \right] - \left[ { - \frac{{4\pi }}{\lambda }\left({\overline {{{\rm{A}}_{\rm{2}}}{\rm{P}}} {\rm{ - }}\overline {{{\rm{A}}_{\rm{2}}}{\rm{C}}} } \right)} \right] \\ = & - \frac{{4\pi }}{\lambda }\left[ {\left({\overline {{{\rm{A}}_{\rm{1}}}{\rm{P}}} {\rm{ - }}\overline {{{\rm{A}}_{\rm{2}}}{\rm{P}}} } \right){\rm{ - }}\left({\overline {{{\rm{A}}_{\rm{1}}}{{\rm{Q}}_{\rm{1}}}} {\rm{ - }}\overline {{{\rm{A}}_{\rm{2}}}{{\rm{Q}}_{\rm{1}}}} } \right)} \right] \\ & - \frac{{4\pi }}{\lambda }\left[ {\left({\overline {{{\rm{A}}_{\rm{1}}}{{\rm{Q}}_{\rm{1}}}} {\rm{ - }}\overline {{{\rm{A}}_{\rm{2}}}{{\rm{Q}}_{\rm{1}}}} } \right){\rm{ - }}\left({\overline {{{\rm{A}}_{\rm{1}}}{\rm{C}}} {\rm{ - }}\overline {{{\rm{A}}_{\rm{2}}}{\rm{C}}} } \right)} \right] \\ =& (\Delta {\varphi _{\rm{P}}} - \Delta {\varphi _{{{\rm{Q}}_1}}}) + (\Delta {\varphi _{{{\rm{Q}}_{\rm{1}}}}} - \Delta {\varphi _{\rm{C}}})\\[-12pt] \end{split} $$ (5)

      在传统InSAR的干涉相位模型中,P和${{\rm{Q}}}_{1}$两点间的干涉相位差与其高程差$ \Delta {h}_{\mathrm{P}{\mathrm{Q}}_{1}} $有关,根据高程相位的表达式[20],该干涉相位差如式(6)所示。

      $$\Delta {\varphi _{\rm{P}}} - \Delta {\varphi _{{{\rm{Q}}_{\rm{1}}}}} = - \frac{{4\pi }}{\lambda }\frac{{B\cos \left({\theta - \alpha } \right)}}{{\overline {{{\rm{A}}_{\rm{1}}}{\rm{P}}} \sin \theta }}\Delta {h_{{\rm{P}}{{\rm{Q}}_{\rm{1}}}}}$$ (6)

      根据图1中的几何关系,${{\rm{Q}}}_{1}$和C点之间的干涉相位差可以表示为$ \Delta {r}_{\mathrm{C}{\mathrm{Q}}_{1}} $引起的平地相位和$ \Delta {h}_{\mathrm{C}{\mathrm{Q}}_{1}} $引起的高程相位之和的形式,如式(7)所示。

      $$\begin{split} \Delta {\varphi _{{Q_{\rm{1}}}}} - \Delta {\varphi _C} =& - \frac{{4\pi }}{\lambda }\frac{{B\cos \left({\theta - \alpha } \right)}}{{\overline {{{\rm{A}}_{\rm{1}}}{\rm{P}}} \sin \theta }}\Delta {h_{{\rm{C}}{{\rm{Q}}_{\rm{1}}}}}\\ & - \frac{{4\pi B\cos (\theta - \alpha )}}{{\lambda \overline {{{\rm{A}}_{\rm{1}}}{{\rm{Q}}_{\rm{1}}}} \tan \theta }}\Delta {r_{{\rm{C}}{{\rm{Q}}_{\rm{1}}}}} \end{split}$$ (7)

      其中,$ \Delta {h}_{\mathrm{C}{\mathrm{Q}}_{1}} $$ \Delta {r}_{\mathrm{C}{\mathrm{Q}}_{1}} $分别是${{\rm{Q}}}_{1}$和C两点之间的高程差与到雷达平台的斜距差。

      图1中,$ \Delta h=\Delta {h}_{\mathrm{P}{\mathrm{Q}}_{1}}+\Delta {h}_{\mathrm{C}{\mathrm{Q}}_{1}} $,所以干涉相位可以写为

      $$\begin{split} \Delta \varphi = & - \frac{{4\pi }}{\lambda }\frac{{B\cos \left({\theta - \alpha } \right)}}{{\overline {{{\rm{A}}_{\rm{1}}}{\rm{P}}} \sin \theta }}\Delta h - \frac{{4\pi B\cos (\theta - \alpha )}}{{\lambda \overline {{{\rm{A}}_{\rm{1}}}{{\rm{Q}}_{\rm{1}}}} \tan \theta }}\Delta {r_{{\rm{C}}{{\rm{Q}}_{\rm{1}}}}} \\ =& \Delta {\varphi _h} + \Delta {\varphi _r}\\[-12pt] \end{split} $$ (8)

      当成像投影面为平面时,即为$ \Delta {h}_{\mathrm{C}{\mathrm{Q}}_{1}}=0 $的特例,式(8)仍然成立。

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

      $$\Delta \varphi \approx - \frac{{4\pi }}{\lambda } \cdot \frac{{B\cos \left({\theta - \alpha } \right)}}{{\overline {{{\rm{A}}_{\rm{1}}}{\rm{P}}} \cdot \sin \theta }} \cdot \Delta h$$ (9)

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

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

      图  2  高程模糊度几何示意图

      Figure 2.  Geometric description of the height of ambiguity

      $$\Delta {h_{2\pi }} = \frac{{\lambda \cdot \overline {{{\rm{A}}_{\rm{1}}}{\rm{P}}} \cdot \sin \theta }}{{2B\cos \left({\theta - \alpha } \right)}}\hspace{45pt}$$ (10)
      $$\begin{split} {h_{2\pi }} =& {h_{\rm{2}}} - {h_1} \\ =& (\Delta {h_2} - \Delta {h_1}) + ({h_{{\rm{e}}2}} - {h_{{\rm{e}}1}}) \\ =& \frac{{\lambda \cdot \overline {{{\rm{A}}_{\rm{1}}}{\rm{P}}} \cdot \sin \theta }}{{2B\cos \left({\theta - \alpha } \right)}} + ({h_{{\rm{e}}2}} - {h_{{\rm{e}}1}}) \end{split} $$ (11)

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

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

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

      图  3  InSAR处理流程

      Figure 3.  InSAR processing chain

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

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

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

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

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

      表 1  仿真参数

      Table 1.  Simulation parameters

      成像仿真参数数值
      载频(GHz)9.6
      距离向带宽(MHz)100
      距离向采样率(MHz)120
      脉冲宽度(μs)3.7
      脉冲重复频率(Hz)300
      平台平均速度(m/s)113.5
      平台平均高度(m)3286.5
      中心下视角(rad)0.8727
      基线长度(m)2.189
      基线倾角(rad)0

      图  4  仿真实验结果图

      Figure 4.  Simulation results

      首先利用后向投影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高程反演方法在解决相位模糊问题的同时简化了干涉处理流程,因此具有更大的优势。

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

      图  5  两种方法的高程反演误差

      Figure 5.  The elevation inversion errors of two methods

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

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

    • (1) 平地区域实验结果

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

      表 2  X波段机载InSAR系统参数

      Table 2.  X-band airborne InSAR system parameters

      机载InSAR参数数值
      载频(GHz)9.6
      距离向带宽(MHz)300
      距离向采样率(MHz)500
      脉冲宽度(μs)15
      脉冲重复频率(Hz)1000
      平台平均速度(m/s)108
      平台平均高度(m)4874
      中心下视角(rad)0.7854
      基线长度(m)1.05
      基线倾角(rad)–0.2358

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

      图  6  强散射目标点成像结果

      Figure 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数据处理结果图

      Figure 7.  Airborne InSAR data processing results

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

      表 3  地面检查点处高程反演误差

      Table 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
    • (2) 丘陵区域实验结果

      本文另外选取了丘陵区域的实验数据,进一步验证所提的基于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高程反演模型的有效性。

      图  8  机载InSAR山区数据处理结果图

      Figure 8.  Airborne InSAR data processing results of mountainous area

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

参考文献 (20)

目录

    /

    返回文章
    返回