一种利用InSAR反演地下流体体积变化和三维地表形变的方法
【技术领域】
[0001] 本发明属于基于遥感影像的大地测量领域，尤其涉及一种利用InSAR反演地下流 体体积变化和三维地表形变的方法。
【背景技术】
[0002] 差分 InSAR (Differential InSAR，D-InSAR)技术是目前国际上在 InSAR 应用上 最为成熟的技术，它最主要的目的就是监测地球表面厘米级甚至毫米级的形变。多时相 InSAR (Multi-Temporal InSAR, MT-InSAR)技术则是近20多年来在D-InSAR技术的基础上 发展起来的一种地表形变监测方法，如PS、SBAS和TCP等。通过对多景SAR影像的联合分 析，MT-InSAR技术可以更好地抑制干涉图中时空失相关和大气噪声的影响，并能提供地表 在SAR影像获取时刻的形变序列。然而，到目前为止，无论是D-InSAR还是MT-InSAR技术都 只能利用单一平台、单一轨道的SAR数据，因此只能监测地表在雷达视线方向（即斜距向） 上的一维形变结果。但是在现实中，地表形变是发生在三维空间框架下的，即所谓的三维形 变场。因此InSAR的斜距向形变测量结果不能反映出真实的地表形变。
[0003] 如何将InSAR技术监测的一维形变拓展为三维，目前国际上已经有一些学者进行 了探索性的研究，大致可以分为以下几类：（1)多方向InSAR斜距向形变测量值融合法，即 利用三个或以上的InSAR斜距向形变测量值反演三维地表形变场，但由于目前的SAR卫星 都是极轨飞行轨道，这种方法只在高炜度地区适用；（2)升降轨InSAR斜距向和方位向形变 测量值融合法，即融合D-InSAR提供的升轨、降轨斜距向形变测量值和Offset-Tracking或 MAI提供的升轨、降轨方位向形变测量值反演三维地表形变场，但受限于Offset-Tracking 或MAI技术的分米到米级的测量精度，这种方法目前只适用于监测地震、火山喷发和冰川 移动等大型地表形变；（3) InSAR和GPS融合法，即融合D-InSAR或MT-InSAR提供的斜距向 测形变量值和GPS提供的离散点上的三维形变测量值反演三维地表形变场，但这种方法的 精度强烈依赖于GPS台站的密度和分布，在GPS台站较少的地区不适用。
[0004] 而在现实中，往往存在着一些缓慢的地表形变需要进行大范围、全方位和长期的 监测，如地下水抽取、油气田开采和岩浆活动等地下流体运动引起的地表形变。这些地表形 变已经成为影响区域经济和社会可持续发展的重要因素。例如，由于长期过度集中开采地 下水资源，我国迄今为止已有70多个城市发生了不同程度的地面沉降，其中上海和天津市 区的最大沉降已超过两米，严重制约了城市社会经济的发展，同时对人们的生活也构成了 极大的威胁。然而根据上述分析，目前InSAR技术难以广泛地用于监测由地下流体运动引 起的三维地表形变。
【发明内容】
[0005] 本发明的目的在于，克服现有InSAR技术在监测地下流体运动引起的地表形变的 不足和局限性，提供一种利用InSAR斜距向形变测量值监测三维地表形变，并同步反演出 地下流体体积变化量的方法。
[0006] 一种利用InSAR反演地下流体体积变化和三维地表形变的方法，包括以下步骤：
[0007] 步骤1 :利用D-InSAR或多时相InSAR技术获取待监测的受地下流体运动影响的 地区的地表在雷达视线方向的形变测量值，并对其进行地理编码；
[0008] 所述雷达视线方向即为斜距向；
[0009] 步骤2 :根据SAR卫星成像几何，按照以下公式构建InSAR斜距向形变测量值I (X1) 和三维地表形变之间的函数模型：
[0010]
[0011] 其中，X1表示地表观测点，i = 1，2,…，M，共有M个地表观测点；
[0012] S1(X1)、S2(X1)及S3(X 1)分另Ij为地表观测点X1的东西、南北和垂直向地表形变在 InSAR 斜距向上的投影系数，S1(Xi) = -cos (α -3 π /2) sin Θ，S2(Xi) = -sin (α -3 π /2) sin θ，S3(Xi) = cos θ ; θ和α分别为雷达局部入射角和卫星飞行方向角；
[0013] Cl1(X1)为地表观测点位置X1的形变量，1 = 1，2, 3分别代表东西、南北和垂直向上 的三个分量；
[0014] n (X1)为 InSAR 观测误差；
[0015] 【该值很小，通过最小二乘平差可以使其对InSAR斜距向形变测量值的影响最 小；】
[0016] 步骤3 :根据弹性半空间理论建立每个观测点上的三维地表形变与地下流体的体 积变化之间的函数模型：
[0017]
[0018] 其中，G1 (X，y)为格林函数，
V为泊松比，S为地下块源y 到地表观测点X之间的距离，共有N个块源(X)和P1 (y)分别为地表观测点X和块源y 的三维空间位置，
为地下块源y到地 面观测点X之间的距离；Dv(y)为地下半空间体积V中块源y的流体体积变化量；Vy表示单 位块源的体积，由单位块源的尺寸和地下流体的厚度计算获得；ε Jx1)表示模型与真实形 变之间的残差；
[0019] 步骤4 :利用SAR卫星成像几何将步骤3得到的模型转化为地表每个观测点上的 InSAR斜距向形变测量值与地下流体的体积变化之间的函数关系：
[0020]
[0021] 步骤5 :令地下流体中的所有块源的体积变化AV(y)相同，利用地下流体场的先 验信息确定泊松比V和地下流体的层数，利用M个地表观测点上的InSAR斜距向形变测量 值，通过非线性方法反演出块源的体积变化AV(y)、地下流体的深度h和厚度t ;
[0022] 步骤6 :构建以所有地表观测点的三维形变量和所有地下块源的体积变化量为未 知参数的联合模型：
[0023]
[0024] 其中Ω为4MX 1的观测矩阵，由M个地理编码后的InSAR斜距向形变测量值和3M 个虚拟观测量组成：Ω = [I (X1)…I (χΜ) O O O ......... O O 0]τ;
[0025] Δ为4ΜΧ 1的残差矩阵，由M个InSAR观测误差和3Μ个模型误差组成：
[0026] A = [ n (X1) ..· η (χΜ) ε ! (X1) ε 2(Χι) ε 3(Xl) ......... E1(Xm) ε2(χΜ) ε 3 (χμ)]
[0027] Γ为（3Μ+Ν) X 1的待求参数矩阵，由M个地面观测点上的三维形变量和N个地下 块源的体积变化量组成：
[0028]
[0029] B为4ΜΧ (3Μ+Ν)的设计矩阵：
[0030]
[0031] 步骤7 :利用稀疏最小二乘算法求解联合模型，计算所有观测点上的三维形变量 和所有块源的体积变化量，即得到整个监测地区的三维地表形变量Cl 1 (Xl)、d2(Xl)和d3(Xl) 以及所有地下流体块源的体积变化量D v (yj)。
[0032] 所述地下块源的数量少于或等于所述地表观测点的数量。
[0033] 【以保证解算联合模型时不出现秩亏情况。】
[0034] 所述泊松比V取值为0. 25,地下流体的层数为1-3层。
[0035] 有益效果
[0036] 本发明提供了一种利用InSAR反演地下流体体积变化和三维地表形变的方法，1) 利用InSAR技术获取该地区的地理编码后的斜距向地表形变场的测量值；2)利用SAR卫星 成像几何建立InSAR斜距向形变测量值与三维地表形变之间的函数模型，并基于弹性半空 间理论建立三维地表形变与地下流体体积变化之间的函数模型；3)设定地下流体块源体 积均一变化，利用InSAR斜距向形变测量值非线性反演出地下流体块源的深度和厚度；4) 利用上述的两个函数模型建立针对所有地面观测点和地下流体块源的联合解算模型；最后 利用最小二乘方法对联合模型进行解算，得到所有地面观测点的三维形变量和所有地下流 体块源的体积变化量。该方法实现简单，不受区域的限制，且不依赖于任何其它大地测量数 据（如GPS)，对地下流体运动引起的三维地表形变而言是一种高精度、大范围、低成本和切 实可行的监测方法。不仅突破目前InSAR难以获得高精度三维地表形变场的技术瓶颈，积 极推动InSAR大地测量技术向实用化和市场化发展，而且还能同时获得地下流体的体积变 化结果，对研究地球内部的地球物理过程具有重要的科学价值和指导意义。
【附图说明】
[0037] 图1是本发明所述方法的流程图；
[0038] 图2是地下流体体积变化引起地表形变的示意图；
[0039] 图3是SAR卫星的成像几何图；
[0040] 图4是模拟的地下流体体积变化；
[0041] 图5是由模拟的地下流体体积变化造成的地表形变，其中，（a)为东西向形变；(b) 为南北向形变；（c)为垂直向形变；（d)为含噪的