一种无需地面辅助的升降轨InSAR监测沉降区绝对地表形变的方法
【专利摘要】本发明公开了一种无需地面辅助的升降轨InSAR监测沉降区绝对地表形变的方法，首先利用升轨和降轨InSAR数据获取地表在两个不同斜距向上的相对地表形变测量值；然后融合升轨和降轨InSAR相对形变场估计地表在垂直向和东西向上的相对形变；进而利用东西向相对形变量级小于一定阈值的地面点计算升轨和降轨InSAR相对形变测量值的绝对偏移量；最后利用经绝对偏移量改正后的升轨和降轨InSAR形变值估计地表在垂直向和东西向上的绝对形变。突破了传统InSAR只能获取相对地表形变的技术瓶颈，无需布设地面控制点或者假设远场形变为零就可以得到地表在垂直向和东西向上的绝对形变结果。
【专利说明】
一种无需地面辅助的升降轨I nSAR监测沉降区绝对地表形变 的方法
技术领域
[0001] 本发明属于基于遥感影像的大地测量领域，尤其涉及一种无需地面辅助的升降轨 InSAR监测沉降区绝对地表形变的方法。
【背景技术】
[0002] 合成孔径雷达干涉测量（InSAR)技术是近几十年发展起来的一种空间大地测量和 遥感手段，具有监测范围大、测量精度高、空间近似连续、全天时全天候等优势，已经被广泛 地用于地震、火山喷发、冰川漂移、滑坡、泥石流等地质灾害引起的地表形变监测。特别是随 着多时相InSAR技术(如永久散射体和短基线集技术）的日益成熟，可以对同一轨道获取的 时序SAR影像展开联合分析，从而更好的抑制了 InSAR干涉图中的地形残差、大气延迟和失 相关噪声等，进一步提高了 InSAR形变监测的精度和可靠性，使得InSAR技术在监测由于地 下流体和矿产资源开采、冻土冻胀和融沉等导致的缓慢长期的地面沉降时也展现出了极大 的潜力，详见参考文献[1]和[2]。但是，由于InSAR技术采用的是相位测量方式，在形变解算 过程中需要已知或假设其中一个监测点的形变量才能够解算其他监测点上的形变结果，因 此InSAR技术得到的地表形变场是相对的，和绝对形变场之间存在一个常数差。
[0003] 如何将InSAR技术监测得到的地表相对形变转化为绝对形变，目前国际上的主流 方法有以下两种：（1)地面控制点法，即在监测区域布设一个或多个地面观测站，利用水准、 GPS等手段获取这些地面离散点上的形变值，然后利用这些点上的实测形变值和InSAR形变 测量值之差将InSAR的相对形变转化为绝对形变。但是对于没有布设地面控制点的地区，该 方法无法适用。此外，InSAR所监测的地面点和实地测量的地面点也往往难以统一，影响了 形变转化的精度。（2)远场形变假设法，即在远离地表形变区域的区域寻找稳定的地面点， 并假设稳定点上的地表形变为零，从而得到InSAR形变场中其他地面点相对于该稳定点的 地表形变结果。很明显，对于形变影响区域较大而导致无法找到稳定点的情况，该方法无法 得到地表的绝对形变结果，详见参考文献[3]。
[0004] 通过上述分析可以看出，对于没有地面控制点的地区，而地表形变的影响范围又 超过了整个SAR影像的覆盖范围，那么现有的方法都无法得到地表的绝对形变结果。而在实 际中，上述情况是非常常见的，例如在面对地下流体和矿产资源开采导致的大范围地面沉 降时，以及冻土由于冻胀和融沉产生季节性垂直形变时，往往会出现难以找到可用的地面 控制点或稳定点的情况，导致所监测到的形变结果不能真实反映地面沉降。特别是随着近 年来高分SAR数据的普及，InSAR形变场的分辨率越来越高，但覆盖范围却越来越小，使得 InSAR技术的这个局限性愈加凸显。
【发明内容】
[0005] 本发明的目的在于，克服现有InSAR技术难以监测沉降区绝对地表形变的局限性， 提供一种无需地面辅助的升降轨InSAR监测沉降区绝对地表形变的方法。
[0006] -种无需地面辅助的升降轨InSAR监测沉降区绝对地表形变的方法，首先利用升 轨和降轨InSAR数据获取沉降区地表在两个不同斜距向上的相对地表形变测量值;然后融 合升轨和降轨InSAR相对形变场估计地表在垂直向和东西向上的相对形变;接着，利用东西 向相对形变量级小于设定阈值的地面点计算升轨和降轨InSAR相对形变测量值的绝对偏移 量;最后利用经绝对偏移量修正后的升轨和降轨InSAR形变值估计地表在垂直向和东西向 上的绝对形变。
[0007] 沉降区地表在两个不同斜距向上的相对地表形变测量值是通过收集所监测沉降 区的升轨和降轨SAR数据，利用差分InSAR或多时相InSAR技术获取该区域地表在两个不同 雷达视线方向上的相对形变测量值，将它们编码至统一的地理坐标系和格网下，并将它们 的时间尺度和地面参考点统一；
[0008] 所述沉降区即地壳运动以垂直向形变为主要特征的区域；
[0009] 所述雷达视线方向即为斜距向；
[0010] ............... .....................-照以下公式计算：
[0011]
[0012 ]其中，yi表示第i个挑选出来的地面观测点，i的取值为I-N，且地表观测点yi的东西 向相对形变不大于设定的形变阈值即(Ji) S ￡，ε为事先设定的形变阈值；
[0013] Afa)代表第i个挑选出来的地面观测点在垂直向上的绝对形变；
[0014] Kas和Kde3s分别代表升轨和降轨InSAR斜距向相对形变与绝对形变测量值之间的常 数差，即升轨和降轨InSAR绝对偏移量；
[0015] 分别代表地面观测点yi上升轨和降轨InSAR斜距向相对形变测量 值；
[0016] 和S=(Ji)分别为地表观测点Yl的垂直向地表形变在升轨和降轨 InSAR斜距 向上的投影系数，-?%) = cos6?as(y;)，= cos&es(%)，0￡^和0(^分别表示地表观测点7 1 上的升轨和降轨雷达局部入射角。
[0017] 利用SAR卫星的极轨特性，解算地表在垂直向和东西向上的相对形变。由于现有 SAR卫星都是在近南北的极轨轨道上飞行，因此无论是升轨还是降轨InSAR斜距向测量值都 对南北向形变极其不敏感。
[0018] 为了能够利用升降轨InSAR斜距向形变测量值估计垂直向和东西向形变，通常忽 略南北向地表形变对InSAR斜距向形变测量值的贡献：
[0019]
[0020] 利用稀疏最小二乘算法解算上权利要求2中所述方程，得到升降轨InSAR绝对偏移 量，以及所挑选出来的N个地面观测点在垂直向上的绝对形变。
[0021] 建立地面观测点Xi上升降轨InSAR斜距向绝对形变测量值与垂直向和东西向绝对 形变之间的函数关系，利用最小二乘平差解算出地面观测点在垂直向和东西向上的绝对形 变β妒⑷和職齡
[0022]
[0023] 其中，⑷分别代表地面观测点Xi上升轨和降轨InSAR斜距向相对形 变测量值；Cp (Xi)、(Xi)分别为地表观测点Xi的垂直向、东西向地表形变在升轨InSAR斜 距向上的投影系数，CCO = ?)成.,;(々），= -cos(aa.s(x,.) - 3:^ aas(Xl)则分别为地表观测点X1上的升轨雷达局部入射角和卫星飞行方向角；
[0024] 以(幻、Sgf (Xi)分别为地表观测点Xi的垂直向、东西向和南北向地表形变在降轨 InSAR斜距向上的投影系数；
[0025]
则分别为地表观测点X1上的升轨雷达局部入射角和卫星飞行方向角。
[0026]所述的形变阈值ε的取值小于或等于3mm，以保证忽略东西向形变的假设具有合理 性。
[0027] 所述沉降区表示地壳运动以垂直向形变为主要特征的区域，包括由于地下水开 采、地下油气开采、地下采矿、地热开采、地下工程、冻土冻胀和融沉等引起地表形变的区 域。
[0028]所述的升轨和降轨InSAR斜距向相对形变测量值由单一 SAR卫星数据或联合不同 SAR卫星数据获取。
[0029]有益效果
[0030] 本发明提供了一种无需地面辅助的升降轨InSAR监测沉降区绝对地表形变的方 法，1)利用InSAR技术获取沉降区地理编码后的升轨和降轨斜距向地表形变场的相对测量 值，并统一它们的格网、时间尺度和地面参考点；2)根据SAR卫星的成像几何建立升降轨 InSAR斜距向相对形变测量值与三维相对形变之间的函数关系；3)根据SAR卫星的极轨特性 忽略南北向形变对InSAR斜距向测量值的贡献，估计垂直向和东西向上的相对形变;4)针对 东西向相对形变小于一定阈值的地面观测点估计升轨和降轨InSAR斜距向相对形变测量值 与绝对形变测量值之间的绝对偏移量;5)利用绝对偏移量将升降轨InSAR相对形变测量值 转化为绝对形变测量值，并在此基础上估计垂直向和东西向上的绝对形变结果。该方法实 现简单，无需在地面布设控制点或者寻找稳定区域，对地下开采、冻土活动等导致的绝对地 表形变监测而言是一种高效率、低成本、高精度和大范围的方法，突破了InSAR技术本身只 能监测相对地表形变的局限，进一步推动了InSAR大地测量技术的市场化和工程化进程。
【附图说明】
[0031] 图1是升轨和降轨SAR卫星的成像几何图；
[0032]图2是本发明所述方法的流程图；
[0033] 图3是模拟的沉降区地表形变，其中，（a)为垂直向形变；（b)为东西向形变；（c)为 南北向形变;单位:cm/yr;
[0034] 图4是含噪的InSAR斜距向相对形变测量值，其中，（a)为升轨测量值；（b)为降轨测 量值;单位:cm/yr;
[0035]图5是垂直向和东西向上的相对形变结果;单位:cm/yr;
[0036] 图6是东西向相对形变小于阈值的地面观测点分布；
[0037] 图7是本发明得到的垂直向和东西向上的绝对形变结果；其中，（a)是垂直向绝对 地表形变，（b)是东西向绝对地表形变，单位:cm/yr;
[0038] 图8为应用本发明得到的绝对形变与模拟形变之间的差值；（a)是垂直向绝对和模 拟的地表形变之间的差值，（b)为是东西向绝对和模拟的地表形变之间的差值;单位:cm/ yr〇
【具体实施方式】
[0039]下面将结合附图和实施例对本发明做进一步的说明。
[0040] 为了便于理解本发明，首先提供本发明的理论基础：
[0041] 众所周知，地质灾害的直接后果是地表发生形变。而对于InSAR技术而言，其只能 监测地面观测Ax1在斜距向上相对于某个参考点的相对形变，即
[0042]
[0043] 其中，Drat(X1)为地面观测点X1上的InSAR斜距向相对形变测量值;D lreal(X1)为地面 观测点Xi上的斜距向绝对形变;K为InSAR斜距向相对形变测量值与斜距向绝对形变之间的 常数差，即绝对偏移量。
[0044] 然而在实际中，地表形变通常是发生在三维框架中的，即会同时发生垂直向、东西 向和南北向形变。因此，斜距向绝对形变只是地表真实三维形变在斜距向上的投影，两者之 间的关系可由下式表达：
[0045]
[0046] 其中，A= fe)、呢?丨(々)和〇=^0^)分别为地表真实三维形变;Su(Xi)、Se( XiWPSn (X1)分别为地表观测点^的垂直向、东西向和南北向地表形变在InSAR斜距向上的投影系 数：
[0047 (3)
[0048]其中，θ(Χι)和α(Χι)则分别为地表观测点^上的雷达局部入射角和卫星飞行方向 角（以北方向为起始顺时针旋转）。升轨和降轨SAR卫星的成像几何图如图1所示。由于目前 的SAR卫星均在太阳同步轨道（即极轨）上飞行，因此Ct(X 1)通常在350度(升轨）和190度（降 轨）左右。而作为侧视成像雷达，SAR的雷达局部入射角一般在20度到50度之间变化。因此， 相应的投影系数变化范围分别大约为:0.7彡S u(Xi)彡0.9、0.3彡Se(Xi)彡0.6、0.08彡S n(Xi) 彡0.10，从而导致InSAR的斜距向形变测量值对垂直向地表形变最为敏感，东西向地表形变 次之，而南北向地表形变最差，详见参考文献[4]和[5]。
[0049] 如图2所示，一种无需地面辅助的升降轨InSAR监测沉降区绝对地表形变的方法， 包括以下步骤：
[0050] (1)利用差分InSAR或多时相InSAR技术获取待监测沉降区地表分别在升轨和降轨 雷达视线方向（即斜距向）上的相对形变测量值，并对它们进行地理编码，并统一它们的时 间尺度和地面参考点；
[0051] (2)利用SAR影像头文件中包含的雷达局部入射角、卫星飞行方位角，按照公式(3) 计算每个地面观测点上的投影系数Su(X1)、Se(XlWPS n(Xl)。然后根据公式⑵就可以构建升 轨和降轨I n S A R斜距向相对形变测量值[0Γ (?), (?)]7与三维相对地表形变
[D#⑷几怒(4)』=丨办)广之间的函数模型：
[0052]
C4)
[0053] 其中[私f (Xi),7?認(χ0Γ为升轨和降轨InSAR斜距向相对形变测量值的误差；
[0054] (3)忽略南北向地表形变对InSAR斜距向形变测量值的贡献，则公式(4)可以写成 如下形忒，
[0055；
(5)
[0056] 根据最小二乘平差原?
nin就可以解算出垂直向和东西向上 的相对形变：
[0057]
[0058] (4)挑选Wife) S ε的地面观测点，其中￡为事先设定的形变阈值，一般不超过 3mm/yr。对于地下开采、冻土活动等导致的地表形变而言，垂直形变较大的地方往往水平形 变较小（如沉降漏斗中心），而垂直形变较小的地方往往水平形变较大（如沉降漏斗边缘）。 因此，可以假设挑选出来的N个地面观测点上的升轨和降轨InSAR斜距向相对形变测量值主 要由垂直向形变所贡献，因此忽略东西向形变对其的贡献，从而建立以下函数模型：
[0059]
[0060] 其中，(凡),…,^^(^,^^(^,…^^(^，为挑选出来的^^个地面观测点上升 轨和降轨InSAR斜距向相对形变测量值，[^y 1),…,Cy1),…,碟丨(3^)Γ为相应的 观测误差;Kas和K de3s分别代表升轨和降轨InSAR斜距向相对形变与绝对形变测量值之间的常 数差，即升轨和降轨InSAR绝对偏移量。
[0061 ] 上式中观测量的个数为2N个，而待求参数的个数为N+2个，因此观测量的个数大于 待求参数的个数，可以通过最小二乘平差对其进行解算。但是由于设计矩阵为一个大型的 稀疏矩阵，因此需要利用稀疏最小二乘对其进行求解，从而得到升轨和降轨InSAR绝对偏移 量，以及所挑选出来的N个地面观测点在垂直向上的绝对形变。
[0062] (5)根据公式(1)将升降轨InSAR斜距向相对形变测量值转化为绝对形变测量值：
[0063]
(8)
[0066] 其中[成& ~Vf为孙轨和隆轨InSAR斜距向绝对形变测量值的误差。同样
[0064] 再根据公式(2)，同样忽略南北向形变对InSAR斜距向形变测量值的贡献，构建地 面观测点Xi上升降轨InSAR斜距向绝对形变测量值与垂直向和东西向绝对形变之间的函数 模型：
[0065] (9)
利用最小二乘平差 rain就可以解算出垂直向和东西向上的绝对 形变：
[0067]
[0068] 在400 X 450的规则格网中模拟沉降区的三维绝对地表形变，格网尺寸为I Om X l〇m，其中垂直向、东西向和南北向的绝对形变分别如图3(a)、3(b)和3(c)所示。然后利用公 式(2)模拟出升轨和降轨InSAR斜距向相对形变测量值，结果分别如图4(a)和4(b)所示，其 中的三角形代表地面参考点。为了让模拟数据具有真实性，升轨数据的局部雷达入射角和 卫星飞行方位角采用AL0S/PALSAR卫星影像头文件中提供的参数，并将均值为零、标准偏差 为2mm的高斯白噪声加入到升轨InSAR斜距向形变测量值;降轨数据的局部雷达入射角和卫 星飞行方位角则采用ENVISAT/ASAR卫星影像头文件中提供的参数，并将均值为零、标准偏 差为1.5_的高斯白噪声加入到降轨InSAR斜距向形变测量值。由于升轨和降轨模拟数据是 直接模拟的同一地理坐标系下的年均形变速率结果，并且采用同一个参考点，因此在这次 试验中无需时间尺度统一、地理编码、参考点统一等步骤。
[0069] 通过本发明所提出的方法处理，就可以利用上述模拟的含噪升降轨InSAR斜距向 相对形变测量值估计出垂直向和东西向上的绝对地表形变。图5(a)和5(b)分别是利用公式 (6)解算出的地表在垂直向和东西向上的相对地表形变。可以看出，虽然在东西向上该结果 与原始模拟的东西向绝对地表形变比较一致，但在垂直向上该结果的变化区间与原始模拟 的垂直向绝对地表形变相差较大，这也说明传统的InSAR方法难以得到可靠、真实的地面沉 降结果。图6是本次试验挑选出来的东西向相对形变小于3mm的地面观测点，总共有48833 个。图7(a)和7(b)分别是利用公式（10)解算出的地表在垂直向和东西向上的绝对地表形 变，总体而言和原始模拟的垂直向和东西向形变非常一致。图8(a)和8(b)给出的分别是该 方法解算出来的垂直向和东西向形变场与模拟的垂直向和东西向形变场之间的差值。为了 定量验证本发明的效果，实施例中分别计算垂直向和东西向形变的均方根误差，分别为 2. Imm和2.6mm，与InSAR斜距测量值中噪声的标准差相当，从而说明本发明是可行的，能够 得到沉降区可靠的地表垂直向和东西向绝对形变结果。
[0070] 参考文献：
[0071 ] [l]Berardino,P.,FornarojG.,Lanari,R.,Sansosti,Ε.,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；
[0072] [2]Ferretti，A.，Prati，C.，Rocca，F.，2001.Permanent scatterers in SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing 39(1),8-20；
[0073] [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 Researchll5,F03023；
[0074] [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;
[0075] [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〇
【主权项】
1. 一种无需地面辅助的升降轨InSAR监测沉降区绝对地表形变的方法，其特征在于，首 先利用升轨和降轨InSAR数据获取沉降区地表在两个不同斜距向上的相对地表形变测量 值;然后融合升轨和降轨InSAR相对形变场估计地表在垂直向和东西向上的相对形变;接 着，利用东西向相对形变量级小于设定阔值的地面点计算升轨和降轨InSAR相对形变测量 值的绝对偏移量;最后利用经绝对偏移量修正后的升轨和降轨InSAR形变值估计地表在垂 直向和东西向上的绝对形变。2. 根据权利要求1所述的方法，其特征在于，所述升轨和降轨InSAR相对形变测量值的 绝对偏移量按照W下公式计算：其中，y康示第i个挑选出来的地面观测点，i的取值为1-N，且地表观测点yi的东西向相 对形变不大于设定的形变阔值即K芯t (乂·）< f，ε为事先设定的形变阔值； 明严从)代表第i个挑选出来的地面观测点在垂直向上的绝对形变； Kas和Kdes分别代表升轨和降轨InSAR斜距向相对形变与绝对形变测量值之间的常数差， 即升轨和降轨InSAR绝对偏移量； 化f从)和0益从)分别代表地面观测点yi上升轨和降轨InSAR斜距向相对形变测量值； 端P (如和端;(知分别为地表观测点yi的垂直向地表形变在升轨和降轨InSAR斜距向上 的投影系数S和0des分别表示地表观测点yi上的 升轨和降轨雷达局部入射角。3. 根据权利要求2所述的方法，其特征在于，利用稀疏最小二乘算法解算上权利要求2 中所述方程，得到升降轨InSAR绝对偏移量，W及所挑选出来的N个地面观测点在垂直向上 的绝对形变。4. 根据权利要求3所述的方法，其特征在于，建立地面观测点XI上升降轨InSAR斜距向绝 对形变测量值与垂直向和东西向绝对形变之间的函数关系，利用最小二乘平差解算出地面 观测点在垂直向和东西向上的绝对形变斯戸(Χ?)和巧盘(Χ?);其中，Kf (切和D栽(切分别代表地面观测点χι上升轨和降轨InSAR斜距向相对形变测 量值;义'f(Yi)、S晋"挺)分别为地表观测点xi的垂直向、东西向地表形变在升轨InSAR斜距向 上的投影系数和口 as (Xi)则分别为地表观测点XI上的升轨雷达局部入射角和卫星飞行方向角； 《左耗)、S器f(如分别为地表观测点XI的垂直向、东西向和南北向地表形变在降轨 InSAR斜距向上的投影系数；和Odes (Xi)则分 别为地表观测点Xi上的升轨雷达局部入射角和卫星飞行方向角。5. 根据权利要求2-4任一项所述的方法，其特征在于，所述的形变阔值ε的取值小于或 等于3mm。6. 根据权利要求1所述的方法，其特征在于，所述沉降区表示地壳运动W垂直向形变为 主要特征的区域，包括由于地下水开采、地下油气开采、地下采矿、地热开采、地下工程、冻 ±冻胀和融沉等引起地表形变的区域。7. 根据权利要求1所述的方法，其特征在于，所述的升轨和降轨InSAR斜距向相对形变 测量值由单一 SAR卫星数据或联合不同SAR卫星数据获取。
【文档编号】G01S13/90GK105938193SQ201610554394
【公开日】2016年9月14日
【申请日】2016年7月14日
【发明人】胡俊, 孙倩, 李志伟, 朱建军
【申请人】中南大学