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

摘要

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

主权项

一种集成多轨道、长条带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为噪声相位；根据已知参数模拟平地相位与地形相位：$\delta {\phi }_{\mathrm{flat}}=\frac{4\pi }{\mathrm{\lambda r}}·\frac{\mathrm{b\Delta r}}{\tan \theta }---\left(2\right)$$\delta {\phi }_{\mathrm{topo}}=\frac{4\pi }{\mathrm{\lambda r}}·\frac{\mathrm{b\Delta h}}{\sin \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}{\Sigma }}{a}_i\right)}^2}{L·\underset{i=1}{\overset{L}{\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}{\phi }_{\mathrm{diff}}=\delta {\phi }_{\mathrm{linear}}+\delta {\phi }_{\mathrm{toperror}}+\delta {\phi }_{\mathrm{nonlinear}}+\delta {\phi }_{\mathrm{atm}}+\delta {\phi }_{\mathrm{noise}}\\ =\frac{4\pi }{\lambda }·v·T+\frac{4\pi }{\mathrm{\lambda r}}·\frac{\mathrm{bϵ}}{\sin \theta }+\delta {\phi }_{\mathrm{nonlinear}}+\delta {\phi }_{\mathrm{atm}}+\delta {\phi }_{\mathrm{noise}}\end{array}---\left(7\right)$式(7)中，λ为雷达波长，v为线性形变速率，T为生成差分干涉纹图的图像对时间基线，r为雷达到目标之间的斜距，b为垂直基线长，ε为高程误差，θ为雷达入射角；由于在整个干涉纹图集中每个相干点的线性速率和DEM误差均为常量，是模型中的线性分量，利用以下二维模型能实现形变速率和高程改正的估计：$\begin{array}{c}\delta {\phi }_{\mathrm{mode}l}(x_m,y_m,x_n,y_n,T_i)=\frac{4\pi }{\lambda }·T_i·[v_{\mathrm{mode}l}(x_m,y_m)-v_{\mathrm{mode}l}(x_n,y_n)]\\ +\frac{4\pi }{\lambda }·\frac{b\left(T_i\right)}{r\left(T_i\right)·\sin \left({\theta }_i\right)}[{ϵ}_{\mathrm{mode}l}(x_m,y_m)-{ϵ}_{\mathrm{mode}l}(x_n,y_n)]\end{array}---\left(8\right)$T_i为第i个干涉像对的时间基线；对于相干目标，由于差分干涉纹图中存在不同的相位补偿，单个像元的相位值在处理过程中并没有实际意义；在确定相位解缠参考点后连接相邻两个相干目标(x_m,y_m)和(x_n,y_n)，通过式(9)求时序模型函数的γ_model最大值得到式(8)中v和ε的估算结果v_estimated和ε_estimated；${\gamma }_{\mathrm{mode}l}(x_m,y_m,x_n,y_n)=\frac{1}{N}·\left|\underset{i=1}{\overset{N-1}{\Sigma }}\exp \right[j·(\delta {\phi }_{\mathrm{diff}}(x_m,y_m,x_n,y_n,T_i)-\delta {\phi }_{\mathrm{mode}l}(x_m,y_m,x_n,y_n,T_i))\left]\right|---\left(9\right)$式(9)中N为干涉纹图集中包含的干涉纹图个数；去除线性分量，模型残余相位用下式表示：$=\delta {\phi }_{\mathrm{diff}}(x,y,T_i)-\frac{4\pi }{\lambda }·T_i·v_{\mathrm{estimated}}-\frac{4\pi }{\lambda }\frac{b}{r·\sin \theta }·{ϵ}_{\mathrm{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)；利用前面估计出的线性分量与非线性分量获得完整的形变演变量；$\rho \left(t_i\right)=v_{\mathrm{estimated}}\times (t_i-t_0)+\frac{\lambda }{4\pi }\times {\phi }_{\mathrm{nonlinear}}\left(t_i\right)---\left(11\right)$依据上述相干目标时序分析方法，通过相干目标时序分析模型迭代补偿DEM误差、通过滤波减小大气相位以及噪声的影响，获取到各相干目标点的地面沉降量；步骤五：长条带相邻分块检验与校正对长条带中的每个frame重复步骤二至步骤四，得到条带下各分块地面沉降信息；长条带的每个frame获取的形变信息图的参考基准位置不同，需要对连续frame的参考基准进行整体偏差补偿；整体偏差估算依赖于重叠部分相干目标形变参数的统计分析，依据同名点按照式(12)统计连续frame之间形变速率的整体偏差，参与统计的相干目标均为满足相关模型的点；在此过程中选择连续frame中的一个作为基准，其余均相对于该frame进行基准偏差补偿，一般情况下，在每个frame时序分析过程中，根据时间基线的变化关系，将在不同时间间隔下稳定区域内的相干目标作为参考基准，以此降低frame之间形变参数的整体偏差；以两块连续frame之间的检验、校正，选择其中一个frame作为基准，依据frame之间重叠区域同名点的统计信息进行基准偏差的整体补偿$\Delta {\rho }_{\mathrm{off}}=\frac{1}{N}\underset{0}{\overset{N-1}{\Sigma }}({\rho }_{\mathrm{mi}}-{\rho }_{\mathrm{si}})---\left(12\right)$${\hat{\rho }}_s={\rho }_s+\Delta {\rho }_{\mathrm{off}}---\left(13\right)$式(12)和(13)中：Δρ_off为frame之间基准偏差；ρ_mi为基准frame重叠区中相干目标i的形变量；ρ_si为需调整frame重叠区中相干目标i的形变量；N表示两个frame包含的同名点个数；为调整的frame中所有相干目标改正后的形变向量；对于重叠区域内的各相干点目标形变量ρ_i的取值采用距离加权平均方法求解：${\rho }_i=P_{\mathrm{mi}}\times {\rho }_{\mathrm{mi}}+P_{\mathrm{si}}\times {\hat{\rho }}_{\mathrm{si}}---\left(14\right)$$P_{\mathrm{mi}}=\frac{D_{\mathrm{mi}}}{D},P_{\mathrm{si}}=\frac{D_{\mathrm{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)的沉降速率值通过下式计算得出：$\rho (x,y)=\frac{\underset{i=1}{\overset{n}{\Sigma }}\frac{{\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进行检验与校正，检验与校正方法同长条带相邻分块的检验与校正。