一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法
【专利摘要】一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法，它有五大步骤：步骤一：高分辨率SAR数据选取；步骤二：初始高程相位估计；步骤三：选取InSAR时序分析干涉像对数据集并提取相干目标候选点；步骤四：高程相位校正和线性形变分量估计；步骤五：相干目标的形变序列估计。本发明能够克服城区高分辨率InSAR形变监测数据处理中无高分辨率DEM去除高层建筑物附加相位，有效地将城市高层建筑物在高分辨率InSAR中的附加相位与形变相位进行分离，从而大大提高高分辨率InSAR在城区地面沉降中的监测精度。
【专利说明】-种高分辨率InSAR时序分析反演地物高程与地面沉降量 的方法
【技术领域】
[0001] 本发明涉及一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法，属 于合成孔径雷达干涉测量技术（InSAR)领域。它能够有效地将城市高层建筑物在高分辨率 InSAR干涉图中的附加相位与形变相位进行分离，可以大大提高高分辨率InSAR在城区地 面沉降中的监测精度。
【背景技术】
[0002] 利用InSAR求解高程（DEM)或者形变量，相位解缠是必须的步骤之一，这一步骤在 差分干涉处理中称为解缠，而在时序分析方法中则称为参数估计，其实质是求解缠绕相位 的整周数，基本过程是求解相邻点间的相位差，然后按照特定的路径或网络以一定的约束 条件进行积分，求解观测范围整体的解缠相位，进而反演形变场或高程。相位解缠的前提是 干涉图连续分布且变化平缓，满足相位差小于η的约束条件，整个相位场为无旋场，解缠 结果不随路径而变。而实际上干涉图受噪声或不连续相位（如不连续的陡坡）的影响，难 以满足解缠条件，因而需要按照给定路径求解后再进行连接。高分辨率条件下城市高层建 筑物所引起的相位变化类似于非连续的陡坡相位。对于求解高程而言，条纹密度随基线的 增大而密集，且存在与形变条纹混合的可能。因而，如何同时解算高程和形变相位，提高形 变量估计的精度是高分辨率InSAR应用面临的主要难题。
[0003] 地形相位补偿是InSAR时序分析处理过程中差分相位求取的基本步骤，同时地物 高程也是InSAR时序分析的主要解算结果。中等分辨率SAR数据中表现为一个或者几个点 目标的像元在高分辨率雷达数据中以数个点目标表现出来，每个散射体的散射中心不同， 如一栋建筑物，在中分SAR图像中表现为一个点目标，而在高分SAR图像中以多个不同的散 射体表现出来。对高分辨率InSAR地形相位补偿的最直接方法是获取每个散射中心高精度 的DEM数据，利用Lidar、TanDEM-X等获取的DEM数据实现近似同等分辨率差分干涉处理的 高程相位补偿。然而这种高分辨率的DEM数据很多情况下难以获得。因此，这种通过外部 DEM(如SRTM)补偿整个像元高程的方法在高分辨率数据处理中将难以适用，需要研究单个 目标高程精确计算的方法。
[0004] 在InSAR干涉图中，干涉条纹的频率与基线的关系如图1所不，同等1?程下，基线 越短，条纹越稀疏，利于相位解缠或参数估计。当地形变化为2 π时对应的高程差为：
[0005]
【权利要求】
1. 一种商分辨率InSAR时序分析反演地物商程与地面沉降量的方法，其特征在于：该 方法具体步骤如下： 步骤一：高分辨率SAR数据选取 以TerraSAR - X和COSMO - Skymed为代表的高分辨率雷达卫星系统为开展InSAR 反演城区地物高程和形变精细化监测提供了数据源，高分辨率InSAR相对于中等分辨率 InSAR技术，其总体优势体现在两个方面，S卩（i )高密度相干点目标和短周期4-16天； (ii )对地面点目标的准确定位；高分辨率InSAR数据处理的数据选取要满足空间上尽可 能覆盖完整城区的轨道数据，在时间上数据能够连续接收； 步骤二：初始高程相位估计 InSAR形变时序分析所构建的二维参数估计模型中，考虑到大气的空间相关性，对相邻 两点求差以削弱大气的影响；相干目标i和j差分干涉相位的互差为：
(3) 上式中，CB为与垂直基线相关的系数，T为时间基线，Λ ε为相对高程误差，Λν为相 对形变速率，μ m为非线性形变量，α为大气相位，η为噪声，k表示干涉图个数，与干涉图 序列的组合有关；构建目标函数如下：
(,1) 将上式从相位互差式（3)中减去，得到残余相位为：
(5) 显然，当目标函数的参数△ ε和Λν在准确估计时，残余相位将最小化； 采取将差分干涉图相位解缠的基础上进行的估计，这时式（3)转换为二维线性函数， 通过建立Delanay三角网或利用冗余网构建更为复杂的连接关系强化对待解算方程组的 约束，利用邻近法则将所有距离满足大气相关距离的相干目标连接起来，在求解完成相邻 点间的互差后，通过最小二乘或加权平均的方法求解每个目标相对于参考点的高程值和形 变速率场； 当干涉图时间间隔和空间基线足够小时，变形相位可忽略，则可以将二维模型式（3) 转为一维模型式（6)，进而解算得到DEM ;
(6) (1)选取超短时间和空间基线干涉像对数据集 依据上述理论，根据建筑物的高度、干涉图基线及地形相位条纹密度的关系，选取研究 区平均建筑物高度作为高程约束，确定计算原始地物高程时需要的干涉像对序列；选取时 间基线和空间基线都较短的干涉像对序列作为初始求解的数据集，设定初始时间和空间基 线阈值分别为1\和&，对数据集中的干涉像对进行干涉处理，无需外部DEM进行地形相位 校正，综合利用快速傅里叶变换估计和多项式拟合估计的方法去除干涉图中的轨道误差和 趋势性干涉条纹，生成超短时空基线干涉图数据集、相干图数据集和强度图数据集；对每一 个干涉图进行相位解缠，按照一维模型多次迭代修正高程误差相位估计得到原始DEM ;需 要说明的是，在地理编码之前所有的操作都是在雷达坐标系下进行的； (2) 利用点目标识别方法提取相干目标候选点 针对上述所有生成的干涉图序列，综合采用幅度离散指数即Amplitude Dispersion Index和相干系数coherence来筛选相干目标候选点； 幅度离散指数的计算公式为：
其中，σ A和％分别为像素幅度值的标准差和均值；给定一适当阈值
DA低于阈值的 像元确定为相干目标候选点； 雷达干涉相位图的相干系数估计公式为：
(8) 根据各像元点在相干图中的相干系数序列、和给定的相干系数阈值f，如果 mean
那么则将该像元确定为相干目标候选点； (3) 多次迭代修正高程误差相位估计地物高程 对于每一相干目标点，利用超短时空基线干涉图数据集并按照式（6)所述的一维模型 多次迭代修正高程误差相位，将多次估计的高程误差修正相位相加，得到相干目标点处的 高程，然后，将点矢量转换为栅格数据，生成研究区的原始DEM; 步骤三：选取InSAR时序分析干涉像对数据集并提取相干目标候选点 将空间基线和时间基线不超过某一设定阈值的所有干涉像对进行干涉处理，将步骤二 生成的原始高程模拟地形相位用于所有干涉图的高程相位补偿，生成差分干涉图序列；针 对单一差分干涉图中出现的轨道误差和趋势性干涉条纹，综合利用快速傅里叶变换估计和 多项式拟合估计的方法予以去除，最终生成用于InSAR时序分析的差分干涉图数据集和相 应的相干图数据集以及强度图数据集；对上述所有生成的差分干涉图序列，综合采用幅度 离散指数和相干系数来筛选相干目标候选点，降低高分辨率SAR条件下相干目标数量冗 余； 步骤四：高程相位校正和线性形变分量估计 解缠 InSAR时序分析差分干涉图序列中的每一幅差分干涉图；选取时间和空间基线阈 值分别为T2和B2的干涉像对数据集，重新利用一维模型二次估计高程残余值dhOTi ;将上述 一维模型所得高程残余值dhOTi作为形变时序分析二维参数估计模型的初始高程，以Λ T、 ΛΒ和Adh作为时间、空间基线和高程修正的步长，将Τ2+Λ?^ΡΒ2+ΛΒ内的干涉图参与到 二维参数估计的序列中，设定初始最大高程修正阈值DH max OTi，并初次估计高程修正值和形 变速率；分别以步长ΛΤ和ΛΒ逐步增加用于迭代估算的干涉图数量，并以步长Adh逐步 减小最大高程修正阈值，通过多次迭代处理逐步修正高程估计值和形变速率；直至满足时 间和空间基线以及最大高程修正阈值的所有差分干涉图都参与计算，估计最终的高程改正 和形变速率；由于逐次迭代求解得到的高程改正是对已经去除高程值的估计，因而最终的 (9) 地物高程估计值为逐次迭代修正值之和，即： 式中Η为迭代终止时的地物高程估计值，u为总迭代次数； 步骤五：相干目标的形变序列估计 对于具有显著非线性形变过程，仍需对残余相位进行更为复杂的处理，以提取非线性 形变量；处理的前提仍是基于不同相位的时空频率特性，从差分相位中去除高程和线性速 率后的残余相位为：
(10) 由于已解算出整周相位，因而残余相位为相位主值，其大小满足：
(11) 这里，首先要实现对单个差分干涉图中残余相位的滤波处理；由于大气表现出空间低 频特性，而非线性变形的空间变化范围较小，相对大气相位而言表现为高通特性，因而，对 残余相位图进行空间低通滤波进一步弱化大气的影响；在短基线集条件下求解的非线性 形变量，不同于PSInSAR的处理策略，不直接对大气进行估计，而是通过滤波减弱大气的影 响；残余相位经过空间高通滤波后，大气分量已经减弱，此时，进一步的处理是求解与雷达 数据对应的不同时刻的非线性变形量，包含噪声影响；根据主辅影像的关系，将相位解缠后 的残余相位分解为：
(12) 式中，P和q表示生成第k张差分干涉图的主辅影像的获取时间；由于采用短基线原 贝1J，在求解每个时刻对应的残余相位时，会出现秩亏方程组，解决这一问题的办法就是奇异 值分解法即SVD ;由此，在奇异值分解的基础上，对Μ张残余相位进行时域低通滤波处理，以 提取最终的非线性形变量，表示为：
(13) 在求解出非线性形变量后，每个相干目标对应的形变序列为：
(14) 其中，vest为解算得到线性形变速率。
【文档编号】G01S13/90GK104123464SQ201410353107
【公开日】2014年10月29日 申请日期:2014年7月23日 优先权日:2014年7月23日
【发明者】葛大庆, 王艳, 刘斌, 张玲, 李曼, 郭小方 申请人:中国国土资源航空物探遥感中心