一种无需地面辅助的升降轨InSAR监测沉降区绝对地表形变的方法

摘要

本发明公开了一种无需地面辅助的升降轨InSAR监测沉降区绝对地表形变的方法，首先利用升轨和降轨InSAR数据获取地表在两个不同斜距向上的相对地表形变测量值；然后融合升轨和降轨InSAR相对形变场估计地表在垂直向和东西向上的相对形变；进而利用东西向相对形变量级小于一定阈值的地面点计算升轨和降轨InSAR相对形变测量值的绝对偏移量；最后利用经绝对偏移量改正后的升轨和降轨InSAR形变值估计地表在垂直向和东西向上的绝对形变。突破了传统InSAR只能获取相对地表形变的技术瓶颈，无需布设地面控制点或者假设远场形变为零就可以得到地表在垂直向和东西向上的绝对形变结果。

说明书

技术领域

本发明属于基于遥感影像的大地测量领域，尤其涉及一种无需地面辅助的升降轨InSAR监测沉降区绝对地表形变的方法。

背景技术

合成孔径雷达干涉测量(InSAR)技术是近几十年发展起来的一种空间大地测量和遥感手段，具有监测范围大、测量精度高、空间近似连续、全天时全天候等优势，已经被广泛地用于地震、火山喷发、冰川漂移、滑坡、泥石流等地质灾害引起的地表形变监测。特别是随着多时相InSAR技术(如永久散射体和短基线集技术)的日益成熟，可以对同一轨道获取的时序SAR影像展开联合分析，从而更好的抑制了InSAR干涉图中的地形残差、大气延迟和失相关噪声等，进一步提高了InSAR形变监测的精度和可靠性，使得InSAR技术在监测由于地下流体和矿产资源开采、冻土冻胀和融沉等导致的缓慢长期的地面沉降时也展现出了极大的潜力，详见参考文献[1]和[2]。但是，由于InSAR技术采用的是相位测量方式，在形变解算过程中需要已知或假设其中一个监测点的形变量才能够解算其他监测点上的形变结果，因此InSAR技术得到的地表形变场是相对的，和绝对形变场之间存在一个常数差。

如何将InSAR技术监测得到的地表相对形变转化为绝对形变，目前国际上的主流方法有以下两种：(1)地面控制点法，即在监测区域布设一个或多个地面观测站，利用水准、GPS等手段获取这些地面离散点上的形变值，然后利用这些点上的实测形变值和InSAR形变测量值之差将InSAR的相对形变转化为绝对形变。但是对于没有布设地面控制点的地区，该方法无法适用。此外，InSAR所监测的地面点和实地测量的地面点也往往难以统一，影响了形变转化的精度。(2)远场形变假设法，即在远离地表形变区域的区域寻找稳定的地面点，并假设稳定点上的地表形变为零，从而得到InSAR形变场中其他地面点相对于该稳定点的地表形变结果。很明显，对于形变影响区域较大而导致无法找到稳定点的情况，该方法无法得到地表的绝对形变结果，详见参考文献[3]。

通过上述分析可以看出，对于没有地面控制点的地区，而地表形变的影响范围又超过了整个SAR影像的覆盖范围，那么现有的方法都无法得到地表的绝对形变结果。而在实际中，上述情况是非常常见的，例如在面对地下流体和矿产资源开采导致的大范围地面沉降时，以及冻土由于冻胀和融沉产生季节性垂直形变时，往往会出现难以找到可用的地面控制点或稳定点的情况，导致所监测到的形变结果不能真实反映地面沉降。特别是随着近年来高分SAR数据的普及，InSAR形变场的分辨率越来越高，但覆盖范围却越来越小，使得InSAR技术的这个局限性愈加凸显。

发明内容

本发明的目的在于，克服现有InSAR技术难以监测沉降区绝对地表形变的局限性，提供一种无需地面辅助的升降轨InSAR监测沉降区绝对地表形变的方法。

一种无需地面辅助的升降轨InSAR监测沉降区绝对地表形变的方法，首先利用升轨和降轨InSAR数据获取沉降区地表在两个不同斜距向上的相对地表形变测量值；然后融合升轨和降轨InSAR相对形变场估计地表在垂直向和东西向上的相对形变；接着，利用东西向相对形变量级小于设定阈值的地面点计算升轨和降轨InSAR相对形变测量值的绝对偏移量；最后利用经绝对偏移量修正后的升轨和降轨InSAR形变值估计地表在垂直向和东西向上的绝对形变。

沉降区地表在两个不同斜距向上的相对地表形变测量值是通过收集所监测沉降区的升轨和降轨SAR数据，利用差分InSAR或多时相InSAR技术获取该区域地表在两个不同雷达视线方向上的相对形变测量值，将它们编码至统一的地理坐标系和格网下，并将它们的时间尺度和地面参考点统一；

所述沉降区即地壳运动以垂直向形变为主要特征的区域；

所述雷达视线方向即为斜距向；

所述升轨和降轨InSAR相对形变测量值的绝对偏移量按照以下公式计算：

其中，y_i表示第i个挑选出来的地面观测点，i的取值为1-N，且地表观测点y_i的东西向相对形变不大于设定的形变阈值即ε为事先设定的形变阈值；

代表第i个挑选出来的地面观测点在垂直向上的绝对形变；

K_as和K_des分别代表升轨和降轨InSAR斜距向相对形变与绝对形变测量值之间的常数差，即升轨和降轨InSAR绝对偏移量；

和分别代表地面观测点y_i上升轨和降轨InSAR斜距向相对形变测量值；

和分别为地表观测点y_i的垂直向地表形变在升轨和降轨InSAR斜距向上的投影系数，θ_as和θ_des分别表示地表观测点y_i上的升轨和降轨雷达局部入射角。

利用SAR卫星的极轨特性，解算地表在垂直向和东西向上的相对形变。由于现有SAR卫星都是在近南北的极轨轨道上飞行，因此无论是升轨还是降轨InSAR斜距向测量值都对南北向形变极其不敏感。

为了能够利用升降轨InSAR斜距向形变测量值估计垂直向和东西向形变，通常忽略南北向地表形变对InSAR斜距向形变测量值的贡献：

利用稀疏最小二乘算法解算上权利要求2中所述方程，得到升降轨InSAR绝对偏移量，以及所挑选出来的N个地面观测点在垂直向上的绝对形变。

建立地面观测点x_i上升降轨InSAR斜距向绝对形变测量值与垂直向和东西向绝对形变之间的函数关系，利用最小二乘平差解算出地面观测点在垂直向和东西向上的绝对形变和

其中，和分别代表地面观测点x_i上升轨和降轨InSAR斜距向相对形变测量值；分别为地表观测点x_i的垂直向、东西向地表形变在升轨InSAR斜距向上的投影系数，θ_as(x_i)和α_as(x_i)则分别为地表观测点x_i上的升轨雷达局部入射角和卫星飞行方向角；

分别为地表观测点x_i的垂直向、东西向和南北向地表形变在降轨InSAR斜距向上的投影系数；

θ_des(x_i)和α_des(x_i)则分别为地表观测点x_i上的升轨雷达局部入射角和卫星飞行方向角。

所述的形变阈值ε的取值小于或等于3mm，以保证忽略东西向形变的假设具有合理性。

所述沉降区表示地壳运动以垂直向形变为主要特征的区域，包括由于地下水开采、地下油气开采、地下采矿、地热开采、地下工程、冻土冻胀和融沉等引起地表形变的区域。

所述的升轨和降轨InSAR斜距向相对形变测量值由单一SAR卫星数据或联合不同SAR卫星数据获取。

有益效果

本发明提供了一种无需地面辅助的升降轨InSAR监测沉降区绝对地表形变的方法，1)利用InSAR技术获取沉降区地理编码后的升轨和降轨斜距向地表形变场的相对测量值，并统一它们的格网、时间尺度和地面参考点；2)根据SAR卫星的成像几何建立升降轨InSAR斜距向相对形变测量值与三维相对形变之间的函数关系；3)根据SAR卫星的极轨特性忽略南北向形变对InSAR斜距向测量值的贡献，估计垂直向和东西向上的相对形变；4)针对东西向相对形变小于一定阈值的地面观测点估计升轨和降轨InSAR斜距向相对形变测量值与绝对形变测量值之间的绝对偏移量；5)利用绝对偏移量将升降轨InSAR相对形变测量值转化为绝对形变测量值，并在此基础上估计垂直向和东西向上的绝对形变结果。该方法实现简单，无需在地面布设控制点或者寻找稳定区域，对地下开采、冻土活动等导致的绝对地表形变监测而言是一种高效率、低成本、高精度和大范围的方法，突破了InSAR技术本身只能监测相对地表形变的局限，进一步推动了InSAR大地测量技术的市场化和工程化进程。

附图说明

图1是升轨和降轨SAR卫星的成像几何图；

图2是本发明所述方法的流程图；

图3是模拟的沉降区地表形变，其中，(a)为垂直向形变；(b)为东西向形变；(c)为南北向形变；单位：cm/yr；

图4是含噪的InSAR斜距向相对形变测量值，其中，(a)为升轨测量值；(b)为降轨测量值；单位：cm/yr；

图5是垂直向和东西向上的相对形变结果；单位：cm/yr；

图6是东西向相对形变小于阈值的地面观测点分布；

图7是本发明得到的垂直向和东西向上的绝对形变结果；其中，(a)是垂直向绝对地表形变，(b)是东西向绝对地表形变，单位：cm/yr；

图8为应用本发明得到的绝对形变与模拟形变之间的差值；(a)是垂直向绝对和模拟的地表形变之间的差值，(b)为是东西向绝对和模拟的地表形变之间的差值；单位：cm/yr。

具体实施方式

下面将结合附图和实施例对本发明做进一步的说明。

为了便于理解本发明，首先提供本发明的理论基础：

众所周知，地质灾害的直接后果是地表发生形变。而对于InSAR技术而言，其只能监测地面观测点x_i在斜距向上相对于某个参考点的相对形变，即

D^rat(x_i)＝D^real(x_i)+K(1)>

其中，D^rat(x_i)为地面观测点x_i上的InSAR斜距向相对形变测量值；D^real(x_i)为地面观测点x_i上的斜距向绝对形变；K为InSAR斜距向相对形变测量值与斜距向绝对形变之间的常数差，即绝对偏移量。

然而在实际中，地表形变通常是发生在三维框架中的，即会同时发生垂直向、东西向和南北向形变。因此，斜距向绝对形变只是地表真实三维形变在斜距向上的投影，两者之间的关系可由下式表达：

其中，和分别为地表真实三维形变；S^u(x_i)、S^e(x_i)和Sⁿ(x_i)分别为地表观测点x_i的垂直向、东西向和南北向地表形变在InSAR斜距向上的投影系数：

其中，θ(x_i)和α(x_i)则分别为地表观测点x_i上的雷达局部入射角和卫星飞行方向角(以北方向为起始顺时针旋转)。升轨和降轨SAR卫星的成像几何图如图1所示。由于目前的SAR卫星均在太阳同步轨道(即极轨)上飞行，因此α(x_i)通常在350度(升轨)和190度(降轨)左右。而作为侧视成像雷达，SAR的雷达局部入射角一般在20度到50度之间变化。因此，相应的投影系数变化范围分别大约为：0.7≤S^u(x_i）≤0.9、0.3≤S^e(x_i)≤0.6、0.08≤Sⁿ(x_i)≤0.10，从而导致InSAR的斜距向形变测量值对垂直向地表形变最为敏感，东西向地表形变次之，而南北向地表形变最差，详见参考文献[4]和[5]。

如图2所示，一种无需地面辅助的升降轨InSAR监测沉降区绝对地表形变的方法，包括以下步骤：

(1)利用差分InSAR或多时相InSAR技术获取待监测沉降区地表分别在升轨和降轨雷达视线方向(即斜距向)上的相对形变测量值，并对它们进行地理编码，并统一它们的时间尺度和地面参考点；

(2)利用SAR影像头文件中包含的雷达局部入射角、卫星飞行方位角，按照公式(3)计算每个地面观测点上的投影系数S^u(x_i)、S^e(x_i)和Sⁿ(x_i)。然后根据公式(2)就可以构建升轨和降轨InSAR斜距向相对形变测量值与三维相对地表形变之间的函数模型：

其中为升轨和降轨InSAR斜距向相对形变测量值的误差；

(3)忽略南北向地表形变对InSAR斜距向形变测量值的贡献，则公式(4)可以写成如下形式：

根据最小二乘平差原理就可以解算出垂直向和东西向上的相对形变：

(4)挑选的地面观测点，其中ε为事先设定的形变阈值，一般不超过3mm/yr。对于地下开采、冻土活动等导致的地表形变而言，垂直形变较大的地方往往水平形变较小(如沉降漏斗中心)，而垂直形变较小的地方往往水平形变较大(如沉降漏斗边缘)。因此，可以假设挑选出来的N个地面观测点上的升轨和降轨InSAR斜距向相对形变测量值主要由垂直向形变所贡献，因此忽略东西向形变对其的贡献，从而建立以下函数模型：

其中，为挑选出来的N个地面观测点上升轨和降轨InSAR斜距向相对形变测量值，为相应的观测误差；K_as和K_des分别代表升轨和降轨InSAR斜距向相对形变与绝对形变测量值之间的常数差，即升轨和降轨InSAR绝对偏移量。

上式中观测量的个数为2N个，而待求参数的个数为N+2个，因此观测量的个数大于待求参数的个数，可以通过最小二乘平差对其进行解算。但是由于设计矩阵为一个大型的稀疏矩阵，因此需要利用稀疏最小二乘对其进行求解，从而得到升轨和降轨InSAR绝对偏移量，以及所挑选出来的N个地面观测点在垂直向上的绝对形变。

(5)根据公式(1)将升降轨InSAR斜距向相对形变测量值转化为绝对形变测量值：

再根据公式(2)，同样忽略南北向形变对InSAR斜距向形变测量值的贡献，构建地面观测点x_i上升降轨InSAR斜距向绝对形变测量值与垂直向和东西向绝对形变之间的函数模型：

其中为升轨和降轨InSAR斜距向绝对形变测量值的误差。同样利用最小二乘平差原理就可以解算出垂直向和东西向上的绝对形变：

在400×450的规则格网中模拟沉降区的三维绝对地表形变，格网尺寸为10m×10m，其中垂直向、东西向和南北向的绝对形变分别如图3(a)、3(b)和3(c)所示。然后利用公式(2)模拟出升轨和降轨InSAR斜距向相对形变测量值，结果分别如图4(a)和4(b)所示，其中的三角形代表地面参考点。为了让模拟数据具有真实性，升轨数据的局部雷达入射角和卫星飞行方位角采用ALOS/PALSAR卫星影像头文件中提供的参数，并将均值为零、标准偏差为2mm的高斯白噪声加入到升轨InSAR斜距向形变测量值；降轨数据的局部雷达入射角和卫星飞行方位角则采用ENVISAT/ASAR卫星影像头文件中提供的参数，并将均值为零、标准偏差为1.5mm的高斯白噪声加入到降轨InSAR斜距向形变测量值。由于升轨和降轨模拟数据是直接模拟的同一地理坐标系下的年均形变速率结果，并且采用同一个参考点，因此在这次试验中无需时间尺度统一、地理编码、参考点统一等步骤。

通过本发明所提出的方法处理，就可以利用上述模拟的含噪升降轨InSAR斜距向相对形变测量值估计出垂直向和东西向上的绝对地表形变。图5(a)和5(b)分别是利用公式(6)解算出的地表在垂直向和东西向上的相对地表形变。可以看出，虽然在东西向上该结果与原始模拟的东西向绝对地表形变比较一致，但在垂直向上该结果的变化区间与原始模拟的垂直向绝对地表形变相差较大，这也说明传统的InSAR方法难以得到可靠、真实的地面沉降结果。图6是本次试验挑选出来的东西向相对形变小于3mm的地面观测点，总共有48833个。图7(a)和7(b)分别是利用公式(10)解算出的地表在垂直向和东西向上的绝对地表形变，总体而言和原始模拟的垂直向和东西向形变非常一致。图8(a)和8(b)给出的分别是该方法解算出来的垂直向和东西向形变场与模拟的垂直向和东西向形变场之间的差值。为了定量验证本发明的效果，实施例中分别计算垂直向和东西向形变的均方根误差，分别为2.1mm和2.6mm，与InSAR斜距测量值中噪声的标准差相当，从而说明本发明是可行的，能够得到沉降区可靠的地表垂直向和东西向绝对形变结果。

参考文献：

[1]Berardino,P.,Fornaro,G.,Lanari,R.,Sansosti,E.,2002.Anew algorithm for surface deformation monitoring based on small baseline differential SAR interferograms.IEEE Transactions on Geoscience and Remote Sensing 40(11),2375-2383；

[2]Ferretti,A.,Prati,C.,Rocca,F.,2001.Permanent scatterers in SAR interferometry.IEEE Transactions on Geoscience and Remote Sensing 39(1),8-20；

[3]Liu,L.,Zhang,T.J.,Wahr,J.,2010.InSAR measurements of surface deformation over permafrost on the North Slope of Alaska.Journal of Geophysical Research115,F03023；

[4]Fialko,Y.,Sandwell,D.,Simons,M.,Rosen,P.,2005.Three-dimensional deformation caused by the Bam,Iran,earthquake and the origin of shallow slip deficit.Nature 435,295-299；

[5]Hu,J.,Li,Z.W.,Ding,X.L.,Zhu,J.J.,Zhang,L.,Sun,Q.,2014.Resolving Three-Dimensional Surface Displacements from InSAR Measurements:AReview.Earth-Science Reviews 133,1-17。