CN104111457B - 一种升降轨PSInSAR地面沉降监测结果的互检验与时序融合方法 - Google Patents

Abstract

一种升降轨PSInSAR地面沉降监测结果的互检验与时序融合方法，它有五大步骤：一、数据选取及对升降轨雷达数据分别进行PSInSAR处理；二、升降轨PSInSAR观测值坐标系的统一；三、升降轨PSInSAR观测值参考基准的统一；四、基准偏差补偿后升降轨PSInSAR观测沉降速率的相关性计算及精度互检验；五、升降轨PSInSAR观测形变即沉降累积量序列的时序融合。本发明通过计算主辅轨道观测序列的整体偏差和因起始时间差引起的形变累积量差，实现主辅轨道PSInSAR观测形变累积量数据的时序融合，加密了研究区PSInSAR的观测形变序列，从而准确、精细地恢复出研究区地面沉降的非线性动态变化过程。

Description

一种升降轨PSInSAR地面沉降监测结果的互检验与时序融合方法

技术领域

本发明涉及一种升降轨PSInSAR地面沉降监测结果的互检验与时序融合方法，属于合成孔径雷达干涉技术领域。它对升降轨雷达数据分别进行PSInSAR数据处理后，将获取的2组视线向形变量转换为沉降量进行沉降监测精度互检验，在检验InSAR观测的内符合精度的基础上，完成对所有PSInSAR观测的相干目标累积沉降量的时序融合，从而实现加密PSInSAR监测地面沉降形变序列的目的，这一方法可以完整刻画出研究区地面沉降的非线性变化过程。

背景技术

PSInSAR，即永久散射体干涉测量，是一种基于雷达图像的相位数据来提取地面目标三维空间信息的重要技术。经过近10年的发展，在华北平原、长江三角洲等地区已规模化应用，这一技术已成为我国平原区地面沉降监测主要的空间技术手段。PSInSAR技术的核心是基于大量的SAR数据(一般大于20甚至30景)，对永久散射体(即PS点)的干涉相位进行时间序列分析，最终获取每个PS点的线性、非线性形变速率以及形变累积量序列。这种方法一般对雷达影像个数有严格要求，一方面，只有雷达影像个数达到一定数量才能筛选出在整个时间跨度内都具有稳定信号的PS点；另一方面，为了提高PS点的估算精度，也需要大量的雷达数据进行统计分析。于是，通常结合研究区地面沉降速率的大小以及目前在轨的雷达卫星的波长，选取重复周期合适的雷达卫星，订购该卫星的雷达数据并积累到PSInSAR技术对数据量的要求，即可真实地获得研究区的地面沉降状况。如果研究区为城区，可获取的同一模式的雷达数据量较少，但由于该区存在较多相干目标且在长时间间隔下保持了很高的相干性，此时进行PSInSAR数据处理也可得到较高的监测精度。

但是，在地面沉降非线性特征较为明显的地区，由于雷达卫星的重复周期较长，即在时间域同一模式雷达数据的采样密度较小，或者在相同时间跨度内同一模式雷达数据量较少，当用PSInSAR技术获取该区的地面沉降状况时，研究区地面沉降的非线性变化特征并不能完整刻画。相同模式雷达数据中，同一相干目标点p在相邻两时间点mT、nT上的差分相位可由式(1)表达：

上式中，λ是雷达载波波长，s(t)是地面沉降在雷达视线向的分量，m＝0,1,2,…,N，是获 取的同一模式雷达影像的景数，n＝0,1,2,…,N-1,且n＜m，T是雷达卫星的重访周期，是一固定值。

在雷达图像获取的相邻时间内，即在mT时刻与nT时刻之间，同一相干目标相位差的绝对值 都小于或等于π，满足PSInSAR技术相位解缠的条件，从理论上说，通过相位解缠可以获得研究区正确的地面沉降结果。但是，研究区地面沉降的非线性变化特征非常明显，地面沉降的动态变化过程难以准确恢复。由于雷达卫星的重复周期固定，同一模式雷达数据的观测采样密度不变，单纯增大同一模式的雷达数据量也无法解决这一问题。但是，通过尽可能多地获取相同时间跨度内研究区升降轨两种模式的雷达数据，在将这两种模式的PSInSAR监测结果统一到同一坐标系和同一参考基准的前提下，研究区PSInSAR监测结果的密度得到加强，成功恢复研究区地面沉降非线性动态变化过程的困难将会得到缓解。

从理论上说，PSInSAR测量地表形变的精度可以达到亚毫米级，但在数据处理过程中，所用雷达图像的个数、PS点与参考点之间的距离、相位解缠的可靠性以及大气影响的大小等都会降低PSInSAR的测量精度，因此，客观地评价PSInSAR测量结果的精度，衡量PSInSAR技术的可靠性是首要任务。一般地，采用精密水准(一、二等水准测量)或GPS观测对PSInSAR技术的精度进行评定，但这些测量点位密度远远低于相干目标点的密度，其与相干目标点的空间位置存在偏差，且它们的工作周期较长，往往需要长达1-2个月，其观测值的时效性将大大降低。相对于地面测量数据(水准、GPS测量等)，升降轨不同模式下获取的PSInSAR监测结果克服了地面测量数据的不足，且在观测时间同步的条件下，由于同一地区升降轨模式下分别获取的PSInSAR监测结果是独立的，也降低了PSInSAR函数模型自身因素的影响，为InSAR观测精度的内部检验提供了基础。

本专利申请针对同一研究区的地面沉降场，获取升轨和降轨两种模式的雷达数据，然后分别进行PSInSAR数据处理并获得相应的沉降速率，以升降轨PSInSAR观测的沉降速率值为比较对象，在确定主辅轨道的基础上对辅轨道观测数据进行基准补偿。在主辅轨道(升降轨)坐标系和参考基准统一的前提下，以主辅轨道观测值互差的均方差为检验指标进行PSInSAR观测值的统计检验和精度评价。最后，在此基础上对辅轨道的形变观测序列进行沉降速率整体偏差和因起始时间差引起的形变累积量差补偿，从而实现主辅轨道观测序列数据的时序融合，融合后的监测数据加密了PSInSAR监测地面沉降的观测序列，有效地凸显出研究区的非线性动态变化特征。

发明内容

1.目的：本发明的目的是提供一种升降轨PSInSAR地面沉降监测结果的互检验与时序融合方法，它克服了研究区同一模式雷达数据量不足的缺陷，通过获取该区升降轨模式的雷达数据，不仅实现了不同轨道InSAR观测的内符合精度检验，还在主辅轨道(升降轨)坐标系和参考基准统一的前提下，通过计算主辅轨道(升降轨)观测序列的整体偏差和因起始时间差引起的形变累积量差，实现主辅轨道(升降轨)PSInSAR观测形变(沉降)累积量数据的时序融合，加密了研究区PSInSAR的观测形变序列，从而准确、精细地恢复出研究区地面沉降的非线性动态变化过程。

2.技术方案：本发明是一种升降轨PSInSAR地面沉降监测结果的互检验与时序融合方法，该方法具体步骤如下：

步骤一：数据选取及对升降轨雷达数据分别进行PSInSAR处理

基于研究区域地面沉降速率的大小，选取雷达波长及重复周期都较为合适的雷达卫星的升轨数据和降轨数据。

针对雷达卫星的升轨数据和降轨数据，要分别进行PSInSAR处理，分别获取升降轨模式下PSInSAR观测的地面沉降速率和形变累积量序列。首先，利用多准则相干目标识别算法，最大化地提取每一轨道下的相干目标，以Delaunay三角网连接该轨道下所有相干目标，根据相邻相干目标的差分干涉相位差序列，利用二维周期图估算相干点间的形变速率和高程误差改正。以此为基础，进行该轨道下残余相位的滤波处理，得到不同时刻相干点的非线性形变和大气估计，最终获取该轨道下每个相干目标的地面沉降速率值和形变序列。

步骤二：升降轨PSInSAR观测值坐标系的统一

为对比升降轨不同模式下PSInSAR观测的结果，必须首先统一它们的坐标系统，使得升降轨模式下雷达数据获取的地表形变观测值(地面沉降速率和形变序列)处于相同的参考坐标系中，具有相同的空间基准。

坐标系统一可通过2种方式来实现，即雷达坐标系下的统一和地面坐标下的统一。雷达坐标系下的统一是指将具有不同角度的SAR图像直接进行精确配准。由于雷达波入射方向和入射角的差异而造成SAR数据不均匀变形，使得这一配准方式的精度受限；地面坐标下的统一则是要先进行正射校正，即地理编码，使得在升降轨下获取的雷达图像位于相同的地面坐标系下，以消除因地形起伏引起的畸变影响，然后进行升降轨雷达图像的精确匹配。这种方式需要获取研究区的DEM数据，DEM数据的精度决定了地理编码的精度。高程精度优于10m的DEM数据保证了平坦地区配准精度优于1个像元。

本发明专利申请采用地面坐标系统一的方式进行升降轨PSInSAR观测结果的坐标转换。将升降轨雷达影像分别进行地理编码，并实现升降轨PSInSAR处理获取的沉降速率图的坐标转换。确定主轨道，在获取位于地面坐标系下主轨道的平均雷达强度图像后，应用多项式纠正模型，即式(2)，完成辅轨道与主轨道间在地面坐标系下的精确配准，以进一步提高主辅轨道间沉降速率图的配准精度。

$\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---\left(2\right)$

式中：Δx、Δy分别为主辅影像x、y方向的相对偏移量；p为多项式阶数，选择3阶多项式；a、b为多项式系数。

步骤三：升降轨PSInSAR观测值参考基准的统一

受制于雷达波入射角和入射方向的影响，同一相干目标分别在升降轨雷达图像中对应的位置略有差别。因而，升降轨PSInSAR观测结果在平面上是“浮动的”，它们之间的绝对变化量取决于参考基准的变化，实现升降轨PSInSAR观测结果参考基准统一的本质是解决不同参考位置间的整体偏差。

进行参考位置的统一，需以主轨道为基准。首先，在选定主轨道的前提下，提取主轨道上相干目标的地面沉降速率；其次，对辅轨道的PSInSAR观测数据进行插值处理，生成连续分布的沉降面；然后，根据主轨道上相干目标的空间位置，提取辅轨道上相应位置的沉降参数；最后，利用统计方法，即式(3)(4)解算主辅轨道之间的整体偏差，进而完成辅轨道与主轨道的统一。

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

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

式中：Δv_off为主辅轨道PSInSAR监测沉降速率间的整体基准偏差；和分别为主辅轨道下相干目标i的形变(沉降)速率；N为统计样本数；为基准偏差补偿后辅影像中相干目标的形变(沉降)速率。

步骤四：基准偏差补偿后升降轨PSInSAR观测沉降速率的相关性计算及精度互检验

在完成主辅轨道(升降轨)PSInSAR观测值坐标系和参考基准统一的基础上，对主辅轨道(升降轨)PSInSAR观测样本进行相关性计算，进一步证明主辅轨道(升降轨)间PSInSAR观测结果的不同是由主辅轨道(升降轨)不同参考位置间的整体偏差造成的。

为检验主辅轨道(升降轨)上PSInSAR观测形变(沉降)速率的精度，直接利用提取的主辅轨道(升降轨)上相同相干目标的PSInSAR观测速率值进行统计比较，以互差的均方差为统计指标，按照式(5)计算：

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

式中：X_i为主轨道第i个相干目标的PSInSAR观测值，Y_i为辅轨道第i个相干目标的经过基准偏差补偿后的PSInSAR观测值，P为统计样本点个数。

步骤五：升降轨PSInSAR观测形变(沉降)累积量序列的时序融合

假设主轨道上M幅、辅轨道上N幅PSInSAR形变累积量图对应的序列值分别为

$D^m=[d_0^m,d_1^m,...,d_{M-1}^m]---\left(6\right)$

$D^s=[d_0^s,d_1^s,...,d_{N-1}^s]---\left(7\right)$

其中，D^m为主轨道上M幅PSInSAR形变观测图中相干目标对应的形变累积量序列，D^s为辅轨道上N幅PSInSAR形变观测图中相干目标对应的形变累积量序列，为主轨道上第i幅PSInSAR形变观测图中相干目标的形变累积量(i＝0,1,...,M-1)，为辅轨道上第j幅PSInSAR形变观测图中相干目标的形变累积量(j＝0,1,...,N-1)，m为主轨道，s为辅轨道。

显然，这2组PSInSAR观测序列存在参考基准整体速率偏差和因起始时间差引起的形变累积量差，需要逐个进行补偿。对辅轨道下的形变序列进行整体速率偏差修正，以求得参考基准补偿后辅轨道下第j幅形变观测图中相干目标的形变(沉降)累积量可表达为

$d_j^{s^{\′}}={\mathrm{\Δv}}_{\mathrm{off}}\·t_j+d_j^s(j=\mathrm{0,1},...,N-1)---\left(8\right)$

式中，为辅轨道下j幅形变观测图中相干目标经整体速率偏差补偿后的形变累积量，为辅轨道下j幅形变观测图中相干目标沉降速率偏差补偿前的形变累积量，Δv_off为主辅轨道间整体速率观测值基准偏差，t_j为辅轨道下生成第j幅形变观测图的两幅雷达图像获取时刻的时间间隔(单位为年)。

在完成速率偏差修正后，需将这2组PSInSAR观测值统一到相同的时间起点。此时，形变序列的补偿量为辅轨道相对于主轨道观测序列起始时刻的形变量偏差，将之加到序列累积量上即可得到辅轨道下第j幅形变观测图中相干目标经整体速率偏差和因起始时间差引起的 形变累积量差补偿后的形变(沉降)累积量，即

$d_j^{s^{\′\′}}=v_m\·\mathrm{\ΔT}+d_j^{s^{\′}}---\left(9\right)$

式中：为辅轨道下第j幅形变观测图中相干目标经整体速率偏差和因起始时间差引起的形变累积量差补偿后的形变(沉降)累积量，为辅轨道下下j幅形变观测图中相干目标经整体速率偏差补偿后的形变累积量，v_m为同名相干目标在主轨道下的沉降速率，ΔT为辅轨道相对主轨道在起始时刻的时间差(单位为年)。

在完成相干目标的整体速率偏差补偿和时间偏差补偿后，主辅轨道下的形变累积量序列统一到相同的时间和空间参考基准上，从而得到相干目标加密后的PSInSAR形变(沉降)观测序列为

D＝[d₀,d₁,…d_M-1,…,d_M+N-1] (10)

式中，D为相干目标加密后PSInSAR形变(沉降)累积量观测序列，d_i(i＝1,2,3,...,M+N-1)为主轨道相干目标或辅轨道经基准补偿后相干目标的PSInSAR形变累积量。

3.优点及功效：本发明提供了一种升降轨PSInSAR地面沉降监测结果的互检验和时序融合方法，其优点为：

(1)本方法充分利用同一研究区的升降轨雷达数据，在分别获取升降轨模式下PSInSAR观测速率值的基础上，进行主辅轨道间观测形变(沉降)序列的基准偏差补偿，从而实现主辅轨道观测序列的时序融合，有效地加密了对同一研究区PSInSAR监测的序列。

(2)利用升降轨PSInSAR观测的时序融合方法，可以获取同一研究区更多的PSInSAR形变观测序列，进而更加精准、有效地恢复研究区地面沉降的非线性动态变化过程。

(3)本方法提供了在参考基准统一下升降轨PSInSAR地面沉降速率值的比较，实现了PSInSAR观测结果的内部检验。结果表明：在雷达数据量为20～30景的条件下，地面沉降速率监测结果的互检验精度优于2mm，进一步印证了InSAR技术监测地表形变结果的准确性。

(4)本方法不仅解决了升降轨雷达数据监测同一研究区的地面沉降问题，有利于进行地表二维甚至三维活动量的提取，还为不同轨道间地表形变监测结果集成提供了理论依据。

附图说明

图1.本发明的流程图。

图2.升轨Track-39、降轨Track-275雷达数据在研究区的分布图，其中，底图中31、32代表北纬纬度数；底图中119、120、121、122代表东经经度数。

图3(a)主轨道(降轨Track-275)下，PSInSAR监测的地面坐标系下的地面沉降速率图。

图3(b)辅轨道(升轨Track-39)下，PSInSAR监测的地面坐标系下的地面沉降速率图。

图4(a)主轨道(降轨Track-275)相干目标样本点地面沉降速率统计直方图。

图4(b)辅轨道(升轨Track-39)相干目标样本点地面沉降速率统计直方图。

图5.辅轨道(升轨Track-39)相对于主轨道(降轨Track-275)基准补偿后PSInSAR监测地面沉降速率图，P1、P2分别为基准补偿后辅轨道中相干目标的两个样本点。

图6(a)为辅轨道(升轨Track-39)相对于主轨道(降轨Track-275)基准补偿前的相关统计图。

图6(b)为辅轨道(升轨Track-39)相对于主轨道(降轨Track-275)基准补偿后的相关统计图。

图7(a)为图5中相干目标样本点P1加密后的PSInSAR监测沉降累积量观测序列。

图7(b)为图5中相干目标样本点P2加密后的PSInSAR监测沉降累积量观测序列。

具体实施方式

以覆盖苏州市西北地区的ENVISAT升降轨雷达数据为例，说明本发明在实际工程应用中的具体操作步骤。见图1，本发明是一种升降轨PSInSAR地面沉降监测结果的互检验与时序融合方法，该方法具体步骤如下：

步骤一：数据选取及对升降轨雷达数据分别进行PSInSAR处理

1.数据选取

选取多幅覆盖工作区的升降轨星载合成孔径雷达SLC图像，对工作区的地面沉降进行监测。为此目的，选取欧空局ENVISAT星载C波段20m分辨率的时间跨度为2006年09月至2010年08月、两条轨道分别为升轨Track-39的24幅和降轨Track-275的27幅SLC数据，见图2。数据幅宽100公里，长100公里，数据覆盖时间跨距四年，卫星重复过境周期为35天。雷达主要数据参数和扫描日期见表1和表2所示：

表1：升轨Track-39轨道雷达数据日期表

20071206	20080214	20080320	20080424	20080529	20080807
						20080911	20081016	20081225	20090129	20090305	20090409
20090514	20090618	20090723	20090827	20091001	20091210
						20100218	20100325	20100429	20100603	20100708	20100812

表2：降轨Track-275轨道雷达数据日期表

20060924	20061029	20061203	20070318	20070422	20070527	20070701
							20070805	20070909	20071014	20071223	20080127	20080511	20080615
20080720	20081207	20090111	20090215	20090322	20090426	20090809
							20100131	20100307	20100411	20100620	20100725	20100829

2.对升降轨雷达数据分别进行PSInSAR处理

首先处理升轨Track-39的24幅雷达图像，选取其中一幅图像为主图像，原则上它应该位于时域(时间基线集)和空间域(空间基线集)的中心，将该轨道的其它图像与主图像进行配准后形成干涉像对，基于一定准则或阈值，选择具有高相干的像元作为点目标(即PS点)，进而生成PS点的干涉纹图及差分干涉纹图。任一PS点的差分干涉相位可由式(11)表达

$\begin{array}{c}{\mathrm{\φ}}_{i,\mathrm{diff}}={\mathrm{\φ}}_{i,\mathrm{topo}-\mathrm{error}}+{\mathrm{\φ}}_{i,\mathrm{linear}}+{\mathrm{\φ}}_{i,\mathrm{non}-\mathrm{linear}}+{\mathrm{\φ}}_{i,\mathrm{atm}}+{\mathrm{\φ}}_{i,\mathrm{nosie}}\\ =\frac{4\mathrm{\π}}{\mathrm{\λ}R}\frac{B_{\⊥}}{\sin \mathrm{\θ}}\·{\mathrm{\ϵ}}_i+\frac{4\mathrm{\π}}{\mathrm{\λ}}{v}_i\·T+{\mathrm{\φ}}_{i,\mathrm{non}-\mathrm{linear}}+{\mathrm{\φ}}_{i,\mathrm{atm}}+{\mathrm{\φ}}_{i,\mathrm{nosie}}\end{array}---\left(11\right)$

式(11)中，φ_i,diff为第i个PS点的差分干涉相位，φ_{i,topo -error}为由于外部DEM数据不准确在第i个PS点上产生的地形误差相位，φ_i,linear为地表线性变化在第i个PS点上引起的形变相位，φ_{i,non -linear}为地表非线性变化在第i个PS点上引起的形变相位，φ_i,latm为大气波动在第i个PS点上引起的误差相位，φ_i,noise为随机噪声在第i个PS点上引起的误差相位，可忽略不计，R为雷达数据覆盖区域中心点与主轨道卫星天线之间的距离，θ为雷达卫星的入射角，λ为雷达波的波长，B_⊥为两景雷达数据获取时刻雷达卫星之间的垂直基线长度，T为时间基线，即获取两景雷达数据的时间差(单位：秒)，ε_i为第i个PS点上的高程误差值，v_i为第i个PS点上 的地表线性形变速率。

为了避免差分干涉相位整周模糊度问题，对相邻两点P_o与P_i间的相位差模型，即式(12)，进行分析，用二维周期图估算相干点间的线性形变速率Δv_io和相对高程误差Δε_io。在忽略噪声相位Δφ_io,nosie的情况下，基于非线性形变在时域呈低频、空间域呈高频变化以及大气波动在时域呈高频、空间域呈相对低频变化，对残余相位Δφ_io,residual分别进行时域和空间域的低频滤波处理，得到PS点间的非线性形变相位Δφ_{io,non-linear}和大气估计相位Δφ_io,atm，最后，通过大气相位的去除和线性、非线性形变的叠加，获取每个PS点处的地面沉降速率和地面沉降累积量序列。

式(12)中，为定值，λ，R，θ，B_⊥见式(10)，为相邻点P_o与P_i间的差分干涉相位差，Δφ_{io,topo -error}、Δε_io分别为相邻点P_o与P_i间的地形误差相位差、地形误差之差，Δφ_io,linear、Δv_io分别为相邻点P_o与P_i间的线性形变相位差、线性速率之差，Δφ_{io,non -linear}、Δφ_io,atm、Δφ_io,noise分别为相邻点P_o与P_i间的非线性形变相位差、大气波动相位差以及噪声相位差，Δφ_io,residual＝Δφ_{io,non -linear}+Δφ_io,atm+Δφ_io,noise为相邻点P_o与P_i间的残差相位之差。

同样地，对降轨Track-275的27幅雷达数据也进行PSInSAR处理，获取轨道下每个PS点的地面沉降速率和形变(沉降)观测序列。

步骤二：升降轨PSInSAR观测值坐标系的统一

根据升轨Track-39、降轨Track-275中各自生成的精准查找表，将升降轨下获得的PSInSAR观测序列，即地面沉降速率图、地面沉降累积量图以及升降轨下经配准的SLC图像的平均雷达强度图分别进行反变换，将其坐标系统从雷达坐标系(距离多普勒坐标系，Range-Doppler Coordinate)转换到地面坐标系(EQA，(Equiangular coordinates(lat/long))。图3(a)、(b)分别为主轨道(降轨Track-275)、辅轨道(升轨Track-39)在地面坐标系下PSInSAR监测的地面沉降速率图。

选取降轨轨道Track-275为主轨道，以Track-275轨道在地面坐标系下的平均雷达强度图为主图像，利用多项式纠正模型，即式(2)，将地面坐标系下主辅轨道的地面沉降速率图与之进行精配准，以进一步提高在地面坐标系下PSInSAR观测的地面沉降速率图的精度。

步骤三：升降轨PSInSAR观测值参考基准的统一

受雷达波入射方向和入射角的影响，升降轨图像中相干目标(PS点)的分布密度和位置不尽相同，绝对提取每个相干目标的观测值难以实现，需采用统计方法分别提取同一组相干目标对应的2组PSInSAR观测值。于是，首先利用邻近点插值法，以100m为半径对地面坐标系下辅轨道(升轨Track-39)下的PSInSAR监测的地面沉降速率进行插值处理，生成连续分布的沉降面。以降轨轨道Track-275为主轨道，利用该轨道覆盖范围内的相干目标提取对应的升轨轨道Track-39下的地面沉降速率值，共计提取出升降轨共有53 361对相干目标点的地面沉降速率值。对这2组数据分别进行直方图统计(图4a、b)以及相关性分析(图6(a))可知，这2组观测值的分布线性函数特征非常明显，它们的主要差异为存在整体偏差。根据提取出的53361对样本点，按照式(3)(4)，计算出辅轨道(升轨Track-39)下PSInSAR监测的地面沉降速率相对于主轨道(降轨Track-275)下的整体均值偏差为1.18mm，将辅轨道(升轨Track-39)下PSInSAR监测的地面沉降速率值统一都加上该整体偏差值，继而实现辅轨道(升轨Track-39)与主轨道(降轨Track-275)下PSInSAR监测的地面沉降速率值参考基准的统一。图5即为辅轨道(升轨Track-39)相对于主轨道(降轨Track-275)基准补偿后，PSInSAR监测的地面沉降速率图。

步骤四：基准偏差补偿后升降轨PSInSAR观测沉降速率的相关性计算及精度互检验

将基准补偿后升降轨的地面沉降速率值进行相关性分析可知(图6(b))，升降轨模式下研究区的地面沉降速率整体仍呈线性特征分布，于是，利用线性函数，将基准偏差补偿后辅轨道(升轨Track-39)的沉降速率与主轨道(降轨Track-275)的沉降速率按照式(13)统计计算，

Y＝A+BX (13)

式(13)中，X为辅轨道(升轨Track-39)基准补偿后的PSInSAR监测速率，Y为主轨道(降轨Track-275)下PSInSAR观测速率，A,B为常数。

经过线性函数拟合得到：A＝0.17mm，为整体偏差；B＝1.008为线性函数的斜率，接近于1，表明修正后升轨样本点地面沉降速率与降轨样本点地面沉降速率值呈线性函数分布，且二者之间的整体偏差已经降到很小，它们的参考基准不再有显著差异，可忽略不计。

为检验主辅轨道(升降轨)上PSInSAR观测沉降速率的精度，为进一步检验基准补偿后这两组数据的相关性，以互差的均方差为统计指标，按照式(14)，对提取的53361个相干目标点进行统计计算，可得均方差m仅为±1.824mm。这从另一方面说明主辅轨道PSInSAR观测的地面沉降速率参考基准已基本统一。

$m=\±\sqrt{\frac{\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{(X_I-Y_i)}^2}{P}}---\left(14\right)$

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

步骤五：主辅轨道PSInSAR观测形变(沉降)累积量序列的时序融合

主轨道(降轨Track-275)上27幅、辅轨道(升轨Track-39)上24幅PSInSAR监测的形变(沉降)累积量序列可表示为

$D^m=[d_0^m,d_1^m,...,d_{26}^m]---\left(15\right)$

$D^s=[d_0^s,d_1^s,...,d_{23}^s]---\left(16\right)$

已知辅轨道(升轨Track-39)下PSInSAR观测序列的时间间隔为35天，通过步骤三计算得到的辅轨道(升轨Track-39)下PSInSAR监测的地面沉降速率相对于主轨道(降轨Track-275)下沉降速率的整体均值偏差为1.18mm，根据基准补偿前辅轨道(升轨Track-39)的累积量序列(16)以及式(8)，即可实现对辅轨道(升轨Track-39)下形变累积量序列整体速率偏差的修正，此时得到参考基准补偿后辅轨道(升轨Track-39)下相干目标的形变(沉降)累积量序列为

$D^{s^{\′}}=[d_0^{s^{\′}},d_1^{s^{\′}},...,d_{23}^{s^{\′}}]---\left(17\right)$

在完成整体速率偏差修正后，主辅轨道(升降轨)的PSInSAR观测值还存在因起始时间不同而引起的形变累积量差。由式(9)，根据主轨道(降轨Track-275)相干目标点上的沉降速率、辅轨道(升轨Track-39)相干目标点整体速率偏差修正后的结果以及辅轨道(升轨Track-39)相对主轨道(降轨Track-275)在起始时刻的时间差，即可将主辅轨道的观测序列统一到相同的时间起点，此时得到辅轨道(升轨Track-39)下PSInSAR观测的形变累积量序列由(18)表达：

$D^{s^{\′\′}}=[d_0^{s^{\′\′}},d_1^{s^{\′\′}},...,d_{23}^{s^{\′\′}}]---\left(18\right)$

完成辅轨道(升轨Track-39)相干目标相对于主轨道(降轨Track-275)的整体速率偏差补偿和时间偏差补偿后，主辅轨道下的形变累积量序列即可统一到相同的时间和空间参考基 准上，此时苏州市西北部加密后的PSInSAR形变(沉降)累积量序列为：

$D=[d_0,d_1,...d_{26},d_0^{s^{\′\′}},d_1^{s^{\′\′}},...,d_{23}^{s^{\′\′}}]---\left(19\right)$

图7(a)、(b)所示分别为图5中样本点P1和P2加密后的PSInSAR形变(沉降)累积量序列。

本发明提出的升降轨PSInSAR地面沉降监测结果的互检验方法，摒弃了地面测量数据(水准测量、GPS观测等)时效性低、点密度少等缺陷，是对同一目标的多角度独立观测。它从InSAR技术本身说明，当雷达数据为20～30景时，PSInSAR观测的沉降速率互检验精度优于2mm，进一步表明PSInSAR技术方法在监测以垂向下沉为主的地面沉降的可靠性。同时，在坐标系统一和参考基准统一的前提下，本发明不仅解决了升降轨雷达数据更为精细地观测同一研究区地面沉降动态变化特征的问题，还为大区域地面沉降监测结果的集成提供了理论依据。

Claims (1)

1.一种升降轨PSInSAR地面沉降监测结果的互检验与时序融合方法，该方法具体步骤如下：

步骤一：数据选取及对升降轨雷达数据分别进行PSInSAR处理

基于研究区域地面沉降速率的大小，选取雷达波长及重复周期都较为合适的雷达卫星的升轨数据和降轨数据；

针对雷达卫星的升轨数据和降轨数据，要分别进行PSInSAR处理，分别获取升降轨模式下PSInSAR观测的地面沉降速率和形变累积量序列；首先，利用多准则相干目标识别算法，最大化地提取每一轨道下的相干目标，以Delaunay三角网连接该轨道下所有相干目标，根据相邻相干目标的差分干涉相位差序列，利用二维周期图估算相干点间的形变速率和高程误差改正，以此为基础，进行该轨道下残余相位的滤波处理，得到不同时刻相干点的非线性形变和大气估计，最终获取该轨道下每个相干目标的地面沉降速率值和形变序列；

步骤二：升降轨PSInSAR观测值坐标系的统一

为对比升降轨不同模式下PSInSAR观测的结果，必须首先统一它们的坐标系统，使得升降轨模式下雷达数据获取的地表形变观测值处于相同的参考坐标系中，具有相同的空间基准；

坐标系统一通过2种方式来实现，即雷达坐标系下的统一和地面坐标下的统一；雷达坐标系下的统一是指将具有不同角度的SAR图像直接进行精确配准；由于雷达波入射方向和入射角的差异而造成SAR数据不均匀变形，使得这一配准方式的精度受限；地面坐标下的统一则是要先进行正射校正，即地理编码，使得在升降轨下获取的雷达图像位于相同的地面坐标系下，以消除因地形起伏引起的畸变影响，然后进行升降轨雷达图像的精确匹配；这种方式需要获取研究区的DEM数据，DEM数据的精度决定了地理编码的精度，高程精度优于10m的DEM数据保证了平坦地区配准精度优于1个像元；

采用地面坐标系统一的方式进行升降轨PSInSAR观测结果的坐标转换；将升降轨雷达影像分别进行地理编码，并实现升降轨PSInSAR处理获取的沉降速率图的坐标转换；确定主轨道，在获取位于地面坐标系下主轨道的平均雷达强度图像后，应用多项式纠正模型，即式(2)，完成辅轨道与主轨道间在地面坐标系下的精确配准，以进一步提高主辅轨道间沉降速率图的配准精度；

$\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---\left(2\right)$

式中：Δx、Δy分别为主辅影像x、y方向的相对偏移量；p为多项式阶数，选择3阶多项式；a、b为多项式系数；

步骤三：升降轨PSInSAR观测值参考基准的统一

受制于雷达波入射角和入射方向的影响，同一相干目标分别在升降轨雷达图像中对应的位置略有差别，因而，升降轨PSInSAR观测结果在平面上是“浮动的”，它们之间的绝对变化量取决于参考基准的变化，实现升降轨PSInSAR观测结果参考基准统一的本质是解决不同参考位置间的整体偏差；

进行参考位置的统一，需以主轨道为基准；首先，在选定主轨道的前提下，提取主轨道上相干目标的地面沉降速率；其次，对辅轨道的PSInSAR观测数据进行插值处理，生成连续分布的沉降面；然后，根据主轨道上相干目标的空间位置，提取辅轨道上相应位置的沉降参数；最后，利用统计方法，即式(3)(4)解算主辅轨道之间的整体偏差，进而完成辅轨道与主轨道的统一；

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

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

式中：Δv_off为主辅轨道PSInSAR监测沉降速率间的整体基准偏差；和分别为主辅轨道下相干目标i的形变即沉降速率；N为统计样本数；为基准偏差补偿后辅影像中相干目标的形变即沉降速率；

步骤四：基准偏差补偿后升降轨PSInSAR观测沉降速率的相关性计算及精度互检验

在完成主辅轨道PSInSAR观测值坐标系和参考基准统一的基础上，对主辅轨道PSInSAR观测样本进行相关性计算，进一步证明主辅轨道间PSInSAR观测结果的不同是由主辅轨道不同参考位置间的整体偏差造成的；

为检验主辅轨道上PSInSAR观测形变即沉降速率的精度，直接利用提取的主辅轨道上相同相干目标的PSInSAR观测速率值进行统计比较，以互差的均方差为统计指标，按照式(5)计算：

$m=\±\sqrt{\frac{\underset{i=1}{\overset{P}{\Σ}}{(X_i-Y_i)}^2}{P}}---\left(5\right)$

式中：X_i为主轨道第i个相干目标的PSInSAR观测值，Y_i为辅轨道第i个相干目标的经过基准偏差补偿后的PSInSAR观测值，P为统计样本点个数；

其特征在于：

步骤五：升降轨PSInSAR观测形变即沉降累积量序列的时序融合

设主轨道上M幅、辅轨道上N幅PSInSAR形变累积量图对应的序列值分别为

$D^m=\[d_0^m,d_1^m,...,d_{M-1}^m\]---\left(6\right)$

$D^s=\[d_0^s,d_1^s,...,d_{N-1}^s\]---\left(7\right)$

其中，D^m为主轨道上M幅PSInSAR形变观测图中相干目标对应的形变累积量序列，D^s为辅轨道上N幅PSInSAR形变观测图中相干目标对应的形变累积量序列，为主轨道上第i幅PSInSAR形变观测图中相干目标的形变累积量(i＝0,1,...,M-1)，为辅轨道上第j幅PSInSAR形变观测图中相干目标的形变累积量(j＝0,1,...,N-1)，m为主轨道，s为辅轨道；

显然，这2组PSInSAR观测序列存在参考基准整体速率偏差和因起始时间差引起的形变累积量差，需要逐个进行补偿；对辅轨道下的形变序列进行整体速率偏差修正，以求得参考基准补偿后辅轨道下第j幅形变观测图中相干目标的形变即沉降累积量表达为

$d_j^{s^{\′}}={\mathrm{\Δv}}_{off}\·t_j+d_j^s,(j=0,1,...,N-1)---\left(8\right)$

式中，为辅轨道下j幅形变观测图中相干目标经整体速率偏差补偿后的形变累积量，为辅轨道下j幅形变观测图中相干目标沉降速率偏差补偿前的形变累积量，Δv_off为主辅轨道间整体速率观测值基准偏差，t_j为辅轨道下生成第j幅形变观测图的两幅雷达图像获取时刻的时间间隔，单位为年；

在完成速率偏差修正后，需将这2组PSInSAR观测值统一到相同的时间起点，此时，形变序列的补偿量为辅轨道相对于主轨道观测序列起始时刻的形变量偏差，将之加到序列累积量上即得到辅轨道下第j幅形变观测图中相干目标经整体速率偏差和因起始时间差引起的形变累积量差补偿后的形变即沉降累积量，即

$d_j^{s^{\′\′}}=v_m\·\mathrm{\Δ}T+d_j^{s^{\′}}---\left(9\right)$

式中：为辅轨道下第j幅形变观测图中相干目标经整体速率偏差和因起始时间差引起的形变累积量差补偿后的形变即沉降累积量，为辅轨道下下j幅形变观测图中相干目标经整体速率偏差补偿后的形变累积量，v_m为同名相干目标在主轨道下的沉降速率，ΔT为辅轨道相对主轨道在起始时刻的时间差，单位为年；

在完成相干目标的整体速率偏差补偿和时间偏差补偿后，主辅轨道下的形变累积量序列统一到相同的时间和空间参考基准上，从而得到相干目标加密后的PSInSAR形变即沉降观测序列为

D＝[d₀,d₁,…d_M-1,…,d_M+N-1] (10)

式中，D为相干目标加密后PSInSAR形变即沉降累积量观测序列，d_i(i＝1,2,3,...,M+N-1)为主轨道相干目标或辅轨道经基准补偿后相干目标的PSInSAR形变累积量。