CN111273293A - 一种顾及地形起伏的InSAR残余运动误差估计方法及装置 - Google Patents

Abstract

本发明公开了一种顾及地形起伏的InSAR残余运动误差估计方法及装置，其中方法包括：获取多景SAR影像并进行预处理，计算包含InSAR残余运动误差的地形相位；利用外部DEM数据计算模拟地形相位；将包含InSAR残余运动误差的地形相位与模拟地形相位求差，并对差值进行相位解缠处理，再将解缠差分干涉相位作为观测值；构建有关于InSAR干涉几何参数、DEM误差以及基线误差分量的差分干涉相位参数化模型；针对每个从影像的每一距离行像素，均将差分干涉相位的观测值与差分干涉相位参数化模型进行匹配构建等式，得到观测方程；解算多基线观测方程组中的DEM误差和基线误差分量。本发明能准确地估计残余运动误差。

Description

一种顾及地形起伏的InSAR残余运动误差估计方法及装置

技术领域

本发明属于合成孔径雷达干涉测量技术领域，具体涉及一种顾及地形起伏的InSAR残余 运动误差估计方法及装置。

背景技术

近些年来，合成孔径雷达干涉测量(InSAR)技术已经具备了大范围高精度地形测绘的 能力，比如目前已经发布了覆盖全球的高精度DEM数据，如STRM DEM、TanDEM-X DEM等，并在大范围灾害监测和资源调查等方面得到了广泛的应用。机载SAR系统(如PI-SAR、E-SAR、F-SAR、UAVSAR等)凭借着机动灵活、重访时间短等特点，而且飞行高度低不易 受大气影响，因此为局部地区更高分辨率、更高精度DEM反演带来了极大便利。

众所周知，机载SAR系统在飞行过程中经常由于受到外界随机因素扰动，如气流影响等， 其运动轨迹不可避免地会偏离理想的飞行航迹，机载SAR平台的飞行轨迹通常严重偏离了理 想航线，是高度非线性的，于是SAR数据的聚焦和成像质量受到严重影响。尽管SAR数据 经过了高质量的成像处理，但由于惯导系统的定位精度(1-10cm)是有限的和成像算法是不 精确的，运动补偿后的航迹中仍残留未补偿的运动误差，即RMEs，这是当前机载干涉SAR系 统中最主要的误差源。在干涉过程中，RMEs除了导致InSAR几何发生扭曲外，还会造成方 位向和距离向产生严重的相位波动和偏差，特别是重轨模式中。就SAR平台搭载的惯导系统 定位误差而言，假设水平基线为30m，目标到传感器的斜距为5000m，入射角为45度的情况 下，仅1cm的RME就会在距离向引起约1.7m的测高误差。而未知基线偏差是机载SAR运 动补偿过程中不可避免的，对应的量级约半波长的整数倍，这种测高影响会更大。由此可见， 要利用机载重轨干涉SAR数据时必须要考虑RMEs的估计和改正。

根据RMEs的特征表现，通常可以将RMEs分为两部分，即高频部分和低频部分。

(1)高频RMEs，主要是由惯导系统的定位误差等引起，主要表现为时频方位偏差。对 此，目前被采用较为广泛的一类方法就是Multi-Squint(MS)方法，其基本原理在于不同子 视干涉相位的差值与方位向位置误差存在一定的关系，可以进一步建立不同子视干涉相位的 差值与基线误差的导数之间的函数关系，以此消除干涉图中RMEs引起的相位误差。还有一 类是自聚焦的方法。这种方法仅仅针对单一SAR影像开展而不需要任何干涉过程,通过探测 高信噪比估计RMEs。这两类方法大大提高了SAR影像的配准精度和干涉相干性，但这两种 方法需要一定的前提，比如前者要求高相干性和后者要求高信噪比目标。然而在植被区，上 述要求难以得到较好的满足，但只要每一距离行至少存在两个独立的高相干像素，通过升级 的MS方法基本可以估计可靠的时频RMEs。

(2)低频RMEs，是未知积分项，主要包含未知基线偏差和线性基线误差等，在干涉图 中体现为相位屏。低频RMEs是上述两方法自身难以解决的，目前估计低频RMEs(导致的相位趋势)的方法可以分成两类，一类是借助外部数据，如控制点或DEM数据。基于控制 点数据的方法被广泛用于拟合平坦地区的基线偏差，而不适合地形起伏区，即便是增加控制点的数量。而且低频RMEs通常沿方位向是一种非线性程度不高的变化趋势，只要存在未知的地形起伏区域，那么RMEs与地形起伏的耦合误差就仍然不能被改正。基于外部DEM的 方法是目前纠正低频RMEs效果较好的一种策略，例如基于地形和孔径相关的补偿方法、基 于小波分解和多项式拟合的方法以及单基线参数化方法等。而这些方法都是针对单基线InSAR数据，最初衷的目的是考虑借助高精度DEM去估计RMEs和相应的相位误差。而在 地表覆盖层(如森林、冰雪和沙漠等)介质中，要找到与机载高分辨率SAR数据相匹配的高 精度DEM是很难的，特别是长波SAR视野下。另一类是多基线相位定标法，这种方法是在 复数域开展RMEs引起的相位误差定标，或是通过探测稳定的点目标或者是通过双定位法迭 代估计目标端和SAR传感器端的误差。然而，当在森林、冰雪区等地表覆盖层介质场景或者 存在残余的高频RMEs时，稳定的点目标探测会受到限制。另外，在复数域开展RMEs估计 不免会存在RMEs的半波长模糊度问题，造成RMEs引起的相位误差不能被完全修正。而且 当前，它们主要服务于森林或冰川场景的垂直结构反演，实际上不希望观测相位中有明显起 伏的地形信号存在。然而，事实上场景中不免存在地形起伏情形，此时一方面RMEs引起的 相位误差特征会被地形相位掩盖，导致RME很难被准确估计，另一方面，由RMEs与地形 相互耦合会进一步扭曲真实的高程信息。

因此，在地形测绘和SAR层析应用中，很有必要设计一种考虑地形起伏而且不依赖外部 高精度DEM数据的RMEs精确估计方法。

发明内容

本发明要解决的技术问题在于，提供一种顾及地形起伏的InSAR残余运动误差估计方法 及装置，有助于提高机载SAR数据用于地形起伏区的地形测绘和森林垂直结构反演的精度。

为实现上述技术目的，本发明采用如下技术方案：

一种顾及地形起伏的InSAR残余运动误差估计方法，包括以下步骤：

估计从影像对应的地形相位：获取研究区域的N+1景SAR影像，依次进行配准、重采样、去平地相位和干涉处理生成N*(N+1)/2个干涉对，假设公共主影像的相位为0，根据干涉得到的N*(N+1)/2个干涉对估计N个从轨迹中包含InSAR残余运动误差的地形相位；

获取差分相位观测值：利用覆盖研究区域的外部DEM数据和InSAR干涉几何参数，计 算模拟地形相位；将包含InSAR残余运动误差的地形相位与模拟地形相位求差，并对差值进 行相位解缠处理，再将解缠差分干涉相位作为观测值；

构建顾及地形起伏的差分干涉相位参数化模型：根据包含InSAR残余运动误差的地形相 位和模拟地形相位分别关于InSAR干涉几何的计算式，计算得到差分干涉相位关于InSAR干 涉几何、DEM误差以及基线误差分量的表达式，并将该表达式作为差分干涉相位参数化模型；

构建多基线观测方程组：针对每个从影像的每一距离行像素，均将差分干涉相位的观测 值与差分干涉相位参数化模型进行匹配构建等式，得到观测方程；所有观测方程构成多基线 观测方程组；

解算残余运动误差：基于已知的InSAR干涉几何参数和观测值，解算多基线观测方程组 中的DEM误差和基线误差分量；其中的基线误差分量即为待估计的残余运动误差。

本技术方案考虑到地形起伏对应的干涉地形相位将会掩盖残余运动误差引起的相位误 差，不利于残余运动误差估计，故引入外部DEM数据进行去地形处理，消除或削弱地形相 位对残余运动误差相位估计的影响；同时通过多基线观测约束DEM误差，确保地形真实高 程不受影响，从而能更准确地估计残余运动误差。

在更优的技术方案中，上述方法还包括按预设迭代次数，迭代计算新的残余运动误差和 DEM误差，将最后一轮迭代计算得到的残余运动误差作为最终的残余动作误差；其中，每一 轮迭代计算方法为：

(1)使用上一轮解算得到的DEM误差，对当前的外部DEM数据的地形高程进行补偿，得到更新的外部DEM数据；

(2)使用更新的外部DEM数据和InSAR干涉几何参数，重新构建观测方程和解算残余 运动误差，得到新一轮的残余运动误差和DEM误差。

本技术方案不依赖外部高精度的DEM数据，而是采用当前全球免费发布的外部DEM数 据，通过多基线约束从而重建更高精度的DEM数据，进而解算得到InSAR残余运动误差，因此本技术方案对外部DEM数据的精度要求较低，同样能准确估计得到InSAR残余运动误差。

在更优的技术方案中，模拟地形相位的计算式为：

式中，

表示第i个模拟地形相位，h代表外部DEM数据中的地形高程，B_⊥,i代表正交于第i个从轨迹的斜距的基线分量，λ表示波长，R_i0代表第i个从轨迹成像几何下的斜距， θ₀为公共主影像的视角；

包含InSAR残余运动误差的地形相位关于InSAR干涉几何参数的计算式为：

式中，

表示第i个从轨迹的干涉地形相位，ΔB_y,i和ΔB_z,i代表第i个从轨迹的基线误差 水平分量和垂直分量，θ_i代表第i个从轨迹的入射角，

代表真实地形高程；

计算得到的差分干涉相位参数化模型为：

式中，

表示第i个从轨迹对应的差分干涉相位，Δh为DEM误差。

在更优的技术方案中，设每个从影像在距离向有n个像素点，故可由N个从影像构建得 到以下观测方程组：

式中，

表示与第i个从轨迹、第j个距离行像素对应的差分干涉相位，Δh^j表示第j 距离行对应的DEM误差；

采用加权截断奇异值分解方法解算残余运动误差和DEM误差，具体方法为：

(1)将N*n个观测方程表示为矩阵形式为：

L＝A·X+δ

其中：

L是N*n的观测列向量，

A是(N*n)×(n+2*N)的系数矩阵，即A＝[sub_A_(N*n)×2*N sub_B_(N*n)×n]， sub_A＝diag[sA1_n×2,...,sAi_n×2,...,sAN_n×2]，sub_B＝[sB1_n×n,...,sBi_n×n,...,sBN_n×n]，

X是n+2*N的未知数列向量，X＝[ΔB_z,1 ΔB_y,1 … ΔB_z,N ΔB_y,N Δh¹ … Δhⁿ]^T；

δ为观测值中未建模的不确定向量；

(2)将未知数列向量X表示为：

X＝A₁ ⁺L₁；

其中，A₁ ⁺表示A₁的广义逆矩阵，A₁、L₁的具体形式为：A₁＝D^-1·A,L₁＝D^-1·L；D表示对权矩阵P的逆阵进行Cholesky分解得到的三角阵，权矩阵P是所有SP相位的相干值的权重系数矩阵；

然后对矩阵A₁进行奇异值分解得：

A₁＝USV^T；

式中U是N*n)×(N*n)(的酉矩阵，V为(n+2*N)×(n+2*N)的酉矩阵，S是 (N*n)×(n+2*N)的非负对角阵，S＝Diag(s₁,s₂,...,s_n+2*N)的对角元素为相应的奇异值；

其中，加权截断奇异值分解的阈值thres为：

thres＝Th·max{s_j} j＝1,2,...,n+2*N；

式中，Th为阈值比例系数，max{s_j}为求取最大值操作；在条件s_r＞thres的限制下，矩 阵S的广义逆为：

从而可以解算得到未知数列向量：

X＝VS^-1U^T·D^-1L。

在更优的技术方案中，计算N个包含InSAR残余运动误差的地形相位的具体方法为：

借助地物目标的统计特性和假设公共主影像的相位为0的先验信息，根据干涉生成的 (N+1)*N/2个干涉对，按以下公式采用极大似然估计法估计N个对应从轨迹的缠绕相位

式中，

符号“T”代表转置；Γ代表所有干涉 对中对应某一像素的相干性构成的(N+1)x(N+1)矩阵；I_Ω是所有干涉对中对应某一像素的复相 干系数构成的(N+1)x(N+1)矩阵。

在更优的技术方案中，所述N个从轨迹对应的地形相位中的地表形变均为0。

在更优的技术方案中，N≥3。

本发明还提供一种顾及地形起伏的InSAR残余运动误差估计装置，包括：

从影像对应的地形相位估计模块，用于：获取研究区域的N+1景SAR影像，依次进行配准、重采样、去平地相位和干涉处理生成N*(N+1)/2个干涉对，假设公共主影像的相位为0，根据干涉得到的N*(N+1)/2个干涉对估计N个从轨迹中包含InSAR残余运动误差的地形相位；

差分相位观测值获取模块，用于：利用覆盖研究区域的外部DEM数据，计算模拟地形 相位；将包含InSAR残余运动误差的地形相位与模拟地形相位求差，并对差值进行相位解缠 处理，再将解缠差分干涉相位作为观测值；

顾及地形起伏的差分干涉相位参数化模型构建模块，用于：根据包含InSAR残余运动误 差的地形相位和模拟地形相位分别关于InSAR干涉几何的计算式，计算得到差分干涉相位关 于InSAR干涉几何、DEM误差以及基线误差分量的表达式，并将该表达式作为差分干涉相 位参数化模型；

多基线观测方程组构建模块，用于：针对每个从影像的每一距离行像素，均将差分干涉 相位的观测值与差分干涉相位参数化模型进行匹配构建等式，得到观测方程；所有观测方程 构成多基线观测方程组；

残余运动误差解算模块，用于：基于已知的InSAR干涉几何参数和观测值，解算多基线 观测方程组中的DEM误差和基线误差分量；其中的基线误差分量即为待估计的残余运动误 差。

在更优的技术方案中，上述装置还包括数据更新模块，用于：按预设迭代次数，迭代计 算新的残余运动误差和DEM误差，将最后一轮迭代计算得到的残余运动误差作为最终的残 余动作误差；其中，每一轮迭代计算方法为：

使用上一轮解算得到的DEM误差，对当前的外部DEM数据的地形高程进行补偿，得到 更新的外部DEM数据；

使用更新的外部DEM数据和InSAR干涉几何参数，重新构建观测方程和解算残余运动 误差，得到新一轮的残余运动误差和DEM误差。

有益效果

本发明的有益效果包括：

(1)本发明考虑地形起伏对应的干涉地形相位将会掩盖残余运动误差引起的相位误差， 故引入外部DEM数据进行去地形处理，消除或削弱地形相位对残余运动误差相位的影响。 这样，不仅大幅削弱了地形相位对RMEs提取的影响，保证了差分干涉相位参数化模型的合 理构建。

(2)本发明不依赖外部高精度的DEM数据，而是采用当前全球免费发布的外部DEM数据，从而通过多基线重建更高精度的DEM数据，进而解算得到InSAR残余运动误差，因 此本技术方案对外部DEM数据的精度要求较低，同样能准确估计得到InSAR残余运动误差。

(3)在计算观测值时，采用的Phase Linking算法，一方面，大幅削弱多基线中难以避 免的失相干问题的影响，以确保观测值中相位跳变的降低和相位解缠的成功性；另一方面， 最大程度上降低了多基线观测方程组的维度和未知参数解算的计算量，提高了计算效率。

(4)在参数解算方面，一方面，多基线观测方程组是按距离行构建，有助于非线性残余 运动误差的提取；另一方面，采用加权截断奇异值分解(weighted TSVD)参数估计方法，避 免了参数解算过程中出现很大的偏差或者解算不稳定的情形。进而保证了模型参数中DEM误 差的准确估计，确保了高精度DEM重建的精度，为提取高精度的RMEs奠定了基础。

(5)本发明原理直观，易于编程实现及普适性应用，可作为一种稳健可靠的通用算法， 以广泛用于高精度、高分辨率数字高程模型(DEM)和地表覆盖层(如森林、冰川)垂直结构的重建和有效提高测量精度。

附图说明

图1为InSAR获取几何下RMEs示意图；

图2为本发明实施例一种顾及地形起伏的InSAR残余运动误差估计方法的流程图；

图3为本发明实施例中研究概况和步骤III中用到TanDEM-X DEM数据；

图4，其中(a)为本实施例估计的DEM误差统计图估计的DEM误差统计图；(b)本实施案 例重建的高精度DEM；(c)为高精度林下LiDAR DEM数据参考图；

图5为本实施例中步骤I得到的解缠的差分干涉相位图和步骤V中重建的高精度DEM 前后估计的RMEs相位误差图：(a-e)解缠的差分干涉相位图；(f-j)重建高精度DEM前估计 的RMEs相位误差图；(k-o)重建高精度DEM后估计的RMEs相位误差图；(p-t)重建高 精度DEM前后两次估计的RMEs相位误差的差值图；

图6为本实施例中步骤V重建高精度DEM前估计得到的基线误差估值。

具体实施方式

为了使本发明的目的、技术方案及其优点更加清楚明白，采用德国宇航局通过机载E-SAR 系统于2008年飞行试验期间获取的P波段多基线数据，具体参数见附表1，对本发明进行深 入详细地说明。SAR影像数据覆盖区的地形起伏为海拔150-380m，且该地区为典型的植被覆 盖区，平均树高为18m，最高的树高为30m。由于P波段具有较强的穿透性，因此DEM误 差中很大一部分是由植被相位中心差异贡献。此处所描述的具体实施案例仅用以解释本发明， 并不用于限定本发明。

基于上述所采用的机载E-SAR系统P波段多基线干涉数据，如图2所示，本发明顾及地 形起伏的InSAR残余运动误差估计方法的具体实施步骤如下：

步骤I，获取差分相位观测值。

步骤1.1，先对机载P波段6景SAR影像进行基于公共主影像(0103)的配准和重采样， 并基于E-SAR系统P波段SAR数据的干涉几何参数，按照公式(1)进行平地相位估计，得到5个从影像的平地相位

i＝1,2,…5；然后按公式(2)从从影像相位ω_i'中去除平地相位

得到从影像去平地后的相位ω_i；再对此时的从影像进行干涉处理，生成15个互不重合的干涉对。

式中，R₀分别表示主影像的SAR传感器天线中心到参考面的成像几何斜距；H代表SAR 平台的标称飞行高度；B_y,i,B_z,i分别代表第i个干涉几何的标称水平基线和垂直基线分量。

步骤1.2，依据地物目标的统计特性和假设公共主影像(0103)的相位为0的先验信息， 根据上一步得到的15个干涉对，然后采用极大似然估计法(即按公式(3))估计5个从影像的缠绕相位，即包含有InSAR残余运动误差的地形相位(Slave-images phase，简写为SP)

式中，

符号“T”代表转置；Γ代表所有干涉 对中对应某一像素的相干性构成的6×6矩阵；I_Ω是所有干涉对中对应某一像素的复相干系数 构成的6×6矩阵。

步骤1.3，对步骤1.2处理后得到的影像进行多视处理，具体按方位向与距离向比为2：1 进行多视处理，使其空间分辨率为1.5米*1.5米。

步骤1.4，基于表1中的InSAR干涉几何参数和外部DEM数据，按以下公式(4)计算研究区域的模拟地形相位

式中，

表示第i个模拟地形相位，h代表外部DEM数据中的地形高程，B_⊥,i代表 正交于第i个从轨迹的斜距的基线分量，λ表示波长，R_i0代表第i个从轨迹成像几何下的斜距， θ₀为公共主影像的视角；

在本实施例中，外部DEM数据采用的是德国宇航局免费发布的覆盖该地区的TanDEM-X DEM数据(分辨率90米*90米)，见附图3，将其过采样到1.5x1.5米使其与步骤1.3得到 的地形相位图的分辨率和行列值均一致，并编码到SAR坐标系，即可按像素点从该SAR坐 标系中获取相应的地形高程h。

步骤1.5，按以下公式(5)对步骤1.2得到的地形相位

进行去平地和去地形操作，得 到未解缠的差分干涉相位

步骤1.6，采用Goldstein滤波算法分别对5个未解缠的差分干涉相位

进行滤波，滤 波窗口设置为32×32像素，滤波因子为0.35，从而得到滤波的差分干涉相位；然后利用最小 网络费用流算法分别对滤波后的5个差分干涉相位进行相位解缠操作，解缠所采用的加权因 子为滤波前或滤波后的相干图。为了避免低相干区域对解缠的影响，本实施例在解缠的过程 中，掩膜掉相干性低于0.3的像元点，最后获得解缠的差分干涉相位

并将该解缠后的 差分干涉相位

作为观测值，差分干涉相位图可见附图5(a-e)所示。

步骤II，构建顾及地形起伏的差分干涉相位参数化模型。

步骤2.1，构建包含InSAR残余运动误差的地形相位关于InSAR干涉几何的计算式为：

式中，

表示第i个从轨迹的干涉地形相位，ΔB_y,i和ΔB_z,i代表第i个从轨迹的基线误差 水平分量和垂直分量，θ_i代表第i个从轨迹的入射角，

代表真实地形高程。

第i个缠绕的地形相位

采用InSAR干涉几何参数的表达式如以下公式(7)所示：

式中，λ表示波长，R_i分别表示第i个从影像的SAR传感器天线中心到地面目标的斜距。

由于主从轨迹中残余运动误差(简写为RMEs)，即ΔR₀,ΔR_i的存在，公式(7)可转化为公式(8)所示：

式中，

分别表示主影像和从影像的SAR传感器天线中心到地面目标的真实斜距， 而其中，参考图1所示可得：

式中，

代表SAR平台的飞行高度；

分别代表第i个SP干涉几何的真实水平基 线和真实垂直基线分量；

代表InSAR获取几何的真实地形高程。

众所周知，当前机载SAR系统搭载的惯导系统精度有限(1-10cm)，故上述公式(9)中真实的几何参数(如

)是不可能已知的，而是含有误差的标称几何参数。因而上述几何参数的误差都有可能对RMEs有一定贡献，但由于缺少基准，这些参数的误差估计实际上是一个秩亏的问题。为了解决此问题，需要假设公共主影像没有RMEs的影响，因此主轨迹被固定，即

ΔR₀＝0，认为RMEs全部来源于从轨迹。因此公式(8)中ΔR_i的产生只与标称基线参数(B_y,i,B_z,i)有关，由公式(9)可得公式(10)为：

式中，ΔB_y,i和ΔB_z,i分别代表第i个SP相位中的水平基线误差分量和垂直基线误差分量； θ_i代表第i个从轨迹中的入射角。

由公式(10)很明显可得，残余运动误差RMEs的待求问题实际上被转移到了基线误差 分量的求解上。当地形很平坦或接近参考面时，RMEs近似只包含平地相位误差项；但从公 式(9)中的中间部分可以看出，当存在地形起伏时，RMEs除了包含残余平地项以外，还包含RMEs与地形存在耦合关系项，即公式(10)就是地形起伏情形下干涉相位中RMEs成分 的参数化表达。

根据相位成分的分析，通常用于地形测绘和植被垂直结构反演的机载SAR飞行模式设置 的时间基线比较短(通常在几个小时以内)，所以可不考虑地表形变信号和时间失相干影响； 而且机载SAR平台高度比较低、多数机载SAR系统飞行试验都选在晴好天气环境开展，因 此大气延迟影响也可以忽略。因此，可将

由公式(8)进一步简化得到公式(6)，即包 含InSAR残余运动误差的地形相位的计算式。

步骤2.2，为了能探测RMEs的影响和为了能准确估计RMEs项做准备，需要消除地形相 位，因此本实施例引入当前全球免费发布的TanDEM-X DEM数据(还可以为Aster DEM、SRTM DEM等外部DEM数据)，并基于InSAR干涉几何参数来构建差分干涉相位。因此在 该步骤2.2中，根据公式(4)所示的模拟地形相位关于InSAR干涉几何的计算式，以及步骤 2.1得到的包含InSAR残余运动误差的地形相位关于InSAR干涉几何的计算式，计算得到差 分干涉相位关于InSAR干涉几何、DEM误差以及基线误差分量的表达式，如公式(11)所 示，并将公式(11)作为差分干涉相位参数化模型；

式中，

表示第i个从轨迹对应的差分干涉相位，Δh为DEM误差，ΔB_y,i为基线误差 水平分量，ΔB_z,i为基线误差垂直分量。

步骤III，构建多基线观测方程组：针对每个从影像的每一距离行像素，均将差分干涉相 位的观测值与差分干涉相位参数化模型进行匹配构建等式，得到观测方程；假设同一像素点 在每个差分干涉相位中的DEM误差相同，所有观测方程构成如公式(12)所示的多基线观 测方程组：

式中，

表示与第i个从轨迹、第j个距离行像素对应的差分干涉相位，Δh^j表示第j 个距离行像素对应的DEM误差，n为设每个从影像在距离向的像素点个数，N为从影像的数 量。在实施例中，N＝5，n＝1501。

步骤IV，解算残余运动误差，基于已知的InSAR干涉几何参数(包括λ、

R_i0、

)和观测值

解算多基线观测方程组中的DEM误差Δh^j和基线误差分量；其中的基 线误差分量即为待估计的残余运动误差。

观测方程组中，观测方程的个数为N*n，未知数包括n个DEM误差、N个水平基线分量误差和N个垂直基线分量误差，只要满足N*n≥n+2N，观测方程组即是超定可解。

步骤4.1，将N*n个观测方程表示为矩阵形式为：

L＝A·X+δ (13)

其中：

L是N*n的观测列向量，

A是(N*n)×(n+2*N)的系数矩阵，即A＝[sub_A_(N*n)×2*N sub_B_(N*n)×n]， sub_A＝diag[sA1_n×2,...,sAi_n×2,...,sAN_n×2]，sub_B＝[sB1_n×n,...,sBi_n×n,...,sBN_n×n]，

X是n+2*N的未知数列向量，X＝[ΔB_z,1 ΔB_y,1 … ΔB_z,N ΔB_y,N Δh¹ … Δhⁿ]^T；

δ为观测值中未建模的不确定向量；

步骤4.2，采用加权截断奇异值分解的方法解算观测方程组中的未知数列向量X；

由于本实施例中的差分相位参数化模型是基于多基线参数化建立得到，因此系数矩阵存 在病态问题，因此可采用截断奇异值分解(TSVD)方法解算方程组。但为了进一步确保模型 解算的稳健性，本实施例采用加权截断奇异值分解策略解算未知数列向量X，表示为：

X＝A₁ ⁺L₁ (14)；

其中，A₁ ⁺表示A₁的广义逆矩阵，A₁、L₁的具体形式为：A₁＝D^-1·A,L₁＝D^-1·L；D表示对权矩阵P的逆阵进行Cholesky分解得到的三角阵，权矩阵P是所有SP相位的相干值的权重系数矩阵；

然后对矩阵A₁进行奇异值分解得：

A₁＝USV^T (15)；

式中U是N*n)×(N*n)(的酉矩阵，V为(n+2*N)×(n+2*N)的酉矩阵，S是 (N*n)×(n+2*N)的非负对角阵，S＝Diag(s₁,s₂,...,s_n+2*N)的对角元素为相应的奇异值；

其中，加权截断奇异值分解的阈值thres为：

thres＝Th·max{s_j} j＝1,2,...,n+2*N (16)；

式中，Th为阈值比例系数，max{s_j}为求取最大值操作；在条件s_r＞thres的限制下，矩 阵S的广义逆为：

从而可以解算得到未知数列向量：

X＝VS^-1U^T·D^-1L (18)。

步骤V，迭代计算并输出最终的残余运动误差：按预设迭代次数，迭代计算新的残余运 动误差和DEM误差，将最后一轮迭代计算得到的残余运动误差作为最终的残余动作误差； 其中，每一轮迭代计算方法为：

步骤5.1，使用上一轮解算得到的DEM误差，对当前的外部DEM数据的地形高程进行补偿，得到更新的外部DEM数据；

在本实施例中，每轮迭代由步骤IV解算得到未知数列向量后，先分析其中的DEM误差， 其中第一轮迭代统计得到的DEM误差见附图4(a)所示；再利用该第一轮迭代得到DEM误差，执行步骤5.1即可重建得到精度更高的DEM数据，可见图4(b)所示。很明显，植被 相位中心高度差异主导了DEM误差。另外，估计的RMEs及其相应的相位误差见附图6的 虚线和图5(f-j)。估计的相位误差结果与原始观测值中的相位趋势很接近，但由于P波段 具有强穿透性，与TanDEM-X DEM的高程相差较大，因此需要进行重建高精度DEM。可以 看到重建后的DEM数据中细节信息的表达更丰富。

步骤5.2，使用更新的外部DEM数据和InSAR干涉几何参数，按步骤III和步骤IV重新 构建观测方程和解算残余运动误差，得到新一轮的残余运动误差和DEM误差。

在本实施例第二轮迭代计算后，得到了P波段5个相位中存在的RMEs相位误差，见附 图6的实线和图5(k-o)。由此可见，两次估计结果存在明显的偏差，而重新估计的RMEs的可靠性实际上取决于步骤5.1重建得到的更高精度的DEM数据。为此，采用高精度林下LiDAR DEM数据作为参考，精度为2.8m，见附图4(c)。这证明了本实施例重建DEM数 据已经足够准确，也证明了本发明实施例的RMEs精确估计方法的可靠性。

由于DEM误差或者观测值中的不确定性，重建后的DEM数据难以再得到较大提高，因 此重新估计后的RMEs也不会有明显改进，故这里通常迭代一次或两次即可。

本发明还提供一种顾及地形起伏的InSAR残余运动误差估计装置，与上述方法实施例方 案相对应，包括：

从影像对应的地形相位估计模块，用于：获取研究区域的N+1景SAR影像，依次进行配准、重采样、去平地相位和干涉处理生成N*(N+1)/2个干涉对，假设公共主影像的相位为0，根据干涉得到的N*(N+1)/2个干涉对估计N个从轨迹中包含InSAR残余运动误差的地形相位；

差分相位观测值获取模块，用于：利用覆盖研究区域的外部DEM数据，计算模拟地形 相位；将包含InSAR残余运动误差的地形相位与模拟地形相位求差，并对差值进行相位解缠 处理，再将解缠差分干涉相位作为观测值；

顾及地形起伏的差分干涉相位参数化模型构建模块，用于：根据包含InSAR残余运动误 差的地形相位和模拟地形相位分别关于InSAR干涉几何的计算式，计算得到差分干涉相位关 于InSAR干涉几何、DEM误差以及基线误差分量的表达式，并将该表达式作为差分干涉相 位参数化模型；

多基线观测方程组构建模块，用于：针对每个从影像的每一距离行像素，均将差分干涉 相位的观测值与差分干涉相位参数化模型进行匹配构建等式，得到观测方程；所有观测方程 构成多基线观测方程组；

残余运动误差解算模块，用于：基于已知的InSAR干涉几何参数和观测值，解算多基线 观测方程组中的DEM误差和基线误差分量；其中的基线误差分量即为待估计的残余运动误 差。

在更优的实施例中，上述装置还包括数据更新模块，用于：按预设迭代次数，迭代计算 新的残余运动误差和DEM误差，将最后一轮迭代计算得到的残余运动误差作为最终的残余 动作误差；其中，每一轮迭代计算方法为：

使用上一轮解算得到的DEM误差，对当前的外部DEM数据的地形高程进行补偿，得到 更新的外部DEM数据；

使用更新的外部DEM数据和InSAR干涉几何参数，重新构建观测方程和解算残余运动 误差，得到新一轮的残余运动误差和DEM误差。

表1机载E-SAR系统P波段多基线数据的干涉几何参数

以上实施例为本申请的优选实施例，本领域的普通技术人员还可以在此基础上进行各种 变换或改进，在不脱离本申请总的构思的前提下，这些变换或改进都应当属于本申请要求保 护的范围之内。

Claims (9)

1.一种顾及地形起伏的InSAR残余运动误差估计方法，其特征在于，包括以下步骤：

估计从影像对应的地形相位：获取研究区域的N+1景SAR影像，依次进行配准、重采样、去平地相位和干涉处理生成N*(N+1)/2个干涉对，假设公共主影像的相位为0，根据干涉得到的N*(N+1)/2个干涉对估计N个从轨迹中包含InSAR残余运动误差的地形相位；

获取差分相位观测值：利用覆盖研究区域的外部DEM数据和InSAR干涉几何参数，计算模拟地形相位；将包含InSAR残余运动误差的地形相位与模拟地形相位求差，并对差值进行相位解缠处理，再将解缠差分干涉相位作为观测值；

构建顾及地形起伏的差分干涉相位参数化模型：根据包含InSAR残余运动误差的地形相位和模拟地形相位分别关于InSAR干涉几何的计算式，计算得到差分干涉相位关于InSAR干涉几何、DEM误差以及基线误差分量的表达式，并将该表达式作为差分干涉相位参数化模型；

构建多基线观测方程组：针对每个从影像的每一距离行像素，均将差分干涉相位的观测值与差分干涉相位参数化模型进行匹配构建等式，得到观测方程；所有观测方程构成多基线观测方程组；

解算残余运动误差：基于已知的InSAR干涉几何参数和观测值，解算多基线观测方程组中的DEM误差和基线误差分量；其中的基线误差分量即为待估计的残余运动误差。

2.根据权利要求1所述的方法，其特征在于，还包括按预设迭代次数，迭代计算新的残余运动误差和DEM误差，将最后一轮迭代计算得到的残余运动误差作为最终的残余动作误差；其中，每一轮迭代计算方法为：

(1)使用上一轮解算得到的DEM误差，对当前的外部DEM数据的地形高程进行补偿，得到更新的外部DEM数据；

(2)使用更新的外部DEM数据和InSAR干涉几何参数，重新构建观测方程和解算残余运动误差，得到新一轮的残余运动误差和DEM误差。

3.根据权利要求1所述的方法，其特征在于，模拟地形相位的计算式为：

式中，

表示第i个模拟地形相位，h代表外部DEM数据中的地形高程，B_⊥,i代表正交于第i个从轨迹的斜距的基线分量，λ表示波长，R_i0代表第i个从轨迹成像几何下的斜距，θ₀为公共主影像的视角；

包含InSAR残余运动误差的地形相位关于InSAR干涉几何参数的计算式为：

式中，

表示第i个从轨迹的干涉地形相位，ΔB_y,i和ΔB_z,i代表第i个从轨迹的基线误差水平分量和垂直分量，θ_i代表第i个从轨迹的入射角，

代表真实地形高程；

计算得到的差分干涉相位参数化模型为：

式中，

表示第i个从轨迹对应的差分干涉相位，Δh为DEM误差。

4.根据权利要求3所述的方法，其特征在于，设每个从影像在距离向有n个像素点，故可由N个从影像构建得到以下观测方程组：

式中，

表示与第i个从轨迹、第j个距离行像素对应的差分干涉相位，Δh^j表示第j距离行对应的DEM误差；

采用加权截断奇异值分解方法解算残余运动误差和DEM误差，具体方法为：

(1)将N*n个观测方程表示为矩阵形式为：

L＝A·X+δ

其中：

L是N*n的观测列向量，

A是(N*n)×(n+2*N)的系数矩阵，即A＝[sub_A_(N*n)×2*N sub_B_(N*n)×n]，

sub_A＝diag[sA1_n×2,...,sAi_n×2,...,sAN_n×2]，sub_B＝[sB1_n×n,...,sBi_n×n,...,sBN_n×n]，

X是n+2*N的未知数列向量，X＝[ΔB_z,1 ΔB_y,1…ΔB_z,N ΔB_y,N Δh¹…Δhⁿ]^T；

δ为观测值中未建模的不确定向量；

(2)将未知数列向量X表示为：

X＝A₁ ⁺L₁；

其中，A₁ ⁺表示A₁的广义逆矩阵，A₁、L₁的具体形式为：A₁＝D^-1·A,L₁＝D^-1·L；D表示对权矩阵P的逆阵进行Cholesky分解得到的三角阵，权矩阵P是所有SP相位的相干值的权重系数矩阵；

然后对矩阵A₁进行奇异值分解得：

A₁＝USV^T；

式中U是N*n)×(N*n)(的酉矩阵，V为(n+2*N)×(n+2*N)的酉矩阵，S是(N*n)×(n+2*N)的非负对角阵，S＝Diag(s₁,s₂,...,s_n+2*N)的对角元素为相应的奇异值；

其中，加权截断奇异值分解的阈值thres为：

thres＝Th·max{s_j}j＝1,2,...,n+2*N；

式中，Th为阈值比例系数，max{s_j}为求取最大值操作；在条件s_r＞thres的限制下，矩阵S的广义逆为：

从而可以解算得到未知数列向量：

X＝VS^-1U^T·D^-1L。

5.根据权利要求1所述的方法，其特征在于，计算N个包含InSAR残余运动误差的地形相位的具体方法为：

借助地物目标的统计特性和假设公共主影像的相位为0的先验信息，根据干涉生成的(N+1)*N/2个干涉对，按以下公式采用极大似然估计法估计N个对应从轨迹的缠绕相位

式中，

符号“T”代表转置；Γ代表所有干涉对中对应某一像素的相干性构成的(N+1)x(N+1)矩阵；I_Ω是所有干涉对中对应某一像素的复相干系数构成的(N+1)x(N+1)矩阵。

6.根据权利要求1所述的方法，其特征在于，所述N个从轨迹对应的地形相位中的地表形变均为0。

7.根据权利要求1所述的方法，其特征在于，N≥3。

8.一种顾及地形起伏的InSAR残余运动误差估计装置，其特征在于，包括：

从影像对应的地形相位估计模块，用于：获取研究区域的N+1景SAR影像，依次进行配准、重采样、去平地相位和干涉处理生成N*(N+1)/2个干涉对，假设公共主影像的相位为0，根据干涉得到的N*(N+1)/2个干涉对估计N个从轨迹中包含InSAR残余运动误差的地形相位；

差分相位观测值获取模块，用于：利用覆盖研究区域的外部DEM数据，计算模拟地形相位；将包含InSAR残余运动误差的地形相位与模拟地形相位求差，并对差值进行相位解缠处理，再将解缠差分干涉相位作为观测值；

顾及地形起伏的差分干涉相位参数化模型构建模块，用于：根据包含InSAR残余运动误差的地形相位和模拟地形相位分别关于InSAR干涉几何的计算式，计算得到差分干涉相位关于InSAR干涉几何、DEM误差以及基线误差分量的表达式，并将该表达式作为差分干涉相位参数化模型；

多基线观测方程组构建模块，用于：针对每个从影像的每一距离行像素，均将差分干涉相位的观测值与差分干涉相位参数化模型进行匹配构建等式，得到观测方程；所有观测方程构成多基线观测方程组；

残余运动误差解算模块，用于：基于已知的InSAR干涉几何参数和观测值，解算多基线观测方程组中的DEM误差和基线误差分量；其中的基线误差分量即为待估计的残余运动误差。

9.根据权利要求8所述的装置，其特征在于，还包括数据更新模块，用于：按预设迭代次数，迭代计算新的残余运动误差和DEM误差，将最后一轮迭代计算得到的残余运动误差作为最终的残余动作误差；其中，每一轮迭代计算方法为：

使用上一轮解算得到的DEM误差，对当前的外部DEM数据的地形高程进行补偿，得到更新的外部DEM数据；

使用更新的外部DEM数据和InSAR干涉几何参数，重新构建观测方程和解算残余运动误差，得到新一轮的残余运动误差和DEM误差。

Images