CN104111456A - 一种高速铁路沿线地表形变高分辨率InSAR监测方法 - Google Patents

Abstract

一种高速铁路沿线地表形变高分辨率InSAR监测方法，该方法有四大步骤：步骤一：高分辨率SAR数据选取；步骤二：高分辨率InSAR高速铁路区域形变信息提取方法；步骤三：多轨道形变速率结果集成方法；步骤四：高速铁路线上目标的识别和形变提取方法。本发明能够有效地解决高分辨率InSAR完整覆盖监测高速铁路，以及高速铁路路基沉降与区域地面沉降的识别和分离等问题，可以大大提升我国高速铁路地表形变InSAR精细监测水平。

Description

一种高速铁路沿线地表形变高分辨率InSAR监测方法

技术领域

本发明涉及一种高速铁路沿线地表形变高分辨率InSAR监测方法，属于合成孔径雷达干涉测量技术(InSAR)领域。它能够有效地解决高分辨率SAR单轨影像无法完整覆盖研究区，以及高速铁路路基沉降与区域地面沉降的识别和分离等问题，可以大大提升我国高速铁路地表形变InSAR精细监测水平。

背景技术

雷达干涉测量技术研究的核心目标是测量地物高程(地形测图)和反演雷达视线方向上的形变量等参数。以往InSAR技术多以中等分辨率(空间分辨率为20-30m)雷达数据为主，如ERS、ENVISAT、ALOS和RADARSAT等，其测量结果是分辨率单元内若干个散射体的综合，因而从精细程度上难以实现对大型工程、独立构(建)筑物等地物局部精细测量。相对于中等分辨率SAR而言，高分辨率InSAR是指利用高分辨率(1-3m，如TerraSAR-X、Cosmo-skymed)雷达卫星数据开展地表形变精细监测的InSAR技术，适用于大型构建筑物和典型地物目标的精细监测，极大地拓展了InSAR技术的应用范围。

对于高速铁路形变监测，中等分辨率InSAR监测结果仅能够准确反映区域地表变形场的分布特征，由于受到分辨率的限制，难以实现目标准确定位，因此对于高速铁路路基、轨道及其附属物等稳定程度的监测精细程度不足。而高分辨率InSAR技术有利于直接提取单个目标的干涉相位，能够提供对工程体和周围环境的共同监测，易于分辨环境变化对工程体的影响。目前，针对高速铁路的高分辨率InSAR形变精细监测技术亟待解决以下难点。

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

(2)高速铁路的分布和走向不规则，往往地跨雷达数据数个轨道、多个图幅。以目前在轨的高分辨率SAR卫星为例，其扫描成像覆盖范围为30-50km。只有利用相邻平行轨道高分辨率InSAR联合观测才能实现高速铁路的完整覆盖监测，这涉及到多个相邻轨道高分辨率InSAR观测结果的集成问题，也就是多个轨道处理结果坐标系和参考基准的统一问题。

(3)受自然和人为因素的影响，高速铁路沿线不可避免的会出现地面沉降；同时，列车的高速运行及其震动对路基和桥梁产生的压力会引起路基、桥梁等构筑物的垂向变形。这些因素的共同作用所产生的差异性沉降直接影响着高速列车的运行安全，因此，有效区分高铁沿线区域性地面沉降和铁路桥梁本身的沉降对于工程应用而言意义重大。中等分辨率InSAR区域性沉降监测能准确提取线路周边沉降区和漏斗，明确沿线主要沉降区的分布。针对高铁沿线差异性沉降，需研究高分辨率InSAR监测结果中路基沉降和区域地面沉降的识别与分离。解决这一问题的关键在于如何利用高分辨率InSAR有效提取铁路线上目标形变量和线路周边其它地物的形变量。

本发明针对高分辨率InSAR精细监测高速铁路形变所面临的难点，提供一种解决高速铁路沿线地表形变的精细监测方法。

发明内容

1.目的：本发明的目的是提供一种高速铁路沿线地表形变高分辨率InSAR监测方法。它能够有效地解决高分辨率InSAR完整覆盖监测高速铁路，以及高速铁路路基沉降与区域地面沉降的识别和分离等问题，可以大大提升我国高速铁路地表形变InSAR精细监测水平。

2.技术方案：本发明是一种高速铁路沿线地表形变高分辨率InSAR监测方法，该方法具体步骤如下：

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

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

步骤二：高分辨率InSAR高速铁路区域形变信息提取方法

对高速铁路沿线形变监测的基本思路为：针对雷达图像覆盖范围先进行整体InSAR时序分析解算，进而在形变结果中提取高铁线路及其周边变形信息，进行路基形变分析。这一思路是相对于先单独提取线路上的目标，进而对其相位进行时序分析解算而言的。前者具有整体性，后者为局部处理，出现不可靠性的可能增加。在上述基本思路的指导下，对每一轨道数据集，本发明将以永久散射体干涉测量(PS InSAR)为主，结合短基线干涉测量(SBAS InSAR)技术的优点，综合利用长短基线干涉相位序列的二维相位回归分析，建立以迭代算法为主的多模型高分InSAR数据处理方法，通过逐次迭代求解地物高程和地表形变信息。其详细步骤包括：

(1)研究区原始DEM生成

地形相位补偿是InSAR时序分析处理过程中求取差分相位的基本步骤，解决的途径主要有两种：(ⅰ)利用Lidar、TanDEM-X等获取的DEM数据实现近似同等分辨率差分干涉处理中的地形相位补偿；(ⅱ)利用时间序列分析方法同时解算单个散射体的高程和形变信息。本发明中利用后者实现干涉图地形相位补偿，而原始DEM的生成是完成这一过程的第一步。

①利用时间和空间基线都较短的干涉像对生成初始DEM。

首先，对数据集中的干涉像对进行干涉处理，本发明中综合利用快速傅里叶变换(FFT)估计和多项式拟合估计的方法去除干涉图中的轨道误差和趋势性干涉条纹。其次，利用②提供的方法识别相干目标候选点。最后，针对每一相干目标点，通过多次迭代得到相干目标点处的高程，并插值生成研究区的原始DEM，

②利用点目标识别方法提取相干目标候选点。

本发明中综合采用幅度离散指数(Amplitude Dispersion Index)和相干系数(coherence)来筛选得到相干目标候选点。

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

$D_A=\frac{{\mathrm{\σ}}_A}{m_A}---\left(1\right)$

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

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

$\widetilde{\mathrm{\γ}}=\left|\frac{\frac{1}{N}{\mathrm{\Σ}}_{i=0}^N{M}_i{S}_i^{*}}{\frac{1}{N}{\mathrm{\Σ}}_{i=0}^N{M}_i{M}_i^{*}\frac{1}{N}{\mathrm{\Σ}}_{i=0}^N{S}_i{S}_i^{*}}\right|---\left(2\right)$

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

(2)相干目标形变速率和形变序列生成

①InSAR时序分析数据集选取。

根据短基线(Small Baseline，SBAS)思想构建干涉像对序列，对数据集中的干涉像对进行干涉处理，并利用上述生成的原始DEM模拟地形相位实现地形相位补偿。针对每一差分干涉图中出现的轨道误差和趋势性干涉条纹，综合利用快速傅里叶变换(FFT)估计和多项式拟合估计的方法予以去除，进而生成用于InSAR时序分析的初始差分干涉图、相干图以及所有高分数据的强度图像，解缠每一差分干涉图。综合采用幅度离散指数(Amplitude Dispersion Index)和相干系数来筛选相干目标候选点，以降低高分辨率SAR条件下相干目标数量冗余。

②迭代修正高程误差相位求解形变参数。

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

${\mathrm{\Δ\φ}}_{i,j}^k=[C_B\·B^{\left(k\right)}\·{\mathrm{\Δ\ϵ}}_{i,j}+\frac{4\mathrm{\π}}{\mathrm{\λ}}\·T^{\left(k\right)}\·{\mathrm{\Δv}}_{i,j}]+{\mathrm{\μ}}_{\mathrm{NL}}^{\left(k\right)}+{\mathrm{\α}}^{\left(k\right)}+n^{\left(k\right)}---\left(3\right)$

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

${\mathrm{\φ}}_{\mathrm{mode}l}(i,j,T^{\left(k\right)})=C_B\·B^{\left(k\right)}\·{\mathrm{\Δ\ϵ}}_{i,j}+\frac{4\mathrm{\π}}{\mathrm{\λ}}\·T^{\left(k\right)}\·{\mathrm{\Δv}}_{i,j}---\left(4\right)$

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

${\mathrm{\Δw}}_{i,j}^k={\mathrm{\Δ\φ}}_{i,j}^k-[C_B\·B^{\left(k\right)}\·{\mathrm{\Δ\ϵ}}_{i,j}+\frac{4\mathrm{\π}}{\mathrm{\λ}}\·T^{\left(k\right)}\·{\mathrm{\Δv}}_{i,j}]+{\mathrm{\μ}}_{\mathrm{NL}}^{\left(k\right)}+{\mathrm{\α}}^{\left(k\right)}+n^{\left(k\right)}---\left(5\right)$

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

本发明的InSAR时序分析是在差分干涉图相位解缠的基础上进行的，这时式(3)转换为二维线性函数，可通过建立Delanay三角网或利用冗余网构建更为复杂的连接关系强化对待解算方程组的约束，利用邻近法则将所有距离满足大气相关距离的相干目标连接起来，在求解完成相邻点间的互差后，通过最小二乘或加权平均的方法求解每个目标相对于参考点的高程值和形变速率场。通过多次迭代修正高程误差，完成最终的形变速率估计。对于具有显著非线性形变过程，仍需对残余相位进行更为复杂的处理，以提取非线性形变量。因此，这就需要在上述处理的基础上对残余相位进行时域和空域滤波处理，并利用奇异值分解方法求解非线性形变序列，将线性和非线性分量相加得到每个相干目标的形变序列，完成形变参数求解。

步骤三：多轨道形变速率结果集成方法

(1)坐标系统一

图1所示为相邻平行轨道合成孔径雷达成像模式。多轨道集成的主要难点是解决不同轨道(入射角)下InSAR形变时序分析结果坐标系和参考基准的统一问题。为此，需将雷达图像进行正射校正，使不同轨道下的雷达图像位于统一的地面坐标系下，消除因地形起伏以及成像几何关系引起的畸变等因素的影响。本发明采用唯一的地面坐标系进行坐标系的统一。首先将每个轨道下雷达强度图像进行地理编码，并完成形变速率图的坐标转换。在获得位于同一地面坐标系下的多轨道雷达强度图像后，应用多项式回归完成相邻主辅轨道重叠区域的精确配准，配准函数如式(6)，从而实现相邻轨道的无缝拼接。在此基础上，最终完成多轨道形变速率图的坐标系的统一。

$\begin{array}{c}\mathrm{\Δx}(x,y)=\underset{i=0}{\overset{p}{\mathrm{\Σ}}}\underset{j=0}{\overset{i}{\mathrm{\Σ}}}{a}_{i-j,j}{x}^{i-j}{y}^j\\ \mathrm{\Δy}(x,y)=\underset{i=0}{\overset{p}{\mathrm{\Σ}}}\underset{j=0}{\overset{i}{\mathrm{\Σ}}}{b}_{i-j,j}{x}^{i-j}{y}^j\end{array}---\left(6\right)$

式中：Δx，Δy分别为主辅轨道x，y方向的像对偏移量；p为多项式阶数；a，b为多项式系数。由于所有处理结果已经做了地理编码，相邻轨道图像在同名点只有略微的偏移，利用式(6)进行精确配准时一般只取常数项即可，本发明在相邻主辅轨道重叠区域的精确配准时只取常数项。

(2)参考基准统一

由于入射角的差异，同一相干目标在相邻轨道中的形变速率有所差别；即使在同一轨道下，由于参考基准的差异也会导致不同分块的处理结果产生偏差，如图2(a)、(b)所示。因而，根据主影像上相干目标的空间位置，提取同一坐标系下经过配准的辅影像上对应的相干目标形变参数，进行整体偏差求解。在统计重叠区域的所有同名点后，可按照式(7)和(8)估计相邻轨道之间形变速率的整体偏差，参与统计的相干目标均为满足相关模型的点。在此过程中选择一个中间轨道作为主轨道，其余均相对于该轨道进行基准偏差补偿。

${\mathrm{\Δv}}_{\mathrm{off}}=\frac{1}{N}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}(v_{m_i}-v_{n_i})---\left(7\right)$

${\hat{v}}_{s_i}=v_{s_i}+{\mathrm{\Δv}}_{\mathrm{off}}---\left(8\right)$

式中：Δv_off为轨道间基准偏差；和分别为主辅轨道上相干目标i的形变速率；为改正后辅影像相干目标形变速率。完成整体偏差估计后，对主辅轨道重叠区内相干目标形变参数用式(9)求解：

${\stackrel{\‾}{v}}_i=\frac{P_{m_i}{v}_{m_i}+P_{m_i}{\hat{v}}_{s_i}}{P_{m_i}+P_{s_i}}---\left(9\right)$

其中，为基准统一后相干象元i的形变速率；和分别为权重因子；由相干目标i距离参考位置的距离而定。由于相邻轨道重叠区域结果具有很高的相关性，可利用相关系数公式(10)进行评估，因此本发明在参考基准统一过程中采用先校正主轨道形变速率，再通过线性关系式(11)校正其它轨道。

$r=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(x_i-\stackrel{\‾}{x})(y_i-\stackrel{\‾}{y})}{\sqrt{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^2\·\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\stackrel{\‾}{y})}^2}},\stackrel{\‾}{x}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i,\stackrel{\‾}{y}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i---\left(10\right)$

Sub_M＝a·Sub_S+b   (11)

(3)多轨道集成成果的精度检验

为检验多轨道集成成果的拼接效果，直接利用相邻轨道重叠区域相干目标点进行统计比较，以互差的均方差为统计指标，按照下式计算：

$\mathrm{\σ}=\sqrt{\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{(X_i-Y_i)}^2}{m}}---\left(12\right)$

式中：X_i为降轨观测值，Y_i为升轨观测值(经过基准偏差补偿)，m为统计样本点个数。

步骤四：高速铁路线上目标的识别和形变提取方法

如图3(a)所示，相干目标的分布具有如下特点：

(ⅰ)线上相干目标成对或多对出现，横向间距为5-10m左右，沿高铁线路方向的间距为50m；

(ⅱ)地物高程较之周围地物有差异，普遍超过周围地物3-5m；

(ⅲ)形变特征与周围地物的沉降相似，具有连续性。在沉降较大的地区，线上相干目标随之较大，在一般地区，线上目标形变量普遍较小，且呈现较为随机的变化特征。

如图3(b)所示，与之对应的地物及其特征为：

(ⅰ)相干目标由高铁线上高压线架的组件(横梁、腕臂和立柱)与桥梁板面共同作用，构成散射单元，其反射特征与pole效应有关。

(ⅱ)一组高压线架横向间距设计值为11.3m，沿线路方向间距为49.5m，与相干目标的空间分布特征相同。

(ⅲ)高压线架的高度为2-3m，桥墩高度5-6m，普遍高于地面10m左右。与高分InSAR解算出的高程估计值基本对应。

根据高速铁路线上目标的空间位置、分布密度和特征，结合高速铁路线上地物特点，可以确定相干目标由高铁线上绗架、横梁和立柱(高压线架)与桥梁板面共同作用反射构成(图4)。由于高压线架与桥梁板面固定，可视为刚体，那么由线上相干目标干涉相位获取的形变信息反映的就是高铁线路整体的变化特征。因而，对于路基形变的识别和提取的关键即提取路基上的相干目标。针对路基形变，可根据线路的分布提取沿线沉降场(形变因素)，进而进行高铁线上目标识别，通过沉降序列的变化和周围地物的差异分析变形的主要诱因。

3.优点及功效：

本发明提出的一种高速铁路沿线地表形变高分辨率InSAR监测方法与其它领域的InSAR形变监测的原理是一样的。但本发明的优点在于：(ⅰ)利用高分辨率SAR数据同时解算高程和形变相位，提高了高速铁路形变监测的精度；(ⅱ)统一多个轨道处理结果的坐标系和参考基准，实现了高分辨率InSAR完整覆盖监测高速铁路形变；(ⅲ)通过对高速铁路线上目标的有效识别，能够监测得到高速铁路路基沉降信息。因此，本发明的提出能够大大提升我国高速铁路地表形变InSAR精细监测水平。

附图说明

图1.相邻平行轨道下同一地物对应的相干目标分布示意图。

图2(a)为InSAR地形测量结果参考基准统一修正前干涉图.

图2(b)为InSAR地形测量结果参考基准统一修正后干涉图.

图3(a)为识别的高铁线上相干目标候选点图。

图3(b)为识别的高铁线上相干目标候选点对应的地物目标图。

图4.高速铁路线上目标构成示意图。高铁线上相干目标是由高压线架与桥梁板面共同

作用构成的雷达目标。

图5.高速铁路线上目标高分辨率InSAR形变监测技术流程图。

图6(a)为Track008数据时间和空间基线分布图。

图6(b)为Track009数据时间和空间基线分布图。

图7(a)为Track008处理得到形变速率图。

图7(b)为Track009处理得到形变速率图。

图8.Track008与009重叠区内同名相干目标沉降速率散点统计图。

图9.TerraSAR两个相邻轨道集成得到形变速率图。

图10.京津高速铁路北京段形变速率图。

图11.北京市石材城附近高速铁路形变速率图，其相干目标点叠加在Googleearth上。

具体实施方式

以高分辨率InSAR监测京津高速铁路线上目标形变为例，说明本发明在实际工程应用中的具体操作步骤。图1.相邻平行轨道下同一地物对应的相干目标分布示意图；图2(a)为InSAR地形测量结果参考基准统一修正前干涉图；图2(b)为InSAR地形测量结果参考基准统一修正后干涉图；图3(a)为识别的高铁线上相干目标候选点图；图3(b)为识别的高铁线上相干目标候选点对应的地物目标图；图4.高速铁路线上目标构成示意图；高铁线上相干目标是由高压线架与桥梁板面共同作用构成的雷达目标。如图5所示，本发明是一种高速铁路沿线地表形变高分辨率InSAR监测方法，该方法具体步骤如下：

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

选取德国TerraSAR-X(TSX)星载X波段3米分辨率的2010年11月至2013年11月Track008和009两个相邻轨道的SLC数据。TerraSAR-X3米分辨率数据幅宽30km×60km，带宽150MHz，数据覆盖时间跨距3年，卫星重访周期为11天。TerraSAR-X雷达主要数据参数如表1所示，Track008和009的扫描日期详见表2和3。

表1：选用TerraSAR-X雷达主要参数表

载波频率	9.6500000+9Hz
		载波波长	3.1cm
脉冲带宽	1.5000000e+8Hz
		入射角	26.2542°
距离向像元大小	0.909403m
		方位向像元大小	1.965121m
距离向采样率	1.6482919+8Hz
		雷达扫描模式	条带模式
数据类型	单视复数(SLC)

表2：Track008雷达数据日期表

20101119	20101211	20110102	20110124	20110215	20110309	20110331
							20110422	20110627	20110901	20110923	20111015	20111106	20111128
20111220	20120111	20120202	20120224	20120408	20120430	20120909
							20121012	20121114	20130119	20130221	20130326	20130509	20130531
20130703	20130805	20130907	20131101

表3：Track009雷达数据日期表

20101108	20101130	20101222	20110113	20110226	20110320	20110411
							20110503	20110525	20110616	20110730	20110912	20111026	20111209
20111231	20120122	20120213	20120306	20120328	20120511	20120920
							20121023	20130304	20130520	20130622	20130714	20130816	20130918
20131010

步骤二：高分辨率InSAR高速铁路区域形变信息提取

(1)选取超短时间和空间基线干涉像对序列生成原始DEM

需要说明的是，在地理编码前的所有数据处理都是在雷达坐标系下进行的。针对每一轨道数据集，首先确定时间基线小于50天，空间基线小于50m的干涉像对序列作为求解原始DEM的数据集，对数据集中的干涉像对进行干涉处理得到干涉图、相干图和强度图数据集。利用幅度离散指数(阈值1.65)和相干系数(阈值0.75)分别筛选得到相干目标候选点，并合并得到最终相干目标候选点。对这些相干目标进行多次迭代处理解算得到原始DEM，并将求解所得的DEM作为原始DEM模拟地形相位用以补偿InSAR时序分析干涉图序列中的高程相位。

(2)求解相干目标的形变速率和形变序列

①InSAR时序分析数据集选取。

根据短基线思想构建干涉像对序列，将时间间隔在1个年度内，空间基线小于300m的干涉像对进行差分干涉处理，利用上述所得模拟地形相位补偿干涉图中的地形相位。如图6(a)和(b)所示分别为Track008和009两个轨道参与计算的干涉图基线分布，其中Track008和009满足条件的干涉像对分别为116和125个。针对每一干涉图中出现的轨道误差和趋势性干涉条纹，综合利用快速傅里叶变换(FFT)估计和多项式拟合估计的方法予以去除，最终得到用于时序分析的初始差分干涉相位图。利用幅度离散指数(阈值1.45)和相干系数(阈值0.72)分别筛选得到相干目标候选点，并合并得到最终候选点。

②迭代修正地物高程求解每一轨道的形变参数。

针对每一相干目标点，选取时间间隔100天内、垂直基线为100m的解缠相位图，重新利用一维模型求解高程残余值，将上述高程残余值作为二维参数估计的初始高程。同时将时间基线大于50天，空间基线大于100m的干涉图加入到二维参数估计的序列中，增加用于高程误差估计和形变速率估计的干涉图数目，最终将垂直基线小于300m，时间间隔在1个年度以内的干涉图参与计算，迭代修正地物高程并最终求解相干目标的形变速率。在上述处理的基础上对残余相位进行时域和空域滤波处理，将形变相位中的非线性部分加回线性部分中，可求解每个相干目标的形变序列。图7(a)和(b)所示分别为求取的Track008和009的形变速率图。

步骤三：京津高铁沿线TSX相邻轨道的集成

(1)不同轨道坐标系的统一

InSAR数据处理过程需顾及处理效率和处理能力，对同一轨道采取分块处理，通过上下分块合并得到完整轨道处理结果。在解算得到完整轨道形变速率后，对不同轨道下的形变结果进行地理编码，并利用相邻轨道重叠区域强度图像检校图幅之间的偏差完成精确坐标校正。

(2)不同轨道参考基准的统一

两个轨道重叠区的范围为4km×20km。以Track008处理获取的地面沉降速率为基准(实际解算过程中与地面测量基准相同)，Track009处理结果为待校正结果，对两个轨道下重叠地区内的同名相干目标沉降量进行统计计算，根据式(7)求取二者之间的整体偏差，并比较两组数据的相关程度。重叠区内同名点的个数为83324个像元，两组数据的相关系数为R＝0.9922，说明二组沉降速率线性关系显著。满足如下方程：Sub_T009＝1.011Sub_T009-5.55623，利用上述线性关系式对Track009进行参考基准整体偏差校正，最终将两个轨道结果拼接生成大范围地面沉降速率图，如图9所示为北京东部平原区地面沉降速率图。

(3)集成结果的精度检验

计算集成后Track008和009重叠区域两组数据的互差均方差SD＝3.4795mm，均方差较小，能够满足工作需要。图8所示为重叠区域两组数据的散点统计结果。

步骤四：京津高铁沿线沉降速率提取

根据高铁线路走向分布，在集成地面沉降速率图上提取高铁沿线沉降速率，得到高铁沿线地面沉降速率图，图10所示为京津高铁沉降最为突出的地区。将高分InSAR监测结果叠加于Google earth上，可以精确定位京津高铁主要沉降发生的位置和范围，如图11所示为北京市石材城附近高铁沿线形变情况。

Claims (1)

1.一种高速铁路沿线地表形变高分辨率InSAR监测方法，该方法具体步骤如下：

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

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

步骤二：高分辨率InSAR高速铁路区域形变信息提取方法

对高速铁路沿线形变监测的基本思路为：针对雷达图像覆盖范围先进行整体InSAR时序分析解算，进而在形变结果中提取高铁线路及其周边变形信息，进行路基形变分析；这一思路是相对于先单独提取线路上的目标，进而对其相位进行时序分析解算而言的，前者具有整体性，后者为局部处理，出现不可靠性的可能增加；在上述基本思路的指导下，对每一轨道数据集，将以永久散射体干涉测量即PS InSAR为主，结合短基线干涉测量即SBAS InSAR技术的优点，综合利用长短基线干涉相位序列的二维相位回归分析，建立以迭代算法为主的多模型高分InSAR数据处理方法，通过逐次迭代求解地物高程和地表形变信息；其详细步骤包括：

(1)研究区原始DEM生成

地形相位补偿是InSAR时序分析处理过程中求取差分相位的基本步骤，解决的途径主要有两种：(ⅰ)利用Lidar、TanDEM-X获取的DEM数据实现近似同等分辨率差分干涉处理中的地形相位补偿；(ⅱ)利用时间序列分析方法同时解算单个散射体的高程和形变信息；利用后者实现干涉图地形相位补偿，而原始DEM的生成是完成这一过程的第一步；

①利用时间和空间基线都较短的干涉像对生成初始DEM；

首先，对数据集中的干涉像对进行干涉处理，综合利用快速傅里叶变换估计和多项式拟合估计的方法去除干涉图中的轨道误差和趋势性干涉条纹；其次，利用②提供的方法识别相干目标候选点；最后，针对每一相干目标点，通过多次迭代得到相干目标点处的高程，并插值生成研究区的原始DEM；

②利用点目标识别方法提取相干目标候选点；

综合采用幅度离散指数Amplitude Dispersion Index和相干系数coherence来筛选得到相干目标候选点；

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

$D_A=\frac{{\mathrm{\σ}}_A}{m_A}---\left(1\right)$

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

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

$\widetilde{\mathrm{\γ}}=\left|\frac{\frac{1}{N}{\mathrm{\Σ}}_{i=0}^N{M}_i{S}_i^{*}}{\frac{1}{N}{\mathrm{\Σ}}_{i=0}^N{M}_i{M}_i^{*}\frac{1}{N}{\mathrm{\Σ}}_{i=0}^N{S}_i{S}_i^{*}}\right|---\left(2\right)$

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

(2)相干目标形变速率和形变序列生成

①InSAR时序分析数据集选取；

根据短基线思想构建干涉像对序列，对数据集中的干涉像对进行干涉处理，并利用上述生成的原始DEM模拟地形相位实现地形相位补偿；针对每一差分干涉图中出现的轨道误差和趋势性干涉条纹，综合利用快速傅里叶变换估计和多项式拟合估计的方法予以去除，进而生成用于InSAR时序分析的初始差分干涉图、相干图以及所有高分数据的强度图像，解缠每一差分干涉图；综合采用幅度离散指数和相干系数来筛选相干目标候选点，以降低高分辨率SAR条件下相干目标数量冗余；

②迭代修正高程误差相位求解形变参数；

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

${\mathrm{\Δ\φ}}_{i,j}^k=[C_B\·B^{\left(k\right)}\·{\mathrm{\Δ\ϵ}}_{i,j}+\frac{4\mathrm{\π}}{\mathrm{\λ}}\·T^{\left(k\right)}\·{\mathrm{\Δv}}_{i,j}]+{\mathrm{\μ}}_{\mathrm{NL}}^{\left(k\right)}+{\mathrm{\α}}^{\left(k\right)}+n^{\left(k\right)}---\left(3\right)$

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

${\mathrm{\φ}}_{\mathrm{mode}l}(i,j,T^{\left(k\right)})=C_B\·B^{\left(k\right)}\·{\mathrm{\Δ\ϵ}}_{i,j}+\frac{4\mathrm{\π}}{\mathrm{\λ}}\·T^{\left(k\right)}\·{\mathrm{\Δv}}_{i,j}---\left(4\right)$

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

${\mathrm{\Δw}}_{i,j}^k={\mathrm{\Δ\φ}}_{i,j}^k-[C_B\·B^{\left(k\right)}\·{\mathrm{\Δ\ϵ}}_{i,j}+\frac{4\mathrm{\π}}{\mathrm{\λ}}\·T^{\left(k\right)}\·{\mathrm{\Δv}}_{i,j}]+{\mathrm{\μ}}_{\mathrm{NL}}^{\left(k\right)}+{\mathrm{\α}}^{\left(k\right)}+n^{\left(k\right)}---\left(5\right)$

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

InSAR时序分析是在差分干涉图相位解缠的基础上进行的，这时式(3)转换为二维线性函数，通过建立Delanay三角网或利用冗余网构建更为复杂的连接关系强化对待解算方程组的约束，利用邻近法则将所有距离满足大气相关距离的相干目标连接起来，在求解完成相邻点间的互差后，通过最小二乘或加权平均的方法求解每个目标相对于参考点的高程值和形变速率场；通过多次迭代修正高程误差，完成最终的形变速率估计；对于具有显著非线性形变过程，仍需对残余相位进行更为复杂的处理，以提取非线性形变量，因此，这就需要在上述处理的基础上对残余相位进行时域和空域滤波处理，并利用奇异值分解方法求解非线性形变序列，将线性和非线性分量相加得到每个相干目标的形变序列，完成形变参数求解；

步骤三：多轨道形变速率结果集成方法

(1)坐标系统一

多轨道集成的主要难点是解决不同轨道下InSAR形变时序分析结果坐标系和参考基准的统一问题，为此，需将雷达图像进行正射校正，使不同轨道下的雷达图像位于统一的地面坐标系下，消除因地形起伏以及成像几何关系引起的畸变等因素的影响；采用唯一的地面坐标系进行坐标系的统一；首先将每个轨道下雷达强度图像进行地理编码，并完成形变速率图的坐标转换；在获得位于同一地面坐标系下的多轨道雷达强度图像后，应用多项式回归完成相邻主辅轨道重叠区域的精确配准，配准函数如式(6)，从而实现相邻轨道的无缝拼接；在此基础上，最终完成多轨道形变速率图的坐标系的统一；

$\begin{array}{c}\mathrm{\Δx}(x,y)=\underset{i=0}{\overset{p}{\mathrm{\Σ}}}\underset{j=0}{\overset{i}{\mathrm{\Σ}}}{a}_{i-j,j}{x}^{i-j}{y}^j\\ \mathrm{\Δy}(x,y)=\underset{i=0}{\overset{p}{\mathrm{\Σ}}}\underset{j=0}{\overset{i}{\mathrm{\Σ}}}{b}_{i-j,j}{x}^{i-j}{y}^j\end{array}---\left(6\right)$

式中：Δx，Δy分别为主辅轨道x，y方向的像对偏移量；p为多项式阶数；a，b为多项式系数；由于所有处理结果已经做了地理编码，相邻轨道图像在同名点只有略微的偏移，利用式(6)进行精确配准时一般取常数项即可，在相邻主辅轨道重叠区域的精确配准时只取常数项；

(2)参考基准统一

由于入射角的差异，同一相干目标在相邻轨道中的形变速率有所差别；即使在同一轨道下，由于参考基准的差异也会导致不同分块的处理结果产生偏差，因而，根据主影像上相干目标的空间位置，提取同一坐标系下经过配准的辅影像上对应的相干目标形变参数，进行整体偏差求解；在统计重叠区域的所有同名点后，按照式(7)和(8)估计相邻轨道之间形变速率的整体偏差，参与统计的相干目标均为满足相关模型的点；在此过程中选择一个中间轨道作为主轨道，其余均相对于该轨道进行基准偏差补偿；

${\mathrm{\Δv}}_{\mathrm{off}}=\frac{1}{N}\underset{i=0}{\overset{N-1}{\mathrm{\Σ}}}(v_{m_i}-v_{n_i})---\left(7\right)$

${\hat{v}}_{s_i}=v_{s_i}+{\mathrm{\Δv}}_{\mathrm{off}}---\left(8\right)$

式中：Δv_off为轨道间基准偏差；和分别为主辅轨道上相干目标i的形变速率；为改正后辅影像相干目标形变速率；完成整体偏差估计后，对主辅轨道重叠区内相干目标形变参数用式(9)求解：

${\stackrel{\‾}{v}}_i=\frac{P_{m_i}{v}_{m_i}+P_{m_i}{\hat{v}}_{s_i}}{P_{m_i}+P_{s_i}}---\left(9\right)$

其中，为基准统一后相干象元i的形变速率；和分别为权重因子；由相干目标i距离参考位置的距离而定；由于相邻轨道重叠区域结果具有很高的相关性，利用相关系数公式(10)进行评估，因此在参考基准统一过程中采用先校正主轨道形变速率，再通过线性关系式(11)校正其它轨道；

$r=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(x_i-\stackrel{\‾}{x})(y_i-\stackrel{\‾}{y})}{\sqrt{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(x_i-\stackrel{\‾}{x})}^2\·\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\stackrel{\‾}{y})}^2}},\stackrel{\‾}{x}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{x}_i,\stackrel{\‾}{y}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i---\left(10\right)$

Sub_M＝a·Sub_S+b   (11)

(3)多轨道集成成果的精度检验

为检验多轨道集成成果的拼接效果，直接利用相邻轨道重叠区域相干目标点进行统计比较，以互差的均方差为统计指标，按照下式计算：

$\mathrm{\σ}=\sqrt{\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{(X_i-Y_i)}^2}{m}}---\left(12\right)$

式中：X_i为降轨观测值，Y_i为经过基准偏差补偿后的升轨观测值，m为统计样本点个数；

步骤四：高速铁路线上目标的识别和形变提取方法

相干目标的分布具有如下特点：

(ⅰ)线上相干目标成对或多对出现，横向间距为5-10m左右，沿高铁线路方向的间距为50m；

(ⅱ)地物高程较之周围地物有差异，普遍超过周围地物3-5m；

(ⅲ)形变特征与周围地物的沉降相似，具有连续性；在沉降较大的地区，线上相干目标随之较大，在一般地区，线上目标形变量普遍较小，且呈现较为随机的变化特征；

与之对应的地物及其特征为：

(ⅰ)相干目标由高铁线上高压线架的组件与桥梁板面共同作用，构成散射单元，其反射特征与pole效应有关；

(ⅱ)一组高压线架横向间距设计值为11.3m，沿线路方向间距为49.5m，与相干目标的空间分布特征相同；

(ⅲ)高压线架的高度为2-3m，桥墩高度5-6m，普遍高于地面10m左右；与高分InSAR解算出的高程估计值基本对应；

根据高速铁路线上目标的空间位置、分布密度和特征，结合高速铁路线上地物特点，能确定相干目标由高铁线上绗架、横梁和立柱与桥梁板面共同作用反射构成；由于高压线架与桥梁板面固定，可视为刚体，那么由线上相干目标干涉相位获取的形变信息反映的就是高铁线路整体的变化特征，因而，对于路基形变的识别和提取的关键即提取路基上的相干目标；针对路基形变，根据线路的分布提取沿线沉降场，进而进行高铁线上目标识别，通过沉降序列的变化和周围地物的差异分析变形的主要诱因。