[发明专利]一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法

摘要

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

主权项

一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法，其特征在于：该方法具体步骤如下：步骤一：高分辨率SAR数据选取以TerraSAR－X和COSMO－Skymed为代表的高分辨率雷达卫星系统为开展InSAR反演城区地物高程和形变精细化监测提供了数据源，高分辨率InSAR相对于中等分辨率InSAR技术，其总体优势体现在两个方面，即(ⅰ)高密度相干点目标和短周期4‑16天；(ⅱ)对地面点目标的准确定位；高分辨率InSAR数据处理的数据选取要满足空间上尽可能覆盖完整城区的轨道数据，在时间上数据能够连续接收；步骤二：初始高程相位估计InSAR形变时序分析所构建的二维参数估计模型中，考虑到大气的空间相关性，对相邻两点求差以削弱大气的影响；相干目标i和j差分干涉相位的互差为：Δφi,jk=[CB·B(k)·Δϵi,j+4πλ·T(k)·Δvi,j]+μNL(k)+α(k)+n(k)---(3)]]> 上式中，CB为与垂直基线相关的系数，T为时间基线，Δε为相对高程误差，Δv为相对形变速率，μNL为非线性形变量，α为大气相位，n为噪声，k表示干涉图个数，与干涉图序列的组合有关；构建目标函数如下：φmodel(i,j,T(k))=CB·B(k)·Δϵi,j+4πλ·T(k)·Δvi,j---(4)]]> 将上式从相位互差式(3)中减去，得到残余相位为：Δwi,jk=Δφi,j(k)-[CB·B(k)·Δϵi,j+4πλ·T(k)·Δvi,j]=μNL(k)+α(k)+n(k)---(5)]]>显然，当目标函数的参数Δε和Δv在准确估计时，残余相位将最小化；采取将差分干涉图相位解缠的基础上进行的估计，这时式(3)转换为二维线性函数，通过建立Delanay三角网或利用冗余网构建更为复杂的连接关系强化对待解算方程组的约束，利用邻近法则将所有距离满足大气相关距离的相干目标连接起来，在求解完成相邻点间的互差后，通过最小二乘或加权平均的方法求解每个目标相对于参考点的高程值和形变速率场；当干涉图时间间隔和空间基线足够小时，变形相位可忽略，则可以将二维模型式(3)转为一维模型式(6)，进而解算得到DEM；Δφi,jk=CB·B(k)·Δϵi,j+α(k)+n(k)---(6)]]>(1)选取超短时间和空间基线干涉像对数据集依据上述理论，根据建筑物的高度、干涉图基线及地形相位条纹密度的关系，选取研究区平均建筑物高度作为高程约束，确定计算原始地物高程时需要的干涉像对序列；选取时间基线和空间基线都较短的干涉像对序列作为初始求解的数据集，设定初始时间和空间基线阈值分别为T1和B1，对数据集中的干涉像对进行干涉处理，无需外部DEM进行地形相位校正，综合利用快速傅里叶变换估计和多项式拟合估计的方法去除干涉图中的轨道误差和趋势性干涉条纹，生成超短时空基线干涉图数据集、相干图数据集和强度图数据集；对每一个干涉图进行相位解缠，按照一维模型多次迭代修正高程误差相位估计得到原始DEM；需要说明的是，在地理编码之前所有的操作都是在雷达坐标系下进行的；(2)利用点目标识别方法提取相干目标候选点针对上述所有生成的干涉图序列，综合采用幅度离散指数即Amplitude Dispersion Index 和相干系数coherence来筛选相干目标候选点；幅度离散指数的计算公式为：DA=σAmA---(7)]]>其中，σA和mA分别为像素幅度值的标准差和均值；给定一适当阈值DA低于阈值的像元确定为相干目标候选点； 雷达干涉相位图的相干系数估计公式为：γ~=|1NΣi=0NMiSi*1NΣi=0NMiMi*1NΣi=0NSiSi*|---(8)]]>根据各像元点在相干图中的相干系数序列γi和给定的相干系数阈值如果mean那么则将该像元确定为相干目标候选点；(3)多次迭代修正高程误差相位估计地物高程对于每一相干目标点，利用超短时空基线干涉图数据集并按照式(6)所述的一维模型多次迭代修正高程误差相位，将多次估计的高程误差修正相位相加，得到相干目标点处的高程，然后，将点矢量转换为栅格数据，生成研究区的原始DEM；步骤三：选取InSAR时序分析干涉像对数据集并提取相干目标候选点将空间基线和时间基线不超过某一设定阈值的所有干涉像对进行干涉处理，将步骤二生成的原始高程模拟地形相位用于所有干涉图的高程相位补偿，生成差分干涉图序列；针对单一差分干涉图中出现的轨道误差和趋势性干涉条纹，综合利用快速傅里叶变换估计和多项式拟合估计的方法予以去除，最终生成用于InSAR时序分析的差分干涉图数据集和相应的相干图数据集以及强度图数据集；对上述所有生成的差分干涉图序列，综合采用幅度离散指数和相干系数来筛选相干目标候选点，降低高分辨率SAR条件下相干目标数量冗余；步骤四：高程相位校正和线性形变分量估计解缠InSAR时序分析差分干涉图序列中的每一幅差分干涉图；选取时间和空间基线阈值分别为T2和B2的干涉像对数据集，重新利用一维模型二次估计高程残余值dhori；将上述一维模型所得高程残余值dhori作为形变时序分析二维参数估计模型的初始高程，以ΔT、ΔB和Δdh作为时间、空间基线和高程修正的步长，将T2+ΔT和B2+ΔB内的干涉图参与到二维参数估计的序列中，设定初始最大高程修正阈值DHmax_ori，并初次估计高程修正值和形变速率；分别以步长ΔT和ΔB逐步增加用于迭代估算的干涉图数量，并以步长Δdh逐步减小最大高程修正阈值，通过多次迭代处理逐步修正高程估计值和形变速率；直至满足时间和空间基线以及最大高程修正阈值的所有差分干涉图都参与计算，估计最终的高程改正和形变速率；由于逐次迭代求解得到的高程改正是对已经去除高程值的估计，因而最终的地物高程估计值为逐次迭代修正值之和，即：H=Σi=1uϵi---(9)]]> 式中H为迭代终止时的地物高程估计值，u为总迭代次数；步骤五：相干目标的形变序列估计对于具有显著非线性形变过程，仍需对残余相位进行更为复杂的处理，以提取非线性形变量；处理的前提仍是基于不同相位的时空频率特性，从差分相位中去除高程和线性速率后的残余相位为：Δwi,jk=μi,jk+ai,jk+ni,jk---(10)]]> 由于已解算出整周相位，因而残余相位为相位主值，其大小满足：|Δwi,jk|<π---(11)]]> 这里，首先要实现对单个差分干涉图中残余相位的滤波处理；由于大气表现出空间低频特性，而非线性变形的空间变化范围较小，相对大气相位而言表现为高通特性，因而，对残余相位图进行空间低通滤波进一步弱化大气的影响；在短基线集条件下求解的非线性形变量，不同于PSInSAR的处理策略，不直接对大气进行估计，而是通过滤波减弱大气的影响；残余相位经过空间高通滤波后，大气分量已经减弱，此时，进一步的处理是求解与雷达数据对应的不同时刻的非线性变形量，包含噪声影响；根据主辅影像的关系，将相位解缠后的残余相位分解为：Δwik=Δwip-Δwiq,∀k=1,···,N---(12)]]> 式中，p和q表示生成第k张差分干涉图的主辅影像的获取时间；由于采用短基线原则，在求解每个时刻对应的残余相位时，会出现秩亏方程组，解决这一问题的办法就是奇异值分解法即SVD；由此，在奇异值分解的基础上，对M张残余相位进行时域低通滤波处理，以提取最终的非线性形变量，表示为：(dip)NL_est=λ4π·{[Δwip]HP_space}TP_time---(13)]]> 在求解出非线性形变量后，每个相干目标对应的形变序列为：dp=vest·(tp-t1)+(dip)NL_est,k=1,2,···,M---(14)]]> 其中，vest为解算得到线性形变速率。