一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法

摘要

一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法，它有五大步骤：步骤一：高分辨率SAR数据选取；步骤二：初始高程相位估计；步骤三：选取InSAR时序分析干涉像对数据集并提取相干目标候选点；步骤四：高程相位校正和线性形变分量估计；步骤五：相干目标的形变序列估计。本发明能够克服城区高分辨率InSAR形变监测数据处理中无高分辨率DEM去除高层建筑物附加相位，有效地将城市高层建筑物在高分辨率InSAR中的附加相位与形变相位进行分离，从而大大提高高分辨率InSAR在城区地面沉降中的监测精度。

说明书

技术领域

本发明涉及一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法，属于合成孔径雷达干涉测量技术(InSAR)领域。它能够有效地将城市高层建筑物在高分辨率InSAR干涉图中的附加相位与形变相位进行分离，可以大大提高高分辨率InSAR在城区地面沉降中的监测精度。

背景技术

利用InSAR求解高程(DEM)或者形变量，相位解缠是必须的步骤之一，这一步骤在差分干涉处理中称为解缠，而在时序分析方法中则称为参数估计，其实质是求解缠绕相位的整周数，基本过程是求解相邻点间的相位差，然后按照特定的路径或网络以一定的约束条件进行积分，求解观测范围整体的解缠相位，进而反演形变场或高程。相位解缠的前提是干涉图连续分布且变化平缓，满足相位差小于π的约束条件，整个相位场为无旋场，解缠结果不随路径而变。而实际上干涉图受噪声或不连续相位(如不连续的陡坡)的影响，难以满足解缠条件，因而需要按照给定路径求解后再进行连接。高分辨率条件下城市高层建筑物所引起的相位变化类似于非连续的陡坡相位。对于求解高程而言，条纹密度随基线的增大而密集，且存在与形变条纹混合的可能。因而，如何同时解算高程和形变相位，提高形变量估计的精度是高分辨率InSAR应用面临的主要难题。

地形相位补偿是InSAR时序分析处理过程中差分相位求取的基本步骤，同时地物高程也是InSAR时序分析的主要解算结果。中等分辨率SAR数据中表现为一个或者几个点目标的像元在高分辨率雷达数据中以数个点目标表现出来，每个散射体的散射中心不同，如一栋建筑物，在中分SAR图像中表现为一个点目标，而在高分SAR图像中以多个不同的散射体表现出来。对高分辨率InSAR地形相位补偿的最直接方法是获取每个散射中心高精度的DEM数据，利用Lidar、TanDEM-X等获取的DEM数据实现近似同等分辨率差分干涉处理的高程相位补偿。然而这种高分辨率的DEM数据很多情况下难以获得。因此，这种通过外部DEM(如SRTM)补偿整个像元高程的方法在高分辨率数据处理中将难以适用，需要研究单个目标高程精确计算的方法。

在InSAR干涉图中，干涉条纹的频率与基线的关系如图1所示，同等高程下，基线越短，条纹越稀疏，利于相位解缠或参数估计。当地形变化为2π时对应的高程差为：

>

这里，Δh_2π为高程模糊度，表征干涉相位对地形高程变化的灵敏度。相对于高程变化的影响，形变相位与高程无关，但受到高程变化的影响。

根据形变相位和高程相位的关系，可知高程变化对形变相位的影响与基线相关：

$>\mathrm{\delta d}=\frac{B_{\perp }·\mathrm{\delta h}}{R\sin \theta }---\left(2\right)$>

以Cosmo-skymed数据为例，图2表征了高程误差对形变量的影响。当垂直基线为50m时，50m的高程差引起的相位变化小于0.41个整周，100m时为0.82个整周。同样的条件下，高程差为100m时，100m基线所引起的条纹变化为1.64个整周，是基线为50m时的2倍。显然，短基线条件下条纹密度降低，这有利于对形变信息的准确提取。此外，顾及到长基线有利于高程反演精度，本发明考虑综合利用长短基线干涉图迭代修正高程估计值，以达到提高形变估计精度的目的。

发明内容

1.目的：本发明的目的是提供一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法。它能够克服城区高分辨率InSAR形变监测数据处理中无高分辨率DEM去除高层建筑物附加相位，有效地将城市高层建筑物在高分辨率InSAR中的附加相位与形变相位进行分离，从而可以大大提高高分辨率InSAR在城区地面沉降中的监测精度。

2.技术方案：本发明是一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法，该方法具体步骤如下：

步骤一：高分辨率SAR数据选取

以TerraSAR－X和COSMO－Skymed为代表的高分辨率雷达卫星系统为开展InSAR反演城区地物高程和形变精细化监测提供了数据源。高分辨率InSAR相对于中等分辨率InSAR技术，其总体优势体现在两个方面，即(ⅰ)高密度相干点目标和短周期(4-16天)；(ⅱ)对地面点目标的准确定位。高分辨率InSAR数据处理的数据选取要满足空间上尽可能覆盖完整城区的轨道数据，在时间上数据能够连续接收。

步骤二：初始高程相位估计

InSAR形变时序分析所构建的二维参数估计模型中，考虑到大气的空间相关性，对相邻两点(小于大气相关距离)求差以削弱大气的影响。相干目标i和j差分干涉相位的互差为：

$>\Delta {\phi }_{i,j}^k=[C_B·B^{\left(k\right)}·{\mathrm{\Delta ϵ}}_{i,j}+\frac{4\pi }{\lambda }·T^{\left(k\right)}·{\mathrm{\Delta v}}_{i,j}]+{\mu }_{\mathrm{NL}}^{\left(k\right)}+{\alpha }^{\left(k\right)}+n^{\left(k\right)}---\left(3\right)$>

上式中，C_B为与垂直基线相关的系数，T为时间基线，Δε为相对高程误差，Δv为相对形变速率，μ_NL为非线性形变量，α为大气相位，n为噪声，k表示干涉图个数(与干涉图序列的组合有关)。构建目标函数如下：

$>{\phi }_{\mathrm{mod}>}(i,j{,T}^{\left(k\right)})=C_B{·B}^{\left(k\right)}·{\mathrm{\Delta ϵ}}_{i,j}+\frac{4\pi }{\lambda }·T^{\left(k\right)}·{\mathrm{\Delta v}}_{i,j}---\left(4\right)$>

将上式从相位互差式(3)中减去，得到残余相位为：

$>{\mathrm{\Delta w}}_{i,j}^k={\mathrm{\Delta \phi }}_{i,j}^{\left(k\right)}-[C_B·B^{\left(k\right)}·{\mathrm{\Delta ϵ}}_{i,j}+\frac{4\pi }{\lambda }·T^{\left(k\right)}·{\mathrm{\Delta v}}_{i,j}]={\mu }_{\mathrm{NL}}^{\left(k\right)}+{\alpha }^{\left(k\right)}+n^{\left(k\right)---\left(5\right)}$>

显然，当目标函数的参数Δε和Δv在准确估计时，残余相位将最小化。

本发明采取将差分干涉图相位解缠的基础上进行的估计，这时式(3)转换为二维线性函数，可通过建立Delanay三角网或利用冗余网构建更为复杂的连接关系强化对待解算方程组的约束，利用邻近法则将所有距离满足大气相关距离的相干目标连接起来，在求解完成相邻点间的互差后，通过最小二乘或加权平均的方法求解每个目标相对于参考点的高程值和形变速率场。

当干涉图时间间隔和空间基线足够小时，变形相位可忽略，则可以将二维模型(式(3))转为一维模型(式(6))，进而解算得到DEM。

$>\Delta {\phi }_{i,j}^k=C_B·B^{\left(k\right)}·{\mathrm{\Delta ϵ}}_{i,j}+{\alpha }^{\left(k\right)}+n^{\left(k\right)}---\left(6\right)$>

(1)选取超短时间和空间基线干涉像对数据集

依据上述理论，根据建筑物的高度、干涉图基线及地形相位条纹密度的关系，选取研究区平均建筑物高度作为高程约束，确定计算原始地物高程时需要的干涉像对序列。选取时间基线和空间基线都较短的干涉像对序列作为初始求解的数据集，设定初始时间和空间基线阈值分别为T₁和B₁，对数据集中的干涉像对进行干涉处理，无需外部DEM进行地形相位校正，本发明中综合利用快速傅里叶变换(FFT)估计和多项式拟合估计的方法去除干涉图中的轨道误差和趋势性干涉条纹。生成超短时空基线干涉图数据集、相干图数据集和强度图数据集。对每一个干涉图进行相位解缠，按照一维模型多次迭代修正高程误差相位估计得到原始DEM。需要说明的是，在地理编码之前所有的操作都是在雷达坐标系下进行的。

(2)利用点目标识别方法提取相干目标候选点

针对上述所有生成的干涉图序列，本发明中综合采用幅度离散指数(AmplitudeDispersion Index)和相干系数(coherence)来筛选相干目标候选点。

幅度离散指数的计算公式为：

$>D_A=\frac{{\sigma }_A}{m_A}---\left(7\right)$>

其中，σ_A和m_A分别为像素幅度值的标准差和均值。给定一适当阈值D_A低于阈值的像元确定为相干目标候选点。

雷达干涉相位图的相干系数估计公式为：

$>\tilde{\gamma }=\left|\frac{\frac{1}{N}{\Sigma }_{i=0}^N{M}_i{S}_i^{*}}{\frac{1}{N}{\Sigma }_{i=0}^N{M}_i{M}_i^{*}\frac{1}{N}{\Sigma }_{i=0}^N{S}_i{S}_i^{*}}\right|---\left(8\right)$>

根据各像元点在相干图中的相干系数序列γ_i和给定的相干系数阈值如果mean那么则将该像元确定为相干目标候选点。

(3)多次迭代修正高程误差相位估计地物高程

对于每一相干目标点，利用超短时空基线干涉图数据集并按照式(6)所述的一维模型多次迭代修正高程误差相位，将多次估计的高程误差修正相位相加，得到相干目标点处的高程。然后，将点矢量转换为栅格数据，生成研究区的原始DEM，

步骤三：选取InSAR时序分析干涉像对数据集并提取相干目标候选点

将空间基线和时间基线不超过某一设定阈值的所有干涉像对进行干涉处理，将步骤二生成的原始高程模拟地形相位用于所有干涉图的高程相位补偿，生成差分干涉图序列。针对单一差分干涉图中出现的轨道误差和趋势性干涉条纹，综合利用快速傅里叶变换(FFT)估计和多项式拟合估计的方法予以去除，最终生成用于InSAR时序分析的差分干涉图数据集和相应的相干图数据集以及强度图数据集。对上述所有生成的差分干涉图序列，综合采用幅度离散指数和相干系数来筛选相干目标候选点，降低高分辨率SAR条件下相干目标数量冗余。

步骤四：高程相位校正和线性形变分量估计

解缠InSAR时序分析差分干涉图序列中的每一幅差分干涉图。选取时间和空间基线阈值分别为T₂和B₂的干涉像对数据集，重新利用一维模型二次估计高程残余值dh_ori。将上述一维模型所得高程残余值dh_ori作为形变时序分析二维参数估计模型的初始高程，以ΔT、ΔB和Δdh作为时间、空间基线和高程修正的步长，将T₂+ΔT和B₂+ΔB内的干涉图参与到二维参数估计的序列中，设定初始最大高程修正阈值DH_{max_ori}，并初次估计高程修正值和形变速率。分别以步长ΔT和ΔB逐步增加用于迭代估算的干涉图数量，并以步长Δdh逐步减小最大高程修正阈值，通过多次迭代处理逐步修正高程估计值和形变速率。直至满足时间和空间基线以及最大高程修正阈值的所有差分干涉图都参与计算，估计最终的高程改正和形变速率。由于逐次迭代求解得到的高程改正是对已经去除高程值的估计，因而最终的地物高程估计值为逐次迭代修正值之和，即：

$>H=\underset{i=1}{\overset{u}{\Sigma }}{ϵ}^i---\left(9\right)$>

式中H为迭代终止时的地物高程估计值，u为总迭代次数。

步骤五：相干目标的形变序列估计

对于具有显著非线性形变过程，仍需对残余相位进行更为复杂的处理，以提取非线性形变量。处理的前提仍是基于不同相位的时空频率特性。从差分相位中去除高程和线性速率后的残余相位为：

$>{\mathrm{\Delta w}}_{i,j}^k={\mu }_{i,j}^k+a_{i,j}^k+n_{i,j}^k---\left(10\right)$>

由于已解算出整周相位，因而残余相位为相位主值，其大小满足：

$>\left|{\mathrm{\Delta w}}_{i,j}^k\right|<\pi ---\left(11\right)$>

这里，首先要实现对单个差分干涉图中残余相位的滤波处理。由于大气表现出空间低频特性，而非线性变形的空间变化范围较小，相对大气相位而言表现为高通特性。因而，对残余相位图进行空间低通滤波可以进一步弱化大气的影响。在短基线集条件下求解的非线性形变量，不同于PSInSAR的处理策略，不直接对大气进行估计，而是通过滤波减弱大气的影响。残余相位经过空间高通滤波后，大气分量已经减弱。此时，进一步的处理是求解与雷达数据对应的不同时刻的非线性变形量，包含噪声等影响。根据主辅影像的关系，可将相位解缠后的残余相位分解为：

$>{\mathrm{\Delta w}}_i^k={\mathrm{\Delta w}}_i^p-{\mathrm{\Delta w}}_i^q,\forall k=1,···,N---\left(12\right)$>

式中，p和q表示生成第k张差分干涉图的主辅影像的获取时间。由于采用短基线原则，在求解每个时刻对应的残余相位时，会出现秩亏方程组，解决这一问题的办法就是奇异值分解法(SVD)。由此，在奇异值分解的基础上，对M张残余相位进行时域低通滤波处理，以提取最终的非线性形变量，可表示为：

$>{\left(d_i^p\right)}_{\mathrm{NL}\_\mathrm{est}}=\frac{\lambda }{4\pi }·{\left\{\right[{\mathrm{\Delta w}}_i^p{]}_{\mathrm{HP}\_\mathrm{space}}\}}_{\mathrm{TP}\_\mathrm{time}}---\left(13\right)$>

在求解出非线性形变量后，每个相干目标对应的形变序列为：

$>d^p=v_{\mathrm{est}}·(t^p-t^1)+{\left(d_i^p\right)}_{\mathrm{NL}\_\mathrm{est}},k=\mathrm{1,2},···,M---\left(14\right)$>

其中，v_est为解算得到线性形变速率。

3.优点及功效：

本发明提出的一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法与其它领域的InSAR形变监测的原理是一样的。但本发明的特点在于，它无需外部DEM参与，直接利用二维模型同时反演得到地物高程和形变速率，从而可以大大提高高分辨率InSAR在城区地面沉降监测中的精度。

附图说明

图1.条纹密度与垂直基线的关系。对于同一高差，垂直基线越长，条纹越密集。

图2.以Cosmo-skymed3m数据为例，讨论高程误差对形变量的影响。短基线条件下条纹密度降低，有利于形变信息的准确提取。

图3.高分辨率InSAR时序分析反演地物高程与地面沉降量的流程图。

图4.参与InSAR时序分析的干涉图序列基线分布图。其中，垂直基线为300m，时间间隔为1年。

图5(a)为一维模型估计得到的初始DEM，图像为雷达坐标系下旋转90度。

图5(b)为一维模型2次估计得到的高程残差，图像为雷达坐标系下旋转90度。

图5(c)为二维模型初次迭代估计得到的高程修正值，图像为雷达坐标系下旋转90度。

图5(d)为二维模型2次迭代估计得到的高程修正值，图像为雷达坐标系下旋转90度。

图5(e)为二维模型5次迭代估计得到的高程修正值，图像为雷达坐标系下旋转90度。

图5(f)为二维模型9次迭代估计得到的高程修正值，图像为雷达坐标系下旋转90度。

图6.迭代得到的地物高程估计值，图像为雷达坐标系下旋转90度。一维模型高程残差与二维模型多次迭代生成的高程修正值相加得到最终的地物高程估计值，最大值超过300m。其中，黑框位置为滨江花园高层建筑所在位置。

图7.滨江花园高层建筑高度估计值，图像为雷达坐标系下旋转90度。

图8.上海市陆家嘴区域地面沉降速率图，时间间隔：2008/12/10－2010/06/23。

图9.豫园附近某相干目标形变序列，时间间隔：2008/12/10－2010/06/23。

具体实施方式

以Cosmo-skymed InSAR监测上海市陆家嘴区域地面沉降为例，说明本发明在实际工程应中的具体操作步骤。如图3所示，本发明是一种高分辨率InSAR时序分析反演地物高程与地沉降量的方法，该方法具体步骤如下：

步骤一：高分辨率SAR数据选取

选取覆盖上海市主城区的Cosmo－skymed3m分辨率雷达数据，其获取时间为2008年12月至2010年6月，时间跨度为一年半，数据集数量为35景。Cosmo－skymed系统是意大利国防部和空间局共同资助，阿莱尼亚航天公司负责研制的由4颗雷达卫星组成的星座。其3米分辨率数据幅宽40km×40km，方位向和距离向带宽分别为3100MHz和73.5MHz，卫星重访周期为16(4)天。雷达主要数据参数如表1所示，扫描日期详见表2。

表1：选用Cosmo－skymed雷达主要参数表

载波频率9.6000000e+009Hz载波波长3.1cm脉冲带宽9.6000000e+007Hz入射角40°距离向像元大小1.249135m方位向像元大小2.246403m距离向采样率1.2000000e+008Hz雷达扫描模式条带模式数据类型单视复数(SLC)

表2：雷达数据日期表

2008121020090111200902122009022820090316200904012009040920090417200906042009061220090714200908072009081520090924200910022009101020091018200910262009110320091205200912132009122120100122201002072010022320100311201003192010032720100404201004122010042820100506201005142010061520100623

步骤二：超短时空基线干涉像对选取以及初始高程相位估计

考虑到建筑物的高度、干涉图基线及地形相位条纹密度的关系，取研究区平均建筑物高度作为高程约束，确定初次计算地物高程时需要的干涉图序列。这里，根据上海市陆家嘴区域高层建筑物分布密集的特点，确定建筑物平均高度为100m，可知基线为50m时，其条纹密度为0.82整周，100m时为1.64整周。为降低相位解缠的误差以提高初次高程估计的准确性，选择条纹在1个整周以内的干涉图参与运算。确定时间基线小于50天，空间基线小于50m内共34个干涉像对序列作为初始求解的数据集，对每一干涉像对进行干涉处理和相位解缠，并生成相应的相干图和强度图。针对每一干涉图中出现的轨道误差和趋势性干涉条纹，综合利用快速傅里叶变换(FFT)估计和多项式拟合估计的方法予以去除。设定幅度离散指数和相干系数阈值分别为1.5和0.675，筛选得到相干目标候选点。按照一维模型进行地形参数估计，提取如图5a所示的地物高程。

步骤三：InSAR时序分析干涉像对选取和相干目标候选点提取

根据短基线思想构建干涉像对序列，将时间间隔在1个年度内，空间基线小于300m的干涉像对进行差分干涉处理，如图4所示为参与计算的干涉图基线分布，共有187个干涉像对，利用步骤2得到的初始DEM模拟地形相位，去除干涉图中的地形相位，并对每一差分干涉图进行相位解缠。针对每一干涉图中出现的轨道误差和趋势性干涉条纹，综合利用快速傅里叶变换(FFT)估计和多项式拟合估计的方法予以去除，最终得到用于时序分析的初始差分干涉相位图和相干图。在此基础上，利用幅度离散指数(阈值1.5)和相干系数(阈值0.675)分别筛选得到相干目标候选点，并合并得到最终候选点。

步骤四：高程相位校正和线性形变分量估计

选取时间间隔100天内、垂直基线为100m的解缠相位图，重新利用一维模型求解高程残余值。将一维模型求解高程残余值作为二维模型估计的初始高程，分步逐次迭代完成高程和形变信息的求解过程：(1)设定空间基线阈值为100m，初始最大高程修正阈值为45m，将时间基线不超过150天干涉图参与完成初始二维参数估计，其结果如图5c所示；(2)以时间基线步长100天逐步增加干涉图数量直至时间基线为350天，并同时以高程修正步长5m逐步减少最大高程修正阈值，迭代完成二维参数估计；(3)在上述步骤完成的基础上，以空间基线步长100m逐步增加空间基线长度直至达到300m，重复执行(2)和(3)。最终将垂直基线小于300m，时间间隔1个年度以内的干涉图(187个)参与估计最终的高程修正和形变速率。对一维模型求解的高程残余值和逐次求解所得的高程修正值进行相加，求得上海市主城区最终的地物高程估计值，图6所示为本发明估计的陆家嘴高层建筑群区高程。

为了检验建筑物高度估计值的准确性，利用Google Earth三维实景测量数据对楼房高进行比较。选取图6黑框所示成片分布的高层建筑群作为高程估计结果的检验对象，该建筑区对应于世茂滨江花园7栋密集分布的高层楼房(图7)，其平均高度达150-170m，楼层数为45-50层，层间距为3-4m。该建筑群中多栋高程建筑较为均匀的分布，且与其他低矮建筑物相连接，其高度的估计在相位的时序回归分析中具有代表意义，即：在整体上与其他建筑物相连，用于估计高程，同时各自之间又相互独立，在比较高度估计结果中具有代表性。选择与Google三维实景中楼房净高和高程值(相对精度优于1m)进行比较，从Google earth和高程估计中分别提取7座建筑物的楼顶、楼底高程值，计算其净高。表3为世茂滨江花园7栋密集分布的高层建筑实际高程与InSAR估计高程所做的比较，两者的高程互差基本控制在10m以内，说明本发明的地物高程估计值较为可靠。

步骤五：地面沉降速率和形变序列估计

如前所述，在完成初次高程求解后，通过二维迭代回归进行高程和形变速率的逐步改进，得到最终的形变速率结果，在此基础上根据形变与视线向关系将其转换为沉降速率并进行地理编码，图8为上海市陆家嘴区域地面沉降速率图(2008/12/10－2010/06/23)，最大年沉降速率达到57.5mm。对去除高程和线性速率后的残余相位进行时间域和空间域滤波处理，提取最终的非线性形变量。将线性和非线性形变量相加得到每个相干目标的形变序列，图9为豫园附近某相干目标的形变序列。图1为纹密度与垂直基线的关系示意图；图2是以Cosmo-skymed3m数据为例，讨论高程误差对形变量的影响。

表3 建筑物顶底点高差比较