CN105158760A - 一种利用InSAR反演地下流体体积变化和三维地表形变的方法 - Google Patents

Abstract

本发明公开了一种利用InSAR反演地下流体体积变化和三维地表形变的方法，针对受地下流体运动影响的地区，首先利用InSAR技术获取该地区的地理编码后的斜距向地表形变场的测量值；随后利用InSAR斜距向形变测量值非线性反演出地下流体块源的深度和厚度；进而利用卫星成像几何和弹性半空间理论建立针对所有地面观测点和地下流体块源的联合解算模型；最后利用最小二乘方法解算出地面观测点的三维形变量和地下流体块源的体积变化量。突破了InSAR难以精确监测地下流体运动引起的三维形变的技术瓶颈，积极推动InSAR技术向实用化发展；而且还能同时获得地下流体体积变化结果，对研究地球内部的地球物理过程具有重要的科学价值。

Description

一种利用InSAR反演地下流体体积变化和三维地表形变的方法

技术领域

本发明属于基于遥感影像的大地测量领域，尤其涉及一种利用InSAR反演地下流体体积变化和三维地表形变的方法。

背景技术

差分InSAR(DifferentialInSAR，D-InSAR)技术是目前国际上在InSAR应用上最为成熟的技术，它最主要的目的就是监测地球表面厘米级甚至毫米级的形变。多时相InSAR(Multi-TemporalInSAR,MT-InSAR)技术则是近20多年来在D-InSAR技术的基础上发展起来的一种地表形变监测方法，如PS、SBAS和TCP等。通过对多景SAR影像的联合分析，MT-InSAR技术可以更好地抑制干涉图中时空失相关和大气噪声的影响，并能提供地表在SAR影像获取时刻的形变序列。然而，到目前为止，无论是D-InSAR还是MT-InSAR技术都只能利用单一平台、单一轨道的SAR数据，因此只能监测地表在雷达视线方向(即斜距向)上的一维形变结果。但是在现实中，地表形变是发生在三维空间框架下的，即所谓的三维形变场。因此InSAR的斜距向形变测量结果不能反映出真实的地表形变。

如何将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台站较少的地区不适用。

而在现实中，往往存在着一些缓慢的地表形变需要进行大范围、全方位和长期的监测，如地下水抽取、油气田开采和岩浆活动等地下流体运动引起的地表形变。这些地表形变已经成为影响区域经济和社会可持续发展的重要因素。例如，由于长期过度集中开采地下水资源，我国迄今为止已有70多个城市发生了不同程度的地面沉降，其中上海和天津市区的最大沉降已超过两米，严重制约了城市社会经济的发展，同时对人们的生活也构成了极大的威胁。然而根据上述分析，目前InSAR技术难以广泛地用于监测由地下流体运动引起的三维地表形变。

发明内容

本发明的目的在于，克服现有InSAR技术在监测地下流体运动引起的地表形变的不足和局限性，提供一种利用InSAR斜距向形变测量值监测三维地表形变，并同步反演出地下流体体积变化量的方法。

一种利用InSAR反演地下流体体积变化和三维地表形变的方法，包括以下步骤：

步骤1：利用D-InSAR或多时相InSAR技术获取待监测的受地下流体运动影响的地区的地表在雷达视线方向的形变测量值，并对其进行地理编码；

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

步骤2：根据SAR卫星成像几何，按照以下公式构建InSAR斜距向形变测量值I(x_i)和三维地表形变之间的函数模型：

I(x_i)＝[S₁(x_i)S₂(x_i)S₃(x_i)][d₁(x_i)d₂(x_i)d₃(x_i)]^T+η(x_i)

其中，x_i表示地表观测点，i＝1,2,...,M，共有M个地表观测点；

S₁(x_i)、S₂(x_i)及S₃(x_i)分别为地表观测点x_i的东西、南北和垂直向地表形变在InSAR斜距向上的投影系数，S₁(x_i)＝-cos(α-3π/2)sinθ，S₂(x_i)＝-sin(α-3π/2)sinθ，S₃(x_i)＝cosθ；θ和α分别为雷达局部入射角和卫星飞行方向角；

d_l(x_i)为地表观测点位置x_i的形变量，l＝1,2,3分别代表东西、南北和垂直向上的三个分量；

η(x_i)为InSAR观测误差；

【该值很小，通过最小二乘平差可以使其对InSAR斜距向形变测量值的影响最小；】

步骤3：根据弹性半空间理论建立每个观测点上的三维地表形变与地下流体的体积变化之间的函数模型：

$\left[\begin{array}{c}{d}_1\left(x_i\right)\\ d_2\left(x_i\right)\\ d_3\left(x_i\right)\end{array}\right]=\left[\begin{array}{cccc}{V}_y\×G_1\left(x_i,y_1\right)& V_y\×G_1\left(x_i,y_2\right)& \mathrm{...}& V_y\×G_1\left(x_i,y_N\right)\\ V_y\×G_2\left(x_i,y_1\right)& V_y\×G_2\left(x_i,y_2\right)& \mathrm{...}& V_y\×G_2\left(x_i,y_N\right)\\ V_y\×G_3\left(x_i,y_1\right)& V_y\×G_3\left(x_i,y_2\right)& \mathrm{...}& V_y\×G_3\left(x_i,y_N\right)\end{array}\right]\left[\begin{array}{c}{D}_V\left(y_1\right)\\ D_V\left(y_2\right)\\ .\\ .\\ .\\ D_V\left(y_N\right)\end{array}\right]-\left[\begin{array}{c}{\mathrm{\ϵ}}_1\left(x_i\right)\\ {\mathrm{\ϵ}}_2\left(x_i\right)\\ {\mathrm{\ϵ}}_3\left(x_i\right)\end{array}\right]$

其中，G_l(x,y)为格林函数，ν为泊松比，S为地下块源y到地表观测点x之间的距离，共有N个块源；P_l(x)和P_l(y)分别为地表观测点x和块源y的三维空间位置， $S=\sqrt{{\left(P_1\right(x)-P_1(y\left)\right)}^2+{\left(P_2\right(x)-P_2(y\left)\right)}^2+{\left(P_3\right(x)-P_3(y\left)\right)}^2}$ 为地下块源y到地面观测点x之间的距离；D_V(y)为地下半空间体积V中块源y的流体体积变化量；V_y表示单位块源的体积，由单位块源的尺寸和地下流体的厚度计算获得；ε_l(x_i)表示模型与真实形变之间的残差；

步骤4：利用SAR卫星成像几何将步骤3得到的模型转化为地表每个观测点上的InSAR斜距向形变测量值与地下流体的体积变化之间的函数关系：

I(x)＝∫_V(S₁G₁(x,y)+S₂G₂(x,y)+S₃G₃(x,y))ΔV(y)dy

步骤5：令地下流体中的所有块源的体积变化ΔV(y)相同，利用地下流体场的先验信息确定泊松比ν和地下流体的层数，利用M个地表观测点上的InSAR斜距向形变测量值，通过非线性方法反演出块源的体积变化ΔV(y)、地下流体的深度h和厚度t；

步骤6：构建以所有地表观测点的三维形变量和所有地下块源的体积变化量为未知参数的联合模型：

Ω＝ΒΓ+Δ

其中Ω为4M×1的观测矩阵，由M个地理编码后的InSAR斜距向形变测量值和3M个虚拟观测量组成：Ω＝[I(x₁)…I(x_M)000………000]^T；

Δ为4M×1的残差矩阵，由M个InSAR观测误差和3M个模型误差组成：

Δ＝[η(x₁)…η(x_M)ε₁(x₁)ε₂(x₁)ε₃(x₁)………ε₁(x_M)ε₂(x_M)ε₃(x_M)]^T

Γ为(3M+N)×1的待求参数矩阵，由M个地面观测点上的三维形变量和N个地下块源的体积变化量组成：

Γ＝[d₁(x₁)d₂(x₁)d₃(x₁)………d₁(x_M)d₂(x_M)d₃(x_M)D_V(y₁)…D_V(y_N)]^T

Β为4M×(3M+N)的设计矩阵：

步骤7：利用稀疏最小二乘算法求解联合模型，计算所有观测点上的三维形变量和所有块源的体积变化量，即得到整个监测地区的三维地表形变量d₁(x_i)、d₂(x_i)和d₃(x_i)以及所有地下流体块源的体积变化量D_V(y_j)。

所述地下块源的数量少于或等于所述地表观测点的数量。

【以保证解算联合模型时不出现秩亏情况。】

所述泊松比ν取值为0.25，地下流体的层数为1-3层。

有益效果

本发明提供了一种利用InSAR反演地下流体体积变化和三维地表形变的方法，1)利用InSAR技术获取该地区的地理编码后的斜距向地表形变场的测量值；2)利用SAR卫星成像几何建立InSAR斜距向形变测量值与三维地表形变之间的函数模型，并基于弹性半空间理论建立三维地表形变与地下流体体积变化之间的函数模型；3)设定地下流体块源体积均一变化，利用InSAR斜距向形变测量值非线性反演出地下流体块源的深度和厚度；4)利用上述的两个函数模型建立针对所有地面观测点和地下流体块源的联合解算模型；最后利用最小二乘方法对联合模型进行解算，得到所有地面观测点的三维形变量和所有地下流体块源的体积变化量。该方法实现简单，不受区域的限制，且不依赖于任何其它大地测量数据(如GPS)，对地下流体运动引起的三维地表形变而言是一种高精度、大范围、低成本和切实可行的监测方法。不仅突破目前InSAR难以获得高精度三维地表形变场的技术瓶颈，积极推动InSAR大地测量技术向实用化和市场化发展，而且还能同时获得地下流体的体积变化结果，对研究地球内部的地球物理过程具有重要的科学价值和指导意义。

附图说明

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

图2是地下流体体积变化引起地表形变的示意图；

图3是SAR卫星的成像几何图；

图4是模拟的地下流体体积变化；

图5是由模拟的地下流体体积变化造成的地表形变，其中，(a)为东西向形变；(b)为南北向形变；(c)为垂直向形变；(d)为含噪的InSAR斜距向形变测量值；单位：cm；

图6是本发明反演得到的地下流体体积变化；

图7为应用本发明得到的三维地表形变与三维地表形变之间的差值；其中，(a)是东西向三维地表形变，(b)是南北向三维地表形变，(c)是垂直向三维地表形变，(d)是东西向反演和模拟的三维地表形变之间的差值，(e)是南北向反演和模拟的三维地表形变之间的差值，(f)为是垂直向反演和模拟的三维地表形变之间的差值；单位：cm)。

具体实施方式

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

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

如图2所示，根据弹性半空间理论，地下流体的体积变化会引起地表发生形变，并且二者的关系满足：

d_l(x)＝∫_VG_l(x,y)D_V(y)dy(1)

其中d_l(x)为地面观测点x的形变量，l＝1,2,3分别代表东西、南北和垂直向上的三个分量；D_V(y)为地下半空间体积V中块源y的流体体积变化量；G_l(x,y)则为Green函数，

$G_l(x,y)=\frac{(v+1)}{3\mathrm{\π}}\frac{\left(P_l\right(x)-P_l(y\left)\right)}{S^3}---\left(2\right)$

其中ν为泊松比，P_l(x)和P_l(y)分别为地面观测点x和块源y的三维空间位置，S＝ $\sqrt{{\left(P_1\right(x)-P_1(y\left)\right)}^2+{\left(P_2\right(x)-P_2(y\left)\right)}^2+{\left(P_3\right(x)-P_3(y\left)\right)}^2}$ 为地下块源y到地面观测点x之间的距离。公式(1)表明地表观测点的地表形变是地下半空间体积V内的所有块源贡献的总和。

而对于地表观测点x而言，根据SAR卫星成像几何(图3)，其上的InSAR斜距向形变测量值I(x)可以写成：

I(x)＝[S₁(x)S₂(x)S₃(x)][d₁(x)d₂(x)d₃(x)]^T+η(x)(3)

其中S₁(x)，S₂(x)和S₃(x)分别为地面观测点x的东西、南北和垂直向地表形变在InSAR斜距向上的投影系数，

$\left\{\begin{array}{c}{S}_1\left(x\right)=-cos(\mathrm{\α}-3\mathrm{\π}/2)sin\mathrm{\θ}\\ S_2\left(x\right)=-sin(\mathrm{\α}-3\mathrm{\π}/2)sin\mathrm{\θ}\\ S_3\left(x\right)=cos\mathrm{\θ}\end{array}\right.---\left(4\right)$

其中θ和α分别为雷达局部入射角和卫星飞行方向角(以北方向为起始顺时针旋转)。

如图1所示，一种利用InSAR反演地下流体体积变化和三维地表形变的方法，包括以下步骤：

(1)利用D-InSAR或多时相InSAR技术获取待监测的受地下流体运动影响的地区的地表在雷达视线方向(即斜距向)的形变测量值，并对其进行地理编码；

(2)利用SAR影像头文件中包含的雷达入射角、卫星飞行方向角，按照公式(4)计算每个观测点上的投影系数S₁(x)，S₂(x)和S₃(x)。假设共有M个地面观测点，那么对于地面观测点x_i(i＝1,2,...,M)而言，可以基于公式(3)构建InSAR斜距向形变测量值I(x_i)和三维地表形变之间的函数模型为

I(x_i)＝[S₁(x_i)S₂(x_i)S₃(x_i)][d₁(x_i)d₂(x_i)d₃(x_i)]^T+η(x_i)(5)

其中η(x_i)为InSAR观测误差；

(3)由公式(1)和(3)可得每个观测点上的InSAR斜距向形变测量值和地下流体的体积变化之间的函数关系：

I(x)＝∫_V(S₁G₁(x,y)+S₂G₂(x,y)+S₃G₃(x,y))ΔV(y)dy(6)

假设地下流体中的所有块源的体积变化ΔV(y)是一致的，利用地下流体场的先验信息确定泊松比(通常取0.25)和地下流体的层数(通常为1-3层)，则公式(6)中的未知参数只有三个：均一化体积变化ΔV，地下流体的深度h和厚度t。此时公式(6)为非线性方程，利用M个地面观测点上的InSAR斜距向形变测量值，通过非线性方法(如遗传算法)就可以反演出上述的三个未知参数。

(4)假设地下流体共包含有N个块源，可以基于公式(1)构建三维地表形变d_l(x_i)(i＝1,2,...,M)与地下流体的体积变化量D_V(y_j)(j＝1,2,...,N)之间的函数模型

$\left[\begin{array}{c}{d}_1\left(x_i\right)\\ d_2\left(x_i\right)\\ d_3\left(x_i\right)\end{array}\right]=\left[\begin{array}{cccc}{V}_y\×G_1\left(x_i,y_1\right)& V_y\×G_1\left(x_i,y_2\right)& \mathrm{...}& V_y\×G_1\left(x_i,y_N\right)\\ V_y\×G_2\left(x_i,y_1\right)& V_y\×G_2\left(x_i,y_2\right)& \mathrm{...}& V_y\×G_2\left(x_i,y_N\right)\\ V_y\×G_3\left(x_i,y_1\right)& V_y\×G_3\left(x_i,y_2\right)& \mathrm{...}& V_y\×G_3\left(x_i,y_N\right)\end{array}\right]\left[\begin{array}{c}{D}_V\left(y_1\right)\\ D_V\left(y_2\right)\\ .\\ .\\ .\\ D_V\left(y_N\right)\end{array}\right]-\left[\begin{array}{c}{\mathrm{\ϵ}}_1\left(x_i\right)\\ {\mathrm{\ϵ}}_2\left(x_i\right)\\ {\mathrm{\ϵ}}_3\left(x_i\right)\end{array}\right]---\left(7\right)$

其中ε_l(x_i)代表模型与真实形变之间的残差，简称模型误差；V_y为单位块源的体积，可由单位块源的尺寸和地下流体的厚度计算。

那么对于M个地面观测点和N个地下块源而言，公式(5)和(7)可以合并构建如下的联合模型：

Ω＝ΒΓ+Δ(8)

其中Ω为4M×1的观测矩阵，由M个地理编码后的InSAR斜距向形变测量值和3M个虚拟观测量组成，

Ω＝[I(x₁)…I(x_M)000………000]^T

Δ为4M×1的残差矩阵，由M个InSAR观测误差和3M个模型误差组成，

Δ＝[η(x₁)…η(x_M)ε₁(x₁)ε₂(x₁)ε₃(x₁)………ε₁(x_M)ε₂(x_M)ε₃(x_M)]^T

Γ为(3M+N)×1的待求参数矩阵，由M个地面观测点上的三维形变量和N个地下块源的体积变化量组成，

Γ＝[d₁(x₁)d₂(x₁)d₃(x₁)………d₁(x_M)d₂(x_M)d₃(x_M)D_V(y₁)…D_V(y_N)]^T

而Β为4M×(3M+N)的设计矩阵：

(5)从公式(1)可以看出，当地下流体的深度h和厚度t确定时，三维地表形变和地下流体体积变化之间是线性关系，而从公式(3)可以InSAR斜距向形变测量值和三维地表形变之间也满足线性关系，因此基于公式(1)和(3)所构建的联合模型(即公式(8))仍然是一个线性方程。当地面观测点的数量大于地下块源的数量时(即M＞N)，就可以利用最小二乘方法就可以求得未知参数Γ的最优解：

‖Ω-ΒΓ‖＝min(8)

其中‖·‖代表L2范准则。由于设计矩阵Β是一个大型的稀疏矩阵，因此需要利用稀疏最小二乘方法求解，最终得到所有地面观测点的三维地表形变量d₁(x_i)、d₂(x_i)和d₃(x_i)以及所有地下流体块源的体积变化量D_V(y_j)。

在400×450的规则格网中模拟出地下流体的体积变化，格网尺寸为10m×10m，如图4所示。其中地下流体的深度和厚度均取为100m，泊松比取为0.25。则根据公式(1)，可计算出该模拟的地下流体体积变化造成的三维地表形变，其中东西向、南北向和垂直向形变分别如图5(a)-5(c)所示。然后通过公式(3)模拟出InSAR斜距向形变测量值，其中雷达入射角和卫星飞行方向角采用升轨的ALOS/PALSAR卫星影像头文件中提供的参数。为了让模拟数据具有真实性，将均值为零、标准偏差为2mm的高斯白噪声加入到InSAR斜距向形变测量值，结果如图5(d)所示。由于该模拟数据是直接模拟在地理坐标系下，因此在这次试验中无需地理编码步骤。

通过本发明所提出的方法处理，就可以利用上述模拟的含噪InSAR斜距向形变测量值同时反演出三维地表形变和地下流体的体积变化。为了保证公式(8)在解算时能够有足够多的冗余观测，本次实验在解算时，将地下流体块源划分到40×45的规则格网中。图6显示的就是本发明反演得到的地下流体体积变化。可以看出，虽然该反演结果在分辨率上比原始模拟的地下流体体积变化降低了100倍，但二者的空间变化趋势和幅度都非常一致。图7(a)-7(c)显示的分别是本发明反演得到的东西向、南北向和垂直向地表形变，总体而言该三维形变与原始模拟的三维形变非常吻合。图7(d)-7(f)给出的则是反演和模拟的三维形变场之间的差值。可以看出，东西向和垂直向上的差值主要呈现为高斯噪声，这主要是由于这两个方向的地表形变对InSAR斜距测量值较为敏感；而南北向上的差值则主要表现为格网式的噪声，这主要是因为南北向形变受模型的影响更为显著，因此对格网的降采样处理更为敏感。为了定量验证本发明的效果，实施例中分别计算东西向、南北向和垂直向形变的均方根误差，均不到1mm，明显小于InSAR斜距测量值中噪声的标准差，从而说明本发明是可行的和可靠的。

参考文献：

[1]Massonnet,D.,Rossi,M.,Carmona,C.,Adragna,F.,Peltzer,G.,Feigl,K.L.,Rabaute,T.,1993.ThedisplacementfieldoftheLandersearthquakemappedbyradarinterferometry.Nature364(6433),138-142；

[2]Ferretti,A.,Prati,C.,Rocca,F.,2001.PermanentscatterersinSARinterferometry.IEEETransactionsonGeoscienceandRemoteSensing39(1),8-20；

[3]Berardino,P.,Fornaro,G.,Lanari,R.,Sansosti,E.,2002.AnewalgorithmforsurfacedeformationmonitoringbasedonsmallbaselinedifferentialSARinterferograms.IEEETransactionsonGeoscienceandRemoteSensing40(11),2375-2383；

[4]Zhang,L.,Lu,Z.,Ding,X.L.,Jung,H.S.,Feng,G.C.,Lee,C.W.,2012.MappinggroundsurfacedeformationusingtemporarilycoherentpointSARinterferometry:ApplicationtoLosAngelesBasin.RemoteSensingofEnvironment117,429-439；

[5]Rocca,F.,2003.3Dmotionrecoveryfrommulti-angleand/orleftrightinterferometry.In:ProceedingsofthethirdInternationalWorkshoponERSSAR；

[6]Wright,T.J.,Parsons,B.E.,Lu,Z.,2004.TowardmappingsurfacedeformationinthreedimensionsusingInSAR.GeophysicalResearchLetters31(1),doi:01610.01029/02003gl018827。

[7]Gray,L.,2011.UsingmultipleRADARSATInSARpairstoestimateafullthree-dimensionalsolutionforglacialicemovement.GeophysicalResearchLetters38(5),doi:10.1029/2010gl046484；

[8]Fialko,Y.,Simons,M.,Agnew,D.,2001.Thecomplete(3-D)surfacedisplacementfieldintheepicentralareaofthe1999M(w)7.1HectorMineearthquake,California,fromspacegeodeticobservations.GeophysicalResearchLetters28(16),3063-3066；

[9]Fialko,Y.,Sandwell,D.,Simons,M.,Rosen,P.,2005.Three-dimensionaldeformationcausedbytheBam,Iran,earthquakeandtheoriginofshallowslipdeficit.Nature435,295-299；

[10]Hu,J.,Li,Z.W.,Ding,X.L.,Zhu,J.J.,Zhang,L.,Sun,Q.,2012.3Dcoseismicdisplacementof2010Darfield,NewZealandearthquakeestimatedfrommulti-apertureInSARandD-InSARmeasurements.JournalofGeodesy86,1029-1041；

[11]Samsonov,S.,Tiampo,K.,Rundle,J.,Li,Z.H.,2007.ApplicationofDInSAR-GPSoptimizationforderivationoffine-scalesurfacemotionmapsofsouthernCalifornia.IEEETransactionsonGeoscienceandRemoteSensing45(2),512-521；

Claims (3)

1.一种利用InSAR反演地下流体体积变化和三维地表形变的方法，其特征在于，包括以下步骤：

步骤1：利用D-InSAR或多时相InSAR技术获取待监测的受地下流体运动影响的地区的地表在雷达视线方向的形变测量值，并对其进行地理编码；

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

步骤2：根据SAR卫星成像几何，按照以下公式构建InSAR斜距向形变测量值I(x_i)和三维地表形变之间的函数模型为：

I(x_i)＝[S₁(x_i)S₂(x_i)S₃(x_i)][d₁(x_i)d₂(x_i)d₃(x_i)]^T+η(x_i)

其中，x_i表示地表观测点，i＝1,2,...,M，共有M个地表观测点；

S₁(x_i)，S₂(x_i)及S₃(x_i)分别为地表观测点x_i的东西、南北和垂直向地表形变在InSAR斜距向上的投影系数，S₁(x_i)＝-cos(α-3π/2)sinθ，S₂(x_i)＝-sin(α-3π/2)sinθ，S₃(x_i)＝cosθ；θ和α分别为雷达局部入射角和卫星飞行方向角；

d_l(x_i)为地表观测点位置x_i的形变量，l＝1,2,3分别代表东西、南北和垂直向上的三个分量；

η(x_i)为InSAR观测误差；

步骤3：根据弹性半空间理论建立每个观测点上的三维地表形变与地下流体的体积变化之间的函数模型：

$\left[\begin{array}{c}{d}_1\left(x_i\right)\\ d_2\left(x_i\right)\\ d_3\left(x_i\right)\end{array}\right]=\left[\begin{array}{cccc}{V}_y\×G_1\left(x_i,y_1\right)& V_y\×G_1\left(x_i,y_2\right)& \mathrm{...}& V_y\×G_1\left(x_i,y_N\right)\\ V_y\×G_2\left(x_i,y_1\right)& V_y\×G_2\left(x_i,y_2\right)& \mathrm{...}& V_y\×G_2\left(x_i,y_N\right)\\ V_y\×G_3\left(x_i,y_1\right)& V_y\×G_3\left(x_i,y_2\right)& \mathrm{...}& V_y\×G_3\left(x_i,y_N\right)\end{array}\right]\left[\begin{array}{c}{D}_V\left(y_1\right)\\ D_V\left(y_2\right)\\ .\\ .\\ .\\ D_V\left(y_N\right)\end{array}\right]-\left[\begin{array}{c}{\mathrm{\ϵ}}_1\left(x_i\right)\\ {\mathrm{\ϵ}}_2\left(x_i\right)\\ {\mathrm{\ϵ}}_3\left(x_i\right)\end{array}\right]$

其中，G_l(x,y)为格林函数，v为泊松比，S为地下块源y到地表观测点x之间的距离，共有N个块源；P_l(x)和P_l(y)分别为地表观测点x和块源y的三维空间位置， $S=\sqrt{{\left(P_1\right(x)-P_1(y\left)\right)}^2+{\left(P_2\right(x)-P_2(y\left)\right)}^2+{\left(P_3\right(x)-P_3(y\left)\right)}^2}$ 为地下块源y到地面观测点x之间的距离；D_V(y)为地下半空间体积V中块源y的流体体积变化量；V_y表示单位块源的体积，由单位块源的尺寸和地下流体的厚度计算获得；ε_l(x_i)表示模型与真实形变之间的残差；

步骤4：利用SAR卫星成像几何将步骤3得到的模型转化为地表每个观测点上的InSAR斜距向形变测量值与地下流体的体积变化之间的函数关系：

I(x)＝∫_V(S₁G₁(x,y)+S₂G₂(x,y)+S₃G₃(x,y))ΔV(y)dy

步骤5：令地下流体中的所有块源的体积变化ΔV(y)相同，利用地下流体场的先验信息确定泊松比v和地下流体的层数，利用M个地表观测点上的InSAR斜距向形变测量值，通过非线性方法反演出块源的体积变化ΔV(y)、地下流体的深度h和厚度t；

步骤6：构建以所有地表观测点的三维形变量和所有地下块源的体积变化量为未知参数的联合模型：

Ω＝ΒΓ+Δ

其中，Ω为4M×1的观测矩阵，由M个地理编码后的InSAR斜距向形变测量值和3M个虚拟观测量组成：Ω＝[I(x₁)…I(x_M)000………000]^T；

Δ为4M×1的残差矩阵，由M个InSAR观测误差和3M个模型误差组成：

Δ＝[η(x₁)…η(x_M)ε₁(x₁)ε₂(x₁)ε₃(x₁)………ε₁(x_M)ε₂(x_M)ζ₃(x_M)]^T

Γ为(3M+N)×1的待求参数矩阵，由M个地面观测点上的三维形变量和N个地下块源的体积变化量组成：

Γ＝[d₁(x₁)d₂(x₁)d₃(x₁)………d₁(x_M)d₂(x_M)d₃(x_M)D_V(y₁)…D_V(y_N)]^T

Β为4M×(3M+N)的设计矩阵：

步骤7：利用稀疏最小二乘算法求解联合模型，计算所有观测点上的三维形变量和所有块源的体积变化量，即得到整个监测地区的三维地表形变量d₁(x_i)、d₂(x_i)和d₃(x_i)以及所有地下流体块源的体积变化量D_V(y_j)。

2.根据权利要求1所述的一种利用InSAR反演地下流体体积变化和三维地表形变的方法，其特征在于，所述地下块源的数量少于或等于所述地表观测点的数量。

3.根据权利要求1或2所述的一种利用InSAR反演地下流体体积变化和三维地表形变的方法，其特征在于，所述泊松比v取值为0.25，地下流体的层数为1-3层。