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

摘要

本发明公开了一种利用InSAR反演地下流体体积变化和三维地表形变的方法，针对受地下流体运动影响的地区，首先利用InSAR技术获取该地区的地理编码后的斜距向地表形变场的测量值；随后利用InSAR斜距向形变测量值非线性反演出地下流体块源的深度和厚度；进而利用卫星成像几何和弹性半空间理论建立针对所有地面观测点和地下流体块源的联合解算模型；最后利用最小二乘方法解算出地面观测点的三维形变量和地下流体块源的体积变化量。突破了InSAR难以精确监测地下流体运动引起的三维形变的技术瓶颈，积极推动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\times G_1\left(x_i,y_1\right)& V_y\times G_1\left(x_i,y_2\right)& \mathrm{...}& V_y\times G_1\left(x_i,y_N\right)\\ V_y\times G_2\left(x_i,y_1\right)& V_y\times G_2\left(x_i,y_2\right)& \mathrm{...}& V_y\times G_2\left(x_i,y_N\right)\\ V_y\times G_3\left(x_i,y_1\right)& V_y\times G_3\left(x_i,y_2\right)& \mathrm{...}& V_y\times 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}{ϵ}_1\left(x_i\right)\\ {ϵ}_2\left(x_i\right)\\ {ϵ}_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\pi }\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(\alpha -3\pi /2)sin\theta \\ S_2\left(x\right)=-sin(\alpha -3\pi /2)sin\theta \\ S_3\left(x\right)=cos\theta \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\times G_1\left(x_i,y_1\right)& V_y\times G_1\left(x_i,y_2\right)& \mathrm{...}& V_y\times G_1\left(x_i,y_N\right)\\ V_y\times G_2\left(x_i,y_1\right)& V_y\times G_2\left(x_i,y_2\right)& \mathrm{...}& V_y\times G_2\left(x_i,y_N\right)\\ V_y\times G_3\left(x_i,y_1\right)& V_y\times G_3\left(x_i,y_2\right)& \mathrm{...}& V_y\times 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}{ϵ}_1\left(x_i\right)\\ {ϵ}_2\left(x_i\right)\\ {ϵ}_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斜距测量值较为敏感；而南北向上的差值则主要表现为格网式的噪声，这主要 是因为南北向形变受模型的影响更为显著，因此对格网的降采样处理更为敏感。为了定量验 证本发明的效果，实施例中分别计算东西向、南北向和垂直向形变的均方根误差，均不到1 mm，明显小于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.IEEE TransactionsonGeoscienceandRemoteSensing39(1),8-20；

[3]Berardino,P.,Fornaro,G.,Lanari,R.,Sansosti,E.,2002.Anewalgorithmforsurface deformationmonitoringbasedonsmallbaselinedifferentialSARinterferograms.IEEETransactions onGeoscienceandRemoteSensing40(11),2375-2383；

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

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

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

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

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

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

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

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