CN104122553B - 一种集成多轨道、长条带CTInSAR的区域性地面沉降监测方法 - Google Patents

Abstract

一种集成多轨道、长条带CTInSAR的区域性地面沉降监测方法，它有七大步骤：步骤一：长条带雷达数据分块与配准；步骤二：差分干涉纹图集生成；步骤三：轨道残余相位去除；步骤四：CTInSAR时序分析；步骤五：长条带相邻分块检验与校正；步骤六：多轨道坐标系统统一；步骤七：多轨道信息检验与参考基准统一。本发明解决了InSAR在大区域形变调查与监测中遇到的两个主要问题，一是如何减弱形变信息提取过程中基线误差、大气波动以及相位解缠噪声等在空间域随范围增大而表现明显的低频特性的误差影响；二是如何将相邻轨道的形变监测结果在空间位置上进行精确配准，并统一不同轨道间形变测量的参考基准。

Description

一种集成多轨道、长条带CTInSAR的区域性地面沉降监测方法

技术领域

本发明涉及一种集成多轨道、长条带CTInSAR(Coherent Target SARInterferometry，相干目标合成孔径雷达干涉测量)的区域性地面沉降监测方法，属于合成孔径雷达干涉测量(InSAR)技术应用领域。它克服了大区域CTInSAR数据处理过程中误差传播与基准偏差的问题，达到了区域性地面沉降的高精度调查与监测。

背景技术

当前我国存在大范围的地面沉降，华北平原、长江三角洲、汾渭谷地三大主要沉降区的影响范围超过数万平方公里，且有进一步扩展连片分布的特点。其他地区如江汉平原、松嫩平原、下辽河以及珠江三角洲等也有不同程度的地面沉降发生。传统的主要地面沉降监测手段(水准测量、基岩标、分层标和GPS监测网等)均面临空间覆盖有限、监测频率低、监测点密度低、观测时间长且难以同步等问题。星载InSAR技术是根据星载雷达图像的相位数据来获取地面形变信息的一门技术，能满足调查监测精度要求的同时，具有大覆盖范围、快速重访、高监测点位密度等特点，且易于组织，人为干扰少。

大区域地面沉降场调查与监测要求雷达数据具备二维方向上的完整覆盖，而单个图幅重复轨InSAR主要解决局部地区(几百至几千平方公里)地表形变场的高精度连续监测，以常用的ENVISAT ASAR扫描模式的数据来说，单个图幅覆盖范围在100km×100km，涉及几万平方公里大区域的地面沉降监测必然需要多轨道、长条带的InSAR成果集成。多轨道集成包含相邻平行轨道集成(如图1所示)以及升降轨模式集成，其中相邻平行轨道集成是解决大范围覆盖的主要技术途径。

多轨道、长条带集成方法基于单图幅CTInSAR时序分析技术，实现大区域地面沉降监测。但随着空间覆盖范围的增大，InSAR数据处理过程中误差因素也随之增大。这些误差主要包括基线误差、大气波动以及相位解缠噪声(相位不连续)等在空间域表现为低频特性的因素。随着覆盖范围的增大，低频信号的累加效应表现明显，影响了沉降场监测结果的准确性，往往使得沉降场观测结果中附加某种趋势信号。

另一个问题是不同轨道(入射角)下InSAR成果具有不同的坐标系和参考基准。无论采用何种时序分析方法，每个轨道下的形变监测结果都存在于独立的坐标系中(雷达坐标系)，其形变场参考基准也不相同。不同轨道下SAR数据均与各自独立的主图像进行精确配准，因而，其观测结果位于不同主图像的坐标系内。斜距成像的特点决定了不同入射的雷达图像存在不同的变形(扭曲、旋转等)，特别是在地形起伏较大的地区。即使同一景雷达图像内，其象元的分辨率也存在差别。此外，InSAR技术所获取的形变量均是相对于某一个特定的参考基准而言的，其参考基准是在时序分析过程中选择的某一特定相干目标，其活动状况往往未知。这使得不同轨道下的形变场存在参考基准偏差，涉及到整体偏差补偿问题。

因此，为解决区域性地面沉降InSAR监测中存在的问题，本专利通过将长条带雷达数据分块，各分块独立完成CTInSAR时序分析，获取各分块高精度地面沉降监测成果，有效地降低了长条带合成孔径雷达数据的干涉处理误差，并通过对同一条带不同分块和多轨道不同条带下形变监测结果的检验与校正，最终达到大区域地面沉降的全覆盖、同时段调查与监测。

发明内容

1、目的

本发明的目的是为了获取全覆盖、同时段的高精度区域性地面沉降信息，解决了InSAR在大区域形变调查与监测中遇到的两个主要问题，一是如何减弱形变信息提取过程中基线误差、大气波动以及相位解缠噪声等在空间域随范围增大而表现明显的低频特性的误差影响；二是如何将相邻轨道的形变监测结果在空间位置上进行精确配准，并统一不同轨道间形变测量的参考基准。

2、技术方案

本发明是一种集成多轨道、长条带CTInSAR的区域性地面沉降监测方法，建立了满足低相干条件下区域性地面沉降监测的CTInSAR技术，方法具体步骤如下：

步骤一：长条带雷达数据分块与配准

符合InSAR测量要求的星载SAR数据通常以长条带方式获取，数据通常为单视复数数据(single look complex digital image，SLC)，记录了目标在影像上的散射强度以及雷达回波的相位信息。

假设同一长条带获得m个重轨数据，在对长条带数据进行处理时，首先选择一幅主影像M，对其进行方位向上的分块裁切，得到n个分块分别为M_frame1、M_frame2、……、M_framen，每个分块称为一个Frame，每个Frame在雷达影像的方位向尽量与距离向保持相同地面距离范围，相邻分块之间保留20％的重叠区域。

将其他SAR数据做为辅影像分别以M_frame1、M_frame2、……、M_framen为主影像进行精密配准，得到n个frame配准后的SLC数据集[M,S₁,S₂,……,S_m-1]。针对每个frame的SLC数据集执行步骤二至步骤四。

步骤二：差分干涉纹图集生成

考虑空间基线与时间基线对SAR图像对造成的失相干影响，选择具有较短空间基线与时间基线的图像对两两进行组合，生成差分干涉纹图集。

每个图像对主副影像共轭相乘得到干涉纹图，干涉纹图的相位组成包括：

φ_int＝δφ_flat+δφ_topo+δφ_mov+δφ_atm+δφ_noise (1)

式中：δφ_flat为斜距变化产生的平地相位，δφ_topo为地形相位，δφ_mov是两次获取SAR图像时地表移动引起的沿雷达视线向形变相位，δφ_atm为大气波动引起的相位变化，δφ_noise为噪声相位。

根据已知参数模拟平地相位与地形相位：

${\mathrm{\δ\φ}}_{flat}=\frac{4\mathrm{\π}}{\mathrm{\λ}r}\·\frac{b\mathrm{\Δ}r}{\tan \mathrm{\θ}}---\left(2\right)$

${\mathrm{\δ\φ}}_{topo}=\frac{4\mathrm{\π}}{\mathrm{\λ}r}\·\frac{b\mathrm{\Δ}h}{sin\mathrm{\θ}}---\left(3\right)$

其中：λ为雷达波长，r为雷达到目标之间的斜距，b为垂直基线长，Δr和Δh分别为像元间的斜距和高程增量。

去除平地、地形相位得到差分干涉相位：

φ_diff＝δφ_mov+δφ_toperror+δφ_atm+δφ_noise (4)

式(4)中δφ_toperror为差分干涉纹图中由于外部DEM误差带来的相位。

步骤三：轨道残余相位去除

地形相位与平地相位的计算均与基线长度相关，而轨道的精度直接影响干涉相对空间基线的计算。由于轨道的不准确性形成的残余相位，即Orbit ramp或Orbit residual。应用精密轨道参数修正初始轨道参数可以有效降低这种因素的影响，但当干涉图范围扩大时，二维平面上轨道残余相位的积累将使得形变图上附加一种趋势，这种趋势随着数据覆盖范围的增大而变化。图2所示为沉降场附加轨道残余相位的剖面。

本发明中基于最小二乘进行趋势面拟合进而去除轨道残余相位。首先依据先验知识挑选出差分干涉纹图集中有明显残余相位的干涉纹图，针对有误差的干涉纹图首先进行相位解缠，通过提取解缠相位图中一定采样间隔的解缠相位，利用二次或者高次函数进行相位回归，进而拟合趋势面。这里选择二次函数(式5)作为拟合函数进行相位回归。

式中：为采样点解缠相位；(x,y)为采样点图像坐标；a_i为多项式系数。基于最小二乘原则，计算出a_i，从差分相位φ_diff中减去拟合的轨道残余相位。

步骤四：CTInSAR时序分析

当SLC数据集中数据量较少时(m<20)，依据相干目标频谱变化小这一特性进行相干目标的选择。点目标的实际几何尺寸往往小于雷达分辨单元，与分布目标具有不同的信号特征，不具备斑点噪声的特性，相位稳定，因而可将其视为相干目标。由于点目标稳定的后向散射特性，在获取原始信号过程中各个子视图像上表现出相似的后向散射特征，不同距离向和方位向的子带宽内的散射强度大致相同。因而可将SLC图像作反变换，生成若干个子视图(sub-looks)，将子视图作谱相关分析，根据其相关性的大小来识别点目标。假设各子视内的某像元的信号x_i(n)＝a_ig(n)，1≤i≤L为视数，g(n)为该像元信号，按照能量比标准衡量L个子视信号的整体相似性为：

$E=\frac{{\left(\underset{i=1}{\overset{L}{\mathrm{\&Sigma;}}}{a}_i\right)}^2}{L\&CenterDot;\underset{i=1}{\overset{L}{\mathrm{\&Sigma;}}}{a_i}^2},(0\&le;E\&le;1)---\left(6\right)$

当a₁＝a₂＝…＝a_L时，E＝1。能量比越大，整体相似性越好，每幅SLC中频谱特征统计值E高于某一给定阈值的像元都被识别为点目标候选点，在候选点中将雷达目标的散射强度值高于某一阈值的点作为相干点。

利用相干点将去除了轨道残余相位的差分干涉纹图集的相位进行提取，得到差分干涉相位的点集。相位组成如式(4)所示，其中形变相位δφ_mov包含线性形变δφ_linear和非线性形变δφ_nonlinear两部分。差分干涉相位可以表达为

$\begin{array}{c}{\mathrm{\&phi;}}_{diff}={\mathrm{\&delta;\&phi;}}_{linear}+{\mathrm{\&delta;\&phi;}}_{toperror}+{\mathrm{\&delta;\&phi;}}_{nonlinear}+{\mathrm{\&delta;\&phi;}}_{atm}+{\mathrm{\&delta;\&phi;}}_{noise}\\ =\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\&CenterDot;v\&CenterDot;T+\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}r}\&CenterDot;\frac{b\mathrm{\&epsiv;}}{\sin \mathrm{\&theta;}}+{\mathrm{\&delta;\&phi;}}_{nonlinear}+{\mathrm{\&delta;\&phi;}}_{atm}+{\mathrm{\&delta;\&phi;}}_{noise}\end{array}---\left(7\right)$

式(7)中，λ为雷达波长，v为线性形变速率，T为生成差分干涉纹图的图像对时间基线，r为雷达到目标之间的斜距，b为垂直基线长，ε为高程误差，θ为雷达入射角。

由于在整个干涉纹图集中每个相干点的线性速率和DEM误差均为常量，是模型中的线性分量，利用以下二维模型可以很好地实现形变速率和高程改正的估计：

$\begin{array}{c}{\mathrm{\&delta;\&phi;}}_{\mathrm{mod}el}(x_m,y_m,x_n,y_n,T_i)=\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\&CenterDot;T_i\&CenterDot;\&lsqb;v_{\mathrm{mod}el}(x_m,y_m)-v_{\mathrm{mod}el}(x_n,y_n)\&rsqb;\\ +\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\&CenterDot;\frac{b\left(T_i\right)}{r\left(T_i\right)\&CenterDot;\sin \left({\mathrm{\&theta;}}_i\right)}\&lsqb;{\mathrm{\&epsiv;}}_{\mathrm{mod}el}(x_m,y_m)-{\mathrm{\&epsiv;}}_{\mathrm{mod}el}(x_n,y_n)\&rsqb;\end{array}---\left(8\right)$

T_i为第i个干涉像对的时间基线。对于相干目标，由于差分干涉纹图中存在不同的相位补偿，单个像元的相位值在处理过程中并没有实际意义。在确定相位解缠参考点后连接相邻两个相干目标(x_m,y_m)和(x_n,y_n)，通过式(9)求时序模型函数的γ_model最大值可以得到式(8)中v和ε的估算结果v_estimated和ε_estimated。

${\mathrm{\&gamma;}}_{\mathrm{mod}el}(x_m,y_m,x_n,y_n)=\frac{1}{N}\&CenterDot;|\underset{i=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\exp \&lsqb;j\&CenterDot;\left({\mathrm{\&delta;\&phi;}}_{diff}\right(x_m,y_m,x_n,y_n,T_i)-{\mathrm{\&delta;\&phi;}}_{\mathrm{mod}el}(x_m,y_m,x_n,y_n,T_i\left)\right)\&rsqb;|---\left(9\right)$

式(9)中N为干涉纹图集中包含的干涉纹图个数。

去除线性分量，模型残余相位可用下式表示：

$\begin{array}{c}{\mathrm{\&delta;\&phi;}}_{residual}(x,y,T_i)={\mathrm{\&delta;\&phi;}}_{nonliear}(x,y,T_i)+{\mathrm{\&delta;\&phi;}}_{atm}(x,y,T_i)+{\mathrm{\&delta;\&phi;}}_{noise}(x,y,T_i)\\ ={\mathrm{\&delta;\&phi;}}_{diff}(x,y,T_i)-{\mathrm{\&delta;\&phi;}}_{linearest\mathrm{mod}el}(x,y,T_i)\\ ={\mathrm{\&delta;\&phi;}}_{diff}(x,y,T_i)-\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\&CenterDot;T_i\&CenterDot;v_{estimated}-\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\frac{b}{r\&CenterDot;\sin \mathrm{\&theta;}}\&CenterDot;{\mathrm{\&epsiv;}}_{estimated}\end{array}---\left(10\right)$

式中：δφ_diff为初始差分相位，δφ_{linearestmodel}为线性模型将线性形变和DEM误差合成得到的相位。模型残余相位中包含大气相位影响δφ_atmos和非线性移动δφ_nonliear(x,y,T_i)以及失相干等因素造成的噪声δφ_noise(x,y,T_i)。

利用1km×1km二维移动平均窗口进行空间低通滤波，滤波完成后，残余相位的噪声δφ_noise(x,y,T_i)被滤除，最后，利用时域高通滤波估计大气贡献δφ_atm，区分出非线性形变相位δφ_nonliear(x,y,T_i)。

利用前面估计出的线性分量与非线性分量获得完整的形变演变量。

$\mathrm{\&rho;}\left(t_i\right)=v_{estimated}\&times;(t_i-t_0)+\frac{\mathrm{\&lambda;}}{4\mathrm{\&pi;}}\&times;{\mathrm{\&phi;}}_{nonlinear}\left(t_i\right)---\left(11\right)$

依据上述相干目标时序分析方法，通过相干目标时序分析模型迭代补偿DEM误差、通过滤波减小大气相位以及噪声的影响，获取到各相干目标点的地面沉降量。

步骤五：长条带相邻分块检验与校正

对长条带中的每个frame重复步骤二至步骤四，得到条带下各分块地面沉降信息。长条带的每个frame获取的形变信息图的参考基准位置不同，需要对连续frame的参考基准进行整体偏差补偿。整体偏差估算依赖于重叠部分相干目标形变参数的统计分析，依据同名点按照式(12)统计连续frame之间形变速率的整体偏差，参与统计的相干目标均为满足相关模型的点。在此过程中选择连续frame中的一个作为基准，其余均相对于该frame进行基准偏差补偿。一般情况下，在每个frame时序分析过程中，根据时间基线的变化关系，将在不同时间间隔下稳定区域内的相干目标作为参考基准，以此降低frame之间形变参数的整体偏差。以两块连续frame之间的检验、校正为例，选择其中一个frame作为基准(以先验知识为基础进行主frame的选择)，依据frame之间重叠区域同名点的统计信息进行基准偏差的整体补偿

${\mathrm{\&Delta;\&rho;}}_{off}=\frac{1}{N}\underset{0}{\overset{N-1}{\mathrm{\&Sigma;}}}({\mathrm{\&rho;}}_{mi}-{\mathrm{\&rho;}}_{si})---\left(12\right)$

${\widehat{\mathrm{\&rho;}}}_s={\mathrm{\&rho;}}_s+{\mathrm{\&Delta;\&rho;}}_{off}---\left(13\right)$

式(12)和(13)中：Δρ_off为frame之间基准偏差；ρ_mi为基准frame重叠区中相干目标i的形变量；ρ_si为需调整frame重叠区中相干目标i的形变量；N表示两个frame包含的同名点个数；为调整的frame中所有相干目标改正后的形变向量。

对于重叠区域内的各相干点目标形变量ρ_i的取值采用距离加权平均方法求解：

${\mathrm{\&rho;}}_i=P_{mi}\&times;{\mathrm{\&rho;}}_{mi}+P_{si}\&times;{\widehat{\mathrm{\&rho;}}}_{si}---\left(14\right)$

$P_{mi}=\frac{D_{mi}}{D},P_{si}=\frac{D_{si}}{D}---\left(15\right)$

式中P_mi和P_si分别为基准frame和调整frame第i个相干目标的权，ρ_mi和分别为重叠区内第i个相干目标在基准frame和调整frame中的形变量，D为重叠区长度，D_mi和D_si分别为基准frame和调整frame中相干目标距离重叠区内各自影像方位向边缘的距离(图3a、b)。

步骤六：多轨道坐标系统统一

通过对不同轨道下的多个长条带InSAR地面沉降信息集成进而获取大覆盖范围的监测成果，首先需要对不同的轨道数据进行坐标系统一，具有相同的空间坐标基准。

首先需要获取研究区的DEM数据，根据外部DEM将雷达图像与监测成果进行正射校正，即地理编码，使得不同轨道下的雷达图像位于相同的地面坐标系下，消除因地形起伏以及成像几何关系引起的畸变等因素的影响。高程精度优于10m的DEM数据保证了平坦地区配准精度优于一个像元。依据卫星参数将DEM模拟至SAR坐标系，建立初始坐标变换表，再将DEM模拟的影像与SAR图像进行精密配准改进坐标变换表。通过反向变换，将雷达坐标系下各轨道地面形变信息地理编码至DEM所在的地面坐标系内。

步骤七：多轨道信息检验与参考基准统一

位于相同地面坐标系的不同轨道形变信息集成过程中，受制于雷达成像入射角的影响，不能像长条带连续frame重叠区内的相干目标在相邻frame中基本一致，相邻轨道下相干目标是完全不同的(图4)。

选择其中一条轨道为主轨道，相邻另一条为辅轨道，首先通过反距离加权插值方法将辅轨道上相干目标分析结果进行空间域插值为面，设辅轨道重叠区内共有n个相干目标点，其插值后平面各点坐标(x,y)的沉降速率值通过下式计算得出：

$\mathrm{\&rho;}(x,y)=\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{{\mathrm{\&rho;}}_i}{d_i^2}}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{1}{d_i^2}}---\left(16\right)$

式中ρ_i为第i个相干目标点的沉降量，d_i为点(x,y)到第i个相干目标点的距离。再利用主轨道上重叠区内相干点目标位置提取其在辅轨道上插值面上的形变值ρ_s。

根据相同相干点目标在主轨道的沉降值ρ_m和在辅轨道的沉降值ρ_s进行检验与校正，检验与校正方法同长条带相邻分块的检验与校正。

3、优点及功效：本发明通过对长条带、多轨道CTInSAR分析处理，达到了大区域地面沉降的高精度调查与监测。

(1)长条带数据的分块处理，减小了雷达数据常规干涉测量中轨道、大气、相位解缠等误差的传播影响。针对每个分块利用相干目标时序分析方法获取高精度形变信息。

(2)利用长条带连续frame之间重叠区内的地面沉降信息，完成多frame之间的检验与校正，保证了长条带地面沉降信息的完整性与正确性。

(3)通过对长条带数据进行地理编码达到了不同轨道数据的空间基准统一。

(4)利用相邻轨道长条带间重叠区内的相干目标点，通过插值的方法，完成了多轨道的形变信息检验与校正，最终实现了对大区域地面沉降的高精度监测。本方法对我国经济快速发展过程中区域性地面沉降调查与监测，具有普适性。

附图说明

图1：大区域地面沉降InSAR多轨道测量。多轨道长条带数据覆盖大区域范围，通过对这些数据的处理达到大区域地面沉降调查与监测的需要。

图2：残余相位对形变场的影响。沿误差传播方向做剖面，可以看到残余相位为形变量求解带来的趋势性影响，当干涉图范围扩大时，二维平面上轨道残余相位的积累将使得形变图上附加一种明显的趋势，这种趋势随着数据覆盖范围的增大而增大。

图3a：为长条带相邻frame重叠区内相干目标确定权重的距离示意图，D_m为相干目标距离基准frame重叠区边缘的距离，D_s为相干目标距离调整frame重叠区边缘的距离。

图3b：为相邻轨道长条带重叠区内相干目标确定权重的距离示意图，D_m为相干目标距离主轨道重叠区边缘的距离，D_s为相干目标距离辅轨道重叠区边缘的距离。

图4：多轨道相邻长条带重叠区内相干目标位置差异示意图。对于相邻的条带，由于雷达成像入射角的影响，相邻轨道重叠区内相干目标位置是完全不同的，图中十字符号为主轨道相干目标，圆圈符号为辅轨道相干目标。

图5：轨道残余相位去除示意图

图6：Track-218与Track-447重叠区相干目标分布示意图。

图7：Track-218与Track-447重叠区相干目标沉降量互差统计直方图。

图8：多轨道、长条带CTInSAR监测成果示意图。图为京津沉降区2008年度InSAR监测地面沉降速率图，图中散点(相干目标)由浅到深表明沉降愈加严重。

图9：本发明数据处理流程图。

具体实施方式

以京津主要沉降区为例，说明本发明在实施中的具体步骤。

实施过程中利用完整覆盖京津主要沉降区的ENVISAT雷达数据(IS2，stripmode)，涉及三个轨道，Track-218、Track-447和Track-175，三个轨道的相邻关系如图1所示，每个条带的处理长度均为200km，宽度为100km，条带间横向重叠度约为40％，获取的雷达数据时间从2007年12月截止2009年1月。所用的雷达数据获取时间如表1：

表1雷达数据列表

	Track-218	Track-447	Track-175
				1	20071219	20080104	20080120
2	20080123	20080208	20080330
				3	20080227	20080314	20080504
4	20080402	20080418	20080713
				5	20080507	20080523	20080817
6	20080611	20080627	20080921
				7	20080716	20080801	20081026
8	20080820	20080905	20090104
				9	20080924	20081010
10	20081029	20081219
				11	20090107	20090123

实施过程利用的外部DEM为SRTM DEM，轨道参数为doris精密轨道。

见图9，本发明是一种集成多轨道、长条带CTInSAR,的区域性地面沉降监测方法，该方法具体步骤如下：

步骤一：长条带雷达数据分块与配准

长条带雷达数据配准与分块以Track-218为例，根据获取数据的实际长度，以Track-218获取的第一景影像20071219为主影像，将其裁切为2个frame M_patch1和M_patch2，两个Frame在雷达影像的方位向与距离向分别保留25000和5000像元，多视后方位向与距离向的宽度均为5000像元，距离约为100km，两个frame之间保留20％的重叠区域。

将Track218下的其他SAR数据分别以M_patch1、M_patch2为主影像进行精密配准，得到2个分块配准后的SLC数据集。

针对三个条带下的每个分块执行步骤二至步骤四。

步骤二：差分干涉纹图集生成

首先利用DORIS精密轨道进行SLC数据轨道修正。选择空间基线小于350m，时间基线大于60天小于400天的图像对进行组合，Track-218、Track-447和Track-175三个轨道的各frame分别形成35个、34个和21个干涉图像对组合，每个图像对主副影像共轭相乘得到干涉纹图集。根据已知参数及式(2)计算平地相位，从干涉纹图中去除平地相位。依据卫星参数将STRM DEM模拟至SAR影像坐标系，建立初始坐标变换表，再将DEM模拟SAR影像与SAR图像通过地形特征进行精密配准，进而改进坐标变换表，将STRM DEM通过坐标变换表变换到SAR影像坐标系，利用式(3)计算地形相位。从干涉纹图减去平地相位与地形相位得到差分干涉纹图。

步骤三：轨道残余相位去除

依据先验知识判断轨道残余相位的影响，确定需要进行残余相位去除的差分干涉纹图。Track-218下的干涉纹图集中有24个干涉纹图进行了残余相位的去除。对于每个需要进行残余相位去除的差分干涉纹图进行相位解缠，通过提取解缠相位图中一定采样间隔的解缠相位，依据式(3)利用二次函数进行相位回归，拟合趋势面。从差分相位图中减去拟合的趋势面，得到修正了轨道残余相位差分干涉纹图，利用去除了残余的差分干涉纹图进行下一步处理将避免产生附加的形变趋势，如图2。如图5为轨道残余相位去除效果示意。

步骤四：CTInSAR时序分析

依据点目标识别方法，将E阈值取0.4，强度值阈值取图像平均强度的10倍，得到点集。用点位提取差分干涉相位得到差分干涉相位点集pdiff，Track-218各分块的pdiff就是包含35个干涉图像对相位的点集，针对pdiff进行相干目标时序分析，获取到相干目标点的地面沉降量(包含线性与非线性沉降)。

对3个轨道下每2个分块执行步骤二至步骤四，得到3个条带中每个分块的所覆盖区域的地面沉降信息。

步骤五：长条带相邻分块检验与校正

对3个轨道下的每2个frame的解算参考基准进行整体偏差补偿。以Track-218为例，其中第一个frame包含北京市区，依照先验知识，北京市中心稳定，在相干目标时序分析时选择市中心一点做为解算起点，将这个frame做为基准frame。参考基准整体偏差估算依据重叠部分相干目标形变参数的统计分析，依据同名点按照式(12)统计两个frame之间形变速率的整体偏差为5.26365mm，将位于下方的frame所有相干目标的形变量调整-5.26365mm，对于重叠区内的各相干点目标利用公式(14)进行计算，各相干点权重确定如图3a，得到Track-218的长条带地面沉降信息。

步骤六：多轨道坐标系统统一

通过坐标变换表将3个轨道下的长条带SAR图像及地面沉降信息反向变换到STRMDEM所在的WGS-84坐标系内。

步骤七：多轨道信息检验与参考基准统一

以Track-218为基准调整Track-447，Track-218的北京地区和Track-447的天津地区北部的重叠部分，覆盖范围约为40km×100km，包含相干目标167851个，如图6，放大到局部，两相邻轨道的相干目标差异如图4所示。依据式(16)对Track-447的相干目标插值生成面数据，再由Track-218位于重叠区的相干目标进行沉降量提取，得到重叠区相干目标沉降量差值的统计分布如图7所示，显示差值均值为-0.09408mm，差值范围为[-13.0496mm,12.760325mm]，位于[-5mm,5mm]以内的相干目标占总体样本的96.63％。将Track-447下所有相干点目标的沉降量值减去0.09408mm，对于Track-218与Track-447重叠区内的相干点目标取值依式(14)进行计算，各相干点权重确定如图3b。

同样的方法以调整后的Track-447为基准调整Track-175轨道下的沉降量信息，最终得到京津主要沉降区2007年12月到2009年01月的全覆盖、相同时间段的沉降量信息，如图8为京津主要沉降多轨道长条带InSAR监测地面沉降结果图。具体实施过程中获取了天津地区和北京地区的地面沉降水准测量值，视水准观测为沉降量的最或然值时，通过对两组观测值的检验，得到InSAR测量获取的形变量的点位中误差为±5mm。

由此，本发明的方法满足区域性地面沉降调查与监测的测量要求，能够为我国区域性地面沉降的防治工作提供可靠的基础数据。

Claims (1)

1.一种集成多轨道、长条带CTInSAR的区域性地面沉降监测方法，其特征在于：该方法具体步骤如下：

步骤一：长条带雷达数据分块与配准

符合InSAR测量要求的星载SAR数据以长条带方式获取，数据为单视复数数据即SLC，记录了目标在影像上的散射强度以及雷达回波的相位信息；

设同一长条带获得m个重轨数据，在对长条带数据进行处理时，首先选择一幅主影像M，对其进行方位向上的分块裁切，得到n个分块分别为M_frame1、M_frame2、……、M_framen，每个分块称为一个Frame，每个Frame在雷达影像的方位向尽量与距离向保持相同地面距离范围，相邻分块之间保留20％的重叠区域；

将其他SAR数据做为辅影像分别以M_frame1、M_frame2、……、M_framen为主影像进行精密配准，得到n个frame配准后的SLC数据集[M,S₁,S₂,……,S_m-1]；针对每个frame的SLC数据集执行步骤二至步骤四；

步骤二：差分干涉纹图集生成

考虑空间基线与时间基线对SAR图像对造成的失相干影响，选择具有较短空间基线与时间基线的图像对两两进行组合，生成差分干涉纹图集；

每个图像对主副影像共轭相乘得到干涉纹图，干涉纹图的相位组成包括：

φ_int＝δφ_flat+δφ_topo+δφ_mov+δφ_atm+δφ_noise (1)

式中：δφ_flat为斜距变化产生的平地相位，δφ_topo为地形相位，δφ_mov是两次获取SAR图像时地表移动引起的沿雷达视线向形变相位，δφ_atm为大气波动引起的相位变化，δφ_noise为噪声相位；

根据已知参数模拟平地相位与地形相位：

${\mathrm{\&delta;\&phi;}}_{flat}=\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}r}\&CenterDot;\frac{b\mathrm{\&Delta;}r}{\tan \mathrm{\&theta;}}---\left(2\right)$

${\mathrm{\&delta;\&phi;}}_{topo}=\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}r}\&CenterDot;\frac{b\mathrm{\&Delta;}h}{\sin \mathrm{\&theta;}}---\left(3\right)$

其中：λ为雷达波长，r为雷达到目标之间的斜距，b为垂直基线长，Δr和Δh分别为像元间的斜距和高程增量；

去除平地、地形相位得到差分干涉相位：

φ_diff＝δφ_mov+δφ_toperror+δφ_atm+δφ_noise (4)

式(4)中δφ_toperror为差分干涉纹图中由于外部DEM误差带来的相位；

步骤三：轨道残余相位去除

地形相位与平地相位的计算均与基线长度相关，而轨道的精度直接影响干涉相对空间基线的计算；由于轨道的不准确性形成的残余相位，即Orbit ramp或Orbit residual；应用精密轨道参数修正初始轨道参数能有效降低这种因素的影响，但当干涉图范围扩大时，二维平面上轨道残余相位的积累将使得形变图上附加一种趋势，这种趋势随着数据覆盖范围的增大而变化；基于最小二乘进行趋势面拟合进而去除轨道残余相位，首先依据先验知识挑选出差分干涉纹图集中有明显残余相位的干涉纹图，针对有误差的干涉纹图首先进行相位解缠，通过提取解缠相位图中一定采样间隔的解缠相位，利用二次或者高次函数进行相位回归，进而拟合趋势面；这里选择二次函数式(5)作为拟合函数进行相位回归；

式中：为采样点解缠相位；(x,y)为采样点图像坐标；a_i为多项式系数；基于最小二乘原则，计算出a_i，从差分相位φ_diff中减去拟合的轨道残余相位；

步骤四：CTInSAR时序分析

当SLC数据集中数据量较少时即m<20，依据相干目标频谱变化小这一特性进行相干目标的选择；点目标的实际几何尺寸往往小于雷达分辨单元，与分布目标具有不同的信号特征，不具备斑点噪声的特性，相位稳定，因而将其视为相干目标；由于点目标稳定的后向散射特性，在获取原始信号过程中各个子视图像上表现出相似的后向散射特征，不同距离向和方位向的子带宽内的散射强度大致相同，因而将SLC图像作反变换，生成复数个子视图，将子视图作谱相关分析，根据其相关性的大小来识别点目标；设各子视图内的某像元的信号x_i(n)＝a_ig(n)，1≤i≤L为视数，g(n)为该像元信号，按照能量比标准衡量L个子视信号的整体相似性为：

$E=\frac{{\left(\underset{i=1}{\overset{L}{\mathrm{\&Sigma;}}}{a}_i\right)}^2}{L\&CenterDot;\underset{i=1}{\overset{L}{\mathrm{\&Sigma;}}}{a_i}^2},0\&le;E\&le;1---\left(6\right)$

当a₁＝a₂＝…＝a_L时，E＝1；能量比越大，整体相似性越好，每幅SLC中频谱特征统计值E高于给定阈值0.4的像元都被识别为点目标候选点，在候选点中将雷达目标的散射强度值高于平均强度值10倍的点作为相干点；

利用相干点将去除了轨道残余相位的差分干涉纹图集的相位进行提取，得到差分干涉相位的点集；相位组成如式(4)所示，其中形变相位δφ_mov包含线性形变δφ_linear和非线性形变δφ_nonlinear两部分；差分干涉相位表达为

$\begin{array}{c}{\mathrm{\&phi;}}_{diff}={\mathrm{\&delta;\&phi;}}_{linear}+{\mathrm{\&delta;\&phi;}}_{toperror}+{\mathrm{\&delta;\&phi;}}_{nonlinear}+{\mathrm{\&delta;\&phi;}}_{atm}+{\mathrm{\&delta;\&phi;}}_{noise}\\ =\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\&CenterDot;v\&CenterDot;T+\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}r}\&CenterDot;\frac{b\mathrm{\&epsiv;}}{\sin \mathrm{\&theta;}}+{\mathrm{\&delta;\&phi;}}_{nonlinear}+{\mathrm{\&delta;\&phi;}}_{atm}+{\mathrm{\&delta;\&phi;}}_{noise}\end{array}---\left(7\right)$

式(7)中，λ为雷达波长，v为线性形变速率，T为生成差分干涉纹图的图像对时间基线，r为雷达到目标之间的斜距，b为垂直基线长，ε为高程误差，θ为雷达入射角；

由于在整个干涉纹图集中每个相干点的线性速率和DEM误差均为常量，是模型中的线性分量，利用以下二维模型能实现形变速率和高程改正的估计：

$\begin{array}{c}{\mathrm{\&delta;\&phi;}}_{\mathrm{mod}el}\left(x_m,y_m,x_n,y_n,T_i\right)=\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\&CenterDot;T_i\&CenterDot;\&lsqb;v_{\mathrm{mod}el}\left(x_m,y_m\right)-v_{\mathrm{mod}el}\left(x_n,y_n\right)\&rsqb;\\ +\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\&CenterDot;\frac{b\left(T_i\right)}{r\left(T_i\right)\&CenterDot;\sin \left({\mathrm{\&theta;}}_i\right)}\&lsqb;{\mathrm{\&epsiv;}}_{\mathrm{mod}el}\left(x_m,y_m\right)-{\mathrm{\&epsiv;}}_{\mathrm{mod}el}\left(x_n,y_n\right)\&rsqb;\end{array}---\left(8\right)$

T_i为第i个干涉像对的时间基线；对于相干目标，由于差分干涉纹图中存在不同的相位补偿，单个像元的相位值在处理过程中并没有实际意义；在确定相位解缠参考点后连接相邻两个相干目标(x_m,y_m)和(x_n,y_n)，通过式(9)求时序模型函数的γ_model最大值得到式(8)中v和ε的估算结果v_estimated和ε_estimated；

${\mathrm{\&gamma;}}_{\mathrm{mod}el}(x_m,y_m,x_n,y_n)=\frac{1}{N}\&CenterDot;|\underset{i=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\exp \&lsqb;j\&CenterDot;\left({\mathrm{\&delta;\&phi;}}_{diff}\right(x_m,y_m,x_n,y_n,T_i)-{\mathrm{\&delta;\&phi;}}_{\mathrm{mod}el}(x_m,y_m,x_n,y_n,T_i\left)\right)\&rsqb;|---\left(9\right)$

式(9)中N为干涉纹图集中包含的干涉纹图个数；

去除线性分量，模型残余相位用下式表示：

δφ_residual(x,y,T_i)＝δφ_nonliear(x,y,T_i)+δφ_atm(x,y,T_i)+δφ_noise(x,y,T_i)

＝δφ_diff(x,y,T_i)-δφ_{linearestmodel}(x,y,T_i)

$={\mathrm{\&delta;\&phi;}}_{diff}(x,y,T_i)-\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\&CenterDot;T_i\&CenterDot;v_{estimated}-\frac{4\mathrm{\&pi;}}{\mathrm{\&lambda;}}\frac{b}{r\&CenterDot;sin\mathrm{\&theta;}}\&CenterDot;{\mathrm{\&epsiv;}}_{estimated}---\left(10\right)$

式中：δφ_diff为初始差分相位，δφ_{linearestmodel}为线性模型将线性形变和DEM误差合成得到的相位；模型残余相位中包含大气相位影响δφ_atmos和非线性移动δφ_nonliear(x,y,T_i)以及失相干因素造成的噪声δφ_noise(x,y,T_i)；

利用1km×1km二维移动平均窗口进行空间低通滤波，滤波完成后，残余相位的噪声δφ_noise(x,y,T_i)被滤除，最后，利用时域高通滤波估计大气贡献δφ_atm，区分出非线性形变相位δφ_nonliear(x,y,T_i)；

利用前面估计出的线性分量与非线性分量获得完整的形变演变量；

$\mathrm{\&rho;}\left(t_i\right)=v_{estimated}\&times;(t_i-t_0)+\frac{\mathrm{\&lambda;}}{4\mathrm{\&pi;}}\&times;{\mathrm{\&phi;}}_{nonlinear}\left(t_i\right)---\left(11\right)$

依据上述相干目标时序分析方法，通过相干目标时序分析模型迭代补偿DEM误差、通过滤波减小大气相位以及噪声的影响，获取到各相干目标点的地面沉降量；

步骤五：长条带相邻分块检验与校正

对长条带中的每个frame重复步骤二至步骤四，得到条带下各分块地面沉降信息；长条带的每个frame获取的形变信息图的参考基准位置不同，需要对连续frame的参考基准进行整体偏差补偿；整体偏差估算依赖于重叠部分相干目标形变参数的统计分析，依据同名点按照式(12)统计连续frame之间形变速率的整体偏差，参与统计的相干目标均为满足相关模型的点；在此过程中选择连续frame中的一个作为基准，其余均相对于该frame进行基准偏差补偿，在每个frame时序分析过程中，根据时间基线的变化关系，将在不同时间间隔下稳定区域内的相干目标作为参考基准，以此降低frame之间形变参数的整体偏差；以两块连续frame之间的检验、校正，选择其中一个frame作为基准，依据frame之间重叠区域同名点的统计信息进行基准偏差的整体补偿

${\mathrm{\&Delta;\&rho;}}_{off}=\frac{1}{N}\underset{0}{\overset{N-1}{\&Sigma;}}({\mathrm{\&rho;}}_{mi}-{\mathrm{\&rho;}}_{si})---\left(12\right)$

${\widehat{\mathrm{\&rho;}}}_s={\mathrm{\&rho;}}_s+{\mathrm{\&Delta;\&rho;}}_{off}---\left(13\right)$

式(12)和(13)中：Δρ_off为frame之间基准偏差；ρ_mi为基准frame重叠区中相干目标i的形变量；ρ_si为需调整frame重叠区中相干目标i的形变量；N表示两个frame包含的同名点个数；为调整的frame中所有相干目标改正后的形变向量；

对于重叠区域内的各相干点目标形变量ρ_i的取值采用距离加权平均方法求解：

${\mathrm{\&rho;}}_i=P_{mi}\&times;{\mathrm{\&rho;}}_{mi}+P_{si}\&times;{\widehat{\mathrm{\&rho;}}}_{si}---\left(14\right)$

$P_{mi}=\frac{D_{mi}}{D},P_{si}=\frac{D_{si}}{D}---\left(15\right)$

式中P_mi和P_si分别为基准frame和调整frame第i个相干目标的权，ρ_mi和分别为重叠区内第i个相干目标在基准frame和调整frame中的形变量，D为重叠区长度，D_mi和D_si分别为基准frame和调整frame中相干目标距离重叠区内各自影像方位向边缘的距离；

步骤六：多轨道坐标系统统一

通过对不同轨道下的多个长条带InSAR地面沉降信息集成进而获取大覆盖范围的监测成果，首先需要对不同的轨道数据进行坐标系统一，具有相同的空间坐标基准；

首先需要获取研究区的DEM数据，根据外部DEM将雷达图像与监测成果进行正射校正，即地理编码，使得不同轨道下的雷达图像位于相同的地面坐标系下，消除因地形起伏以及成像几何关系引起的畸变因素的影响；高程精度优于10m的DEM数据保证了平坦地区配准精度优于一个像元，依据卫星参数将DEM模拟至SAR坐标系，建立初始坐标变换表，再将DEM模拟的影像与SAR图像进行精密配准改进坐标变换表，通过反向变换，将雷达坐标系下各轨道地面形变信息地理编码至DEM所在的地面坐标系内；

步骤七：多轨道信息检验与参考基准统一

位于相同地面坐标系的不同轨道形变信息集成过程中，受制于雷达成像入射角的影响，不能像长条带连续frame重叠区内的相干目标在相邻frame中基本一致，相邻轨道下相干目标是完全不同的，选择其中一条轨道为主轨道，相邻另一条为辅轨道，首先通过反距离加权插值方法将辅轨道上相干目标分析结果进行空间域插值为面，设辅轨道重叠区内共有n个相干目标点，其插值后平面各点坐标(x,y)的沉降速率值通过下式计算得出：

$\mathrm{\&rho;}(x,y)=\frac{\underset{i=1}{\overset{n}{\&Sigma;}}\frac{{\mathrm{\&rho;}}_i}{d_i^2}}{\underset{i=1}{\overset{n}{\&Sigma;}}\frac{1}{d_i^2}}---\left(16\right)$

式中ρ_i为第i个相干目标点的沉降量，d_i为点(x,y)到第i个相干目标点的距离；再利用主轨道上重叠区内相干点目标位置提取其在辅轨道上插值面上的形变值ρ_s；根据相同相干点目标在主轨道的沉降值ρ_m和在辅轨道的沉降值ρ_s进行检验与校正，检验与校正方法同长条带相邻分块的检验与校正。