一种基于高精度DEM提高InSAR技术监测地表形变精度的方法

摘要

一种基于高精度DEM提高InSAR技术监测地表形变精度的方法，它有五大步骤：步骤一、由雷达数据生成干涉纹图；步骤二、差分干涉相位图的生成；步骤三、误差相位的构成及特征分析；步骤四、误差相位最优函数校正模型的建立；步骤五、根据步骤二和步骤四的结果恢复监测区的地表形变信息。该方法通过提取研究区不同区域的误差相位和高程值或误差相位、高程以及沿距离/方位向的坐标值，基于最小二乘法，分别建立相应区域误差相位的最优函数校正模型，最后，将模拟的误差相位从差分干涉图中去除，进而恢复监测区沿雷达视线向的形变信息。本发明在星载合成孔径雷达监测地表形变技术应用领域具有实用价值和广阔的应用前景。

说明书

技术领域

本发明涉及一种基于高精度数字高程模型（DEM）提高合成孔径雷达干 涉测量（InSAR）技术监测地表形变精度的方法，属于星载合成孔径雷达监 测地表形变技术领域，它适用于基于高分辨率雷达数据和外部高精度DEM， 根据误差相位与研究区高程之间的某种函数关系，建立研究区的区域最优误 差相位校正模型，恢复研究区地面目标沿雷达视线向的形变信息，从而进一 步提高InSAR技术监测地表形变的能力。

背景技术

InSAR技术监测地表形变的核心是通过获取雷达与地面目标之间的相位 信息恢复地面监测目标的几何及形变特征，对地表微小形变具有很强的敏感 性。由于雷达两次观测条件的变化，包括雷达轨道偏差、地物后向散射特性 变化以及噪声的影响，使得两次获取的地面目标的雷达信号相干性降低，难 以完成空间离散分布的雷达目标的相位解缠。此外，雷达卫星的基线误差、 引入的外部DEM数据不准确以及雷达卫星两次观测时刻大气波动的影响等， 这些因素都会降低InSAR技术监测地表形变结果的精度。针对这一问题，当 前常用的解决方法有两种：一种是传统合成孔径雷达差分干涉技术 （D-InSAR）；一种是基于具有稳定散射特性的相干点目标时间序列分析技 术。

传统D-InSAR技术的基本思路是，在保证主辅影像配准误差在1/8个像 元以下，并对主辅影像进行带通滤波处理后，通过主辅影像共轭相乘及差分 处理去除平地相位和地形相位，最终得到的差分干涉图中的误差相位主要是 由雷达数据两次获取期间较大的大气波动引起的。在可获取外部大气数据的 情况下，将大气层中的水汽含量转换为大气相位，将其从差分相位图中消除， 进而得到监测区地表的形变信息。这种方法假设轨道基线及外部DEM数据 完全精确，数据在处理过程中产生的误差忽略不计，单纯依靠外部地面气象 参数进行水汽延迟校正，但现有的大气水汽监测网络以及物理模型并不能对 大气中的水汽含量沿微波信号的路径进行足够准确的外推；基于具有稳定散 射特性的相干点目标时间序列分析技术的基本思路是，基于具有稳定散射特 性的相干点目标，根据其干涉相位的构成及各相位分量的时空特征，对相干 目标的干涉相位进行时间序列分析，估算DEM误差、大气波动与轨道误差 以及噪声等，逐步将这些非形变相位分离，从而准确获取地表的形变信息。 这些方法包括永久散射体干涉测量技术（PSI）、短基线干涉测量法（SBAS）、 点目标干涉测量分析（IPTA）等，它们对雷达数据量的要求较高，一般需要 25景甚至更多的SAR影像，当数据量较少且监测区非线性形变梯度较大时， 应用中存在一定的局限。

发明内容

本发明的目的是，针对现有方法没有考虑轨道基线及外部DEM数据不 准确引起的误差，使得InSAR技术监测地表形变结果的精度不高，甚至有时 难以恢复短时间段内地表形变信息的问题，提供了一种基于高精度DEM提 高InSAR技术监测地表形变精度的方法，它是一种通过建立差分干涉图中误 差相位与研究区高程之间的最优函数校正模型去除误差相位提高InSAR监 测地表形变精度的方法。通过推导、分析发现，差分干涉图中的误差相位主 要由两部分组成：一是轨道基线和外部DEM数据不准确引起的误差相位； 一是大气波动，即大气湿分量变化引起的延迟相位，两者都与研究区的高程 密切相关。另外，在地形起伏较大、降水丰沛的高山峡谷区的局部区域，大 气延迟相位也是地理位置的函数。于是根据研究区误差相位的分布特征不同， 分别提取相应区域内的误差相位、高程信息以及沿距离向/方位向的坐标，基 于最小二乘法，对误差相位与高程、沿距离向/方位向的坐标值进行函数拟合， 在保证函数模型简单、残差平方和最小的前提下，建立相应区域误差相位的 最优多项式校正模型。最后，将模型模拟的误差相位从差分干涉图中去除， 进而较为准确地恢复出形变区域沿雷达视线向的形变信息。

技术解决方案

本发明的技术解决方案见附图1，包括以下五个步骤：

（1）由雷达数据生成干涉纹图；

（2）差分干涉相位图的生成；

（3）误差相位的构成及特征分析；

（4）误差相位最优函数校正模型的建立；

（5）根据步骤二和步骤四的结果恢复监测区的地表形变信息。

本发明一种基于高精度DEM提高InSAR技术监测地表形变精度的方法，该 方法具体步骤如下：

步骤一：由雷达数据生成干涉纹图

通过星载SAR传感器获取的雷达数据，基于相干性进行主辅影像选择。将 主辅影像进行精确配准，根据配准多项式将辅影像与主影像的配准误差控制在 1/8个像元以下，并对两者进行带通滤波处理，然后将它们逐像素共轭相乘生成 干涉纹图。

步骤二：差分干涉相位图的生成

基于卫星轨道参数及干涉纹图的条纹频率，分别估算平行基线和垂直基线； 根据SAR影像的坐标系统及投影等信息，将外部高精度DEM转换到SAR坐标 系下，结合卫星轨道信息模拟研究区的地形相位、平地相位，将其从干涉纹图 中去除，得到差分干涉相位图。

步骤三：误差相位的构成及特征分析

经过差分处理并对差分干涉相位图进行相位解缠，解缠后的差分干涉相位 图中任一像素x的差分干涉相位△φ_diff(x)可由形变相位△φ_def(x)和误差相位 △φ_error(x)两部分组成，其中误差相位△φ_error(x)可由如下公式表示。

△φ_error(x)=△φ_ε-h(x)+△φ_b(x)+△φ_atm(x)+△w(x)   （1）

其中，△φ_ε-h(x)为垂直基线及外部DEM数据不准确造成的地形误差相位， △φ_b(x)为平行基线不准确引起的平地误差相位，△φ_atm(x)为大气湿延迟引起的误 差相位，△w(x)为随机的噪声误差相位。

垂直基线及外部DEM数据不准确造成的地形误差相位的计算公式如下所 示。

$\Delta {\phi }_{ϵ-h}\left(x\right)=-\frac{4\pi }{\lambda R\sin {\theta }_0}(\mathrm{h\Delta }{B}_{\perp }+\mathrm{\Delta h}{B}_{\perp })---\left(2\right)$

其中，h为雷达影像覆盖区任一点x处的高程值，△h为由于外部的DEM数 据不准确，雷达影像覆盖区相应点x处的高程误差值，B_⊥为两次雷达数据获取 时刻雷达卫星之间的垂直基线长度，△B_⊥为垂直基线误差，R为雷达影像覆盖区 域中任一点x与第一副天线之间的距离，θ₀为雷达卫星的入射角，λ为雷达波的 波长。

平行基线不准确引起的平地误差相位的计算公式如下所示。

$\Delta {\phi }_b\left(x\right)=-\frac{4\pi }{\lambda \sin ({\theta }_0-\alpha )}\times \frac{[\frac{B_{//}}{\sin ({\theta }_0-\alpha )}-R\sin (\arccos \frac{R^2+H^2-R^{\prime 2}}{2\mathrm{RH}}-\alpha )]}{\sqrt{\frac{{B_{//}}^2}{{\sin }^2({\theta }_0-\alpha )}+R^2-\frac{2B_{//}R}{\sin ({\theta }_0-\alpha )}\sin (\arccos \frac{R^2+H^2-R^{\prime 2}}{2\mathrm{RH}}-\alpha )}}\Delta B_{//}---\left(3\right)$

其中，B_//为雷达数据两次获取时刻雷达卫星之间的平行基线长度，△B_//为平 行基线误差，H为雷达卫星至地心的距离，R为雷达影像覆盖区域中任一点x与 第一副天线之间的距离，R'为雷达影像覆盖范围内地面任一点至地心的距离，α 为基线与水平面的夹角，θ₀为雷达卫星的入射角，λ为雷达波的波长。

由于雷达卫星在两次获取雷达影像时，雷达天线之间的距离B近似等于雷达 影像覆盖区域中任一点x到第一副天线之间的距离R与B和R之间夹角余弦的 乘积，其计算公式可表达如下所示

$B=\frac{B_{//}}{\sin ({\theta }_0-\alpha )}\approx R\cos (\frac{\pi }{2}+\alpha -\arccos \frac{R^2+H^2-R^{\prime 2}}{2\mathrm{RH}})=R\sin (\arccos \frac{R^2+H^2-R^{\prime 2}}{2\mathrm{RH}}-\alpha )---\left(4\right)$

于是，△φ_b(x)≈0

大气湿延迟引起的误差相位△φ_atm(x)，在空间域呈低频变化，但在空气湿度 变化较大的高山峡谷区（如水库库区），大气湿分量在空间域的分布特征不尽相 同，不仅与研究区高程h有关，在某些局部区域还与研究区的位置x密切相关。

雷达数据处理过程中产生的随机噪声误差相位△w(x)，通过滤波处理可基本 消除。

可见，差分干涉相位图中的误差相位主要由轨道基线和外部DEM数据不准 确引起的误差相位和大气波动引起的误差相位两部分构成，因此，误差相位可 进一步用下列公式表达。

$\left(\begin{array}{c}\Delta {\phi }_{\mathrm{error}}\left(x\right)=\Delta {\phi }_{ϵ-h}\left(x\right)+\Delta {\phi }_{\mathrm{atm}}\left(x\right)\\ =-\frac{4\pi }{\lambda R\sin {\theta }_0}(\mathrm{h\Delta }{B}_{\perp }+\mathrm{\Delta h}{B}_{\perp })+\Delta {\phi }_{\mathrm{atm}}(x,h)\end{array}\right)---\left(5\right)$

步骤四：误差相位最优函数校正模型的建立

差分干涉相位图中的误差相位与高程密切相关，通过提取相应区域（除形 变区）的误差相位△φ_error(x)和高程值h，基于最小二乘法，尝试采用如下多种类 型函数，分别进行曲线拟合。

多项式函数：△φ_error(x)=a₀+a₁×h+a₂×h²+a₃×h³+…   （6）

指数函数：△φ_error(x)=a×e^b×h   （7）

高斯函数：$\Delta {\phi }_{\mathrm{error}}\left(x\right)=a\times e^{-{\left(\frac{h-b}{c}\right)}^2}---\left(8\right)$

幂函数：△φ_error(x)=a×h^b+c   （9）

其中，a₀，a₁，a₂，a₃，…，a，b，c是各拟合函数的系数。

通过多次试验比较发现，二次多项式函数刻画的误差相位校正模型，在研 究区高程范围内没有明显的转折点，残差相位相对较小，且随着多项式最高次 数继续增加，其最高次数的系数小于10^-9，这使得拟合多项式的误差平方和减 小的程度不明显，甚至误差平方和基本不变。其他类型的函数模型，指数函数 在研究区高程范围内有明显的转折点，不符合误差相位随高程递减变化的规律； 高斯函数模型不但有明显转折点，而且误差相位的残留误差平方和较大；幂函 数模型描述的误差相位随高程的变化趋势以及残差相位的最小平方和与多项式 函数模型基本一致，但其表达方式不如多项式简单。因此，误差相位的最优函 数校正模型可以用关于高程h的二次多项式函数较为准确地表达。

△φ_error(x)=a₀+a₁×h+a₂×h²   （10）

其中，a₀，a₁，a₂是二次多项式函数的系数。

在研究区的某些局部区域，虽然大气湿分量在空间域呈低频变化，但在地 形起伏较大的高山峡谷区，空气湿度变化较大，误差相位的变化特征不仅与高 程有关，还与研究区的位置密切相关，此局部区域误差相位的最优函数校正模 型应用如下关于高程h的二次和关于位置x的一次多项式函数进行表达。

△φ_error(x)=a₀+a₁×h+a₂×h²+b₁×x+b₂×x×h   （11）

其中，a₀，a₁，a₂，b₁，b₂是二次多项式函数的系数。

步骤五：根据步骤二和步骤四的结果恢复监测区的地表形变信息

根据建立的不同区域的最优函数校正模型，模拟出研究区整个区域的误差 相位分布图，将其从差分干涉图中去除。最后，由卫星轨道参数将去除误差相 位后的相位信息转换为沿雷达视线向的形变信息，进而恢复出地表形变区的范 围及形变量。

优点及功效：本发明一种基于高精度DEM提高InSAR技术监测地表形 变精度的方法，其优点是：

(1)本方法基于高精度DEM，利用传统合成孔径雷达差分干涉（D-InSAR） 技术，根据误差相位与高程之间的最优函数模型消除差分干涉图中的误 差相位，从而提高干涉测量的精度。

(2)本方法在无需订购多景雷达数据进行回归分析，只需两景雷达数据，且 无需获取研究区外部大气数据的情况下，通过建立研究区误差相位关于 高程的最优函数校正模型即可获取研究区误差相位的分布特征。

(3)本方法利用误差相位与高程之间的函数模型消除差分干涉图中的误差 相位，这在地形起伏较大、大气湿延迟现象严重的高山地区极为有效， 可以准确地获取地表短时间内的非线性形变，从而极大地扩大了传统差 分干涉雷达的测量范围。

附图说明

图1为本发明的流程图。

图2为有严重误差相位存在的解缠后的差分干涉相位图。在树坪滑坡体区 域，误差相位将滑坡体的形变信息完全掩盖。其中，等值线上数值代表差分干 涉图中不同位置相位值的大小，单位为弧度；粗虚线范围内为树坪滑坡体，细 实线和细虚线囊括的范围分别代表区域I和区域II。整个研究区中的误差相位与 高程密切相关,且随着高程的增加误差相位值不断减小，但树坪滑坡体区域和区 域II中的误差相位不仅与高程密切相关，还是研究区沿雷达距离向坐标的函数， 区域A为区域I、II的共同区域。Az代表雷达方位向，Rg代表雷达距离向。

图3为由最优函数校正模型模拟的研究区的误差相位分布图，该误差相位 的分布特征与差分干涉图(即图2)中系统误差相位的分布规律基本一致。其中， 等值线上数值代表研究区不同位置误差相位值的大小，单位为弧度。Az代表雷 达方位向，Rg代表雷达距离向。

图4为树坪滑坡体在11天中沿雷达视线向的形变信息图。其中，粗虚线范 围内为树坪滑坡体区域，等值线上数值代表获得的树坪滑坡体不同区域沿雷达 视线向上的形变值，两个形变中心明显的表示出来，单位为毫米；树坪滑坡体 区域外为稳定区域，等值线上数值代表残留的误差值，单位为毫米，此时整个 研究区残留的误差不但其数值较小，且只在局部区域集中存在。Az代表雷达方 位向，Rg代表雷达距离向。

具体实施方式

见附图，为了更好的说明本发明的方法与步骤，利用三峡水库南岸树坪滑 坡区两景分辨率为1m的间隔时间为11天的TerraSAR-X雷达数据与精度皆为 1m的Lidar数字高程模型（DEM）数据为例，进行误差相位去除，恢复树坪滑 坡体滑动区形变信息的试验。

（1）试验所用的设备为图形工作站，规格型号为Dell Precision T7400。

（2）见图1，本发明一种基于高精度DEM提高InSAR技术监测地表形变精度 的方法，该方法具体步骤如下：

步骤一：由雷达数据生成干涉纹图

基于传统差分干涉测量方法（D-InSAR），以2012年1月2日获取的分辨 率为1m的TerraSAR-X雷达数据为主影像，将2012年1月13日获取的同样 分辨率的TerraSAR-X雷达数据与之进行精确配准，在保证主辅影像的配准误 差控制在1/8个像元以下，并对配准后的这两幅影像进行带通滤波处理，最后 将两者逐像素共轭相乘生成干涉纹图。

步骤二：差分干涉相位图的生成

基于卫星轨道参数及干涉纹图的条纹频率，分别进行空间基线的估算，由 空间基线计算干涉图的平地相位；根据SAR影像的坐标系统及投影等信息，对 分辨率为1m的数字高程模型（DEM）数据进行坐标转换，即地理编码，通过 进行SAR图像和外部DEM之间坐标系统的转换，使外部DEM处于SAR坐标 系下，结合卫星轨道信息模拟出研究区正确的地形相位，最后将平地相位、地 形相位从干涉纹图中去除，得到差分干涉相位图。

步骤三：误差相位的构成及特征分析

为准确恢复树坪滑坡体的形变信息，首先以某一稳定点为相位解缠基准点， 采用最小费用流算法对差分干涉图进行相位解缠，解缠后的差分干涉相位图见 图2。

从图2中可以看出，差分干涉图中存在严重的误差相位，树坪滑坡体滑动 区域的形变相位完全被误差相位“淹没”，无法将滑坡体的形变信息和周围稳 定区域的误差相位信息区分开来。由于误差相位主要由轨道基线和外部DEM数 据不准确引起的误差相位和大气波动引起的误差相位两部分构成，差分干涉相 位图中误差相位是研究区高程的函数，且随着高程的增加，误差相位值不断降 低。但II区域及树坪滑坡区，它们的误差相位分布特征基本一致，不仅与高程 密切相关，还是研究区沿雷达距离向坐标的函数。

于是，为恢复整幅差分干涉相位图中误差相位的分布状态，针对I、II两区 域的误差相位分布特征不同，分别建立相应区域误差相位的最优函数校正模型。

步骤四：误差相位最优函数校正模型的建立

在I、II区域中，分别提取高相干点上的误差相位和相应的高程 值h^I、h^II，并提取II区域中相应点处沿距离向的坐标值x。基于最小二乘法，采 用式（10）和式（11）分别建立I、II区域的误差相位校正模型如下所示。

$\Delta {\phi }_{\mathrm{error}}^I=2.48-8.04\times {10}^{-3}\times h^I+4.31\times {10}^{-6}\times {\left(h^I\right)}^2---\left(12\right)$

$\Delta {\phi }_{\mathrm{error}}^{\mathrm{II}}=3.63-6.32\times {10}^{-3}\times h^{\mathrm{II}}+4.01\times {10}^{-8}\times {\left(h^{\mathrm{II}}\right)}^2-4.47\times {10}^{-3}\times x-3.46\times {10}^{-8}\times x\times h^{\mathrm{II}}---\left(13\right)$

它们的残差相位都介于[-2.0，1.0]之间。由于树坪滑坡体区域的误差相位分布特 征与II区域基本一致，则树坪滑坡区的误差相位分布特征也遵循式（13）的表 达规律。

由提取的I、II区域及树坪滑坡体区域高相干点上的误差相位和 相应的高程值h^I、h^II以及II区域和树坪滑坡体区域中相应点处沿距离向的坐标 值x，根据建立的I、II区域误差相位的最优函数校正模型，即式（12）、（13）， 计算整个研究区高相干点上的误差相位值，其中，I、II区域的共同区域A中高 相干点上的误差相位值为由式（12）、（13）分别计算得到的误差相位的平均值。 最后，通过最邻近插值法，将计算得到的整个区域高相干点上的误差相位值进 行插值，进而得到研究区误差相位的分布特征图，见图3。在树坪滑坡体区域， 误差相位的分布特征也形象地表达了出来。

步骤五：根据步骤二和步骤四的结果恢复监测区的地表形变信息

将由最优函数校正模型式（12）、（13）模拟出的研究区的误差相位，从解 缠后的差分干涉相位图中消除，此时，在树坪滑坡体区域残留的相位只是由滑 坡体的滑动形变引起的。最后基于卫星轨道参数将去除误差相位后的相位信息 转换为沿雷达视线向的形变信息，树坪滑坡体在2012年1月2日至2012年1月13 日间的滑动区域及滑动量就准确、鲜明地表达了出来，见图4。