[发明专利]一种高速铁路沿线地表形变高分辨率InSAR监测方法

权利要求书

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

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

选取德国TerraSAR-X(TSX)星载X波段3米分辨率的2010年11月至2013年11月Track008和009两个相邻轨道的SLC数据；TerraSAR-X 3米分辨率数据幅宽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雷达数据日期表

2010111920101211201101022011012420110215201103092011033120110422201106272011090120110923201110152011110620111128201112202012011120120202201202242012040820120430201209092012101220121114201301192013022120130326201305092013053120130703201308052013090720131101

表3：Track009雷达数据日期表

2010110820101130201012222011011320110226201103202011041120110503201105252011061620110730201109122011102620111209201112312012012220120213201203062012032820120511201209202012102320130304201305202013062220130714201308162013091820131010

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

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

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

①利用时间和空间基线都较短的干涉像对生成初始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}{\Σ}_{i=0}^N{M}_i{S}_i^{*}}{\frac{1}{N}{\Σ}_{i=0}^N{M}_i{M}_i^{*}\frac{1}{N}{\Σ}_{i=0}^N{S}_i{S}_i^{*}}\right|---\left(2\right)$

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

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

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

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

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

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

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{\μ}}_{NL}^{\left(k\right)}+{\mathrm{\α}}^{\left(k\right)}+n^{\left(k\right)}---\left(3\right)$

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

${\mathrm{\φ}}_{\mathrm{mod}el}(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}^{\left(k\right)}-\[C_B\·B^{\left(k\right)}\·{\mathrm{\Δ\ϵ}}_{i,j}+\frac{4\mathrm{\π}}{\mathrm{\λ}}\·T^{\left(k\right)}\·{\mathrm{\Δv}}_{i,j}\]={\mathrm{\μ}}_{NL}^{\left(k\right)}+{\mathrm{\α}}^{\left(k\right)}+n^{\left(k\right)}---\left(5\right)$

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

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

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

(1)坐标系统一

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

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

$\begin{array}{c}\mathrm{\Δ}x(x,y)=\underset{i=0}{\overset{p}{\Σ}}\underset{j=0}{\overset{i}{\Σ}}{a}_{i-j,j}{x}^{i-j}{y}^j\\ \mathrm{\Δ}y(x,y)=\underset{i=0}{\overset{p}{\Σ}}\underset{j=0}{\overset{i}{\Σ}}{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}}_{off}=\frac{1}{N}\underset{i=0}{\overset{N-1}{\Σ}}(v_{m_i}-v_{n_i})---\left(7\right)$

${\hat{v}}_{s_i}=v_{s_i}+{\mathrm{\Δv}}_{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}{\Σ}}(x_i-\stackrel{\‾}{x})(y_i-\stackrel{\‾}{y})}{\sqrt{\underset{i=1}{\overset{n}{\Σ}}{(x_i-\stackrel{\‾}{x})}^2\·\underset{i=1}{\overset{n}{\Σ}}{(y_i-\stackrel{\‾}{y})}^2}},\stackrel{\‾}{x}=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}{x}_i,\stackrel{\‾}{y}=\frac{1}{n}\underset{i=1}{\overset{n}{\Σ}}{y}_i---\left(10\right)$

Sub_M＝a·Sub_S+b (11)

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

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

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

$\mathrm{\σ}=\sqrt{\frac{\underset{i=1}{\overset{m}{\Σ}}{(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解算出的高程估计值基本对应；

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