InSAR斜距向形变测量值；单位：cm ;
[0042] 图6是本发明反演得到的地下流体体积变化；
[0043] 图7为应用本发明得到的三维地表形变与三维地表形变之间的差值；其中，（a)是 东西向三维地表形变，（b)是南北向三维地表形变，（c)是垂直向三维地表形变，（d)是东西 向反演和模拟的三维地表形变之间的差值，（e)是南北向反演和模拟的三维地表形变之间 的差值，（f)为是垂直向反演和模拟的三维地表形变之间的差值；单位：cm)。
【具体实施方式】
[0044] 下面将结合附图和实施例对本发明做进一步的说明。
[0045] 为了便于理解本发明，首先提供本发明的理论基础：
[0046] 如图2所示，根据弹性半空间理论，地下流体的体积变化会引起地表发生形变，并 且二者的关系满足：
[0047] (I1 (x) = I vGi (x, y) Dv (y) dy ⑴
[0048] 其中Cl1 (X)为地面观测点x的形变量，I = 1，2, 3分别代表东西、南北和垂直向上 的三个分量；DV (y)为地下半空间体积V中块源y的流体体积变化量！G1 (X，y)则为Green函 数，
[0049]
(2)
[0050] 其中V为泊松比，P1(X)和？七）分别为地面观测点X和块源y的三维空间位置，
为地下块源y到地面观测点X之间 的距离。公式（1)表明地表观测点的地表形变是地下半空间体积V内的所有块源贡献的总 和。
[0051] 而对于地表观测点X而言，根据SAR卫星成像几何（图3)，其上的InSAR斜距向形 变测量值I(X)可以写成：
[0052]
[0053] 其中S1(X)A2(X)和&〇〇分别为地面观测点X的东西、南北和垂直向地表形变在 InSAR斜距向上的投影系数，
[0054]
(4)
[0055] 其中Θ和α分别为雷达局部入射角和卫星飞行方向角（以北方向为起始顺时针 旋转）。
[0056] 如图1所示，一种利用InSAR反演地下流体体积变化和三维地表形变的方法，包括 以下步骤：
[0057] (1)利用D-InSAR或多时相InSAR技术获取待监测的受地下流体运动影响的地区 的地表在雷达视线方向（即斜距向）的形变测量值，并对其进行地理编码；
[0058] (2)利用SAR影像头文件中包含的雷达入射角、卫星飞行方向角，按照公式⑷计 算每个观测点上的投影系数S 1 (X)，S2(X)和S3(X)。假设共有M个地面观测点，那么对于地 面观测点X 1 (i = 1，2，...，M)而言，可以基于公式（3)构建InSAR斜距向形变测量值I (X1) 和三维地表形变之间的函数模型为
[0059]
[0060] 其中n (X1)为InSAR观测误差；
[0061] (3)由公式（1)和（3)可得每个观测点上的InSAR斜距向形变测量值和地下流体 的体积变化之间的函数关系：
[0062]
[0063] 假设地下流体中的所有块源的体积变化AV(y)是一致的，利用地下流体场的先 验信息确定泊松比（通常取0.25)和地下流体的层数（通常为1-3层），则公式（6)中的 未知参数只有三个：均一化体积变化A V，地下流体的深度h和厚度t。此时公式（6)为非 线性方程，利用M个地面观测点上的InSAR斜距向形变测量值，通过非线性方法（如遗传算 法）就可以反演出上述的三个未知参数。
[0064] (4)假设地下流体共包含有N个块源，可以基于公式（1)构建三维地表形变Cl1(X 1) (i = 1，2,. . .，M)与地下流体的体积变化量Dv(yj) (j = 1，2,. . .，N)之间的函数模型
[0066] 其中ε i (X1)代表模型与真实形变之间的残差，简称模型误差；'为单位块源的体 积，可由单位块源的尺寸和地下流体的厚度计算。
[0067] 那么对于M个地面观测点和N个地下块源而言，公式（5)和（7)可以合并构建如 下的联合模型：
[0068] Ω = Β Γ + Δ (8)
[0069] 其中Ω为4MX 1的观测矩阵，由M个地理编码后的InSAR斜距向形变测量值和3M 个虚拟观测量组成，
[0070]
[0071] Δ为4MX 1的残差矩阵，由M个InSAR观测误差和3M个模型误差组成，
[0073] Γ为（3M+N) X 1的待求参数矩阵，由M个地面观测点上的三维形变量和N个地下 块源的体积变化量组成，
[0074]
[0075] 而B为4MX (3M+N)的设计矩阵：
[0076] CN 105158760 A 说明书 6/7 页
[0077] (5)从公式（1)可以看出，当地下流体的深度h和厚度t确定时，三维地表形变和 地下流体体积变化之间是线性关系，而从公式（3)可以InSAR斜距向形变测量值和三维地 表形变之间也满足线性关系，因此基于公式（1)和（3)所构建的联合模型（即公式（8))仍 然是一个线性方程。当地面观测点的数量大于地下块源的数量时（即M>N)，就可以利用 最小二乘方法就可以求得未知参数Γ的最优解：
[0078] Il Ω-B Γ Il = min (8)
[0079] 其中I卜Il代表L2范准则。由于设计矩阵B是一个大型的稀疏矩阵，因此需要 利用稀疏最小二乘方法求解，最终得到所有地面观测点的三维地表形变量Cl 1 (X1)、d2 (X1)和 d3 (X1)以及所有地下流体块源的体积变化量Dv (yj)。
[0080] 在400 X 450的规则格网中模拟出地下流体的体积变化，格网尺寸为IOmX 10m，如 图4所示。其中地下流体的深度和厚度均取为100m，泊松比取为0. 25。则根据公式（1)，可 计算出该模拟的地下流体体积变化造成的三维地表形变，其中东西向、南北向和垂直向形 变分别如图5(a)-5(c)所示。然后通过公式（3)模拟出InSAR斜距向形变测量值，其中雷 达入射角和卫星飞行方向角采用升轨的AL0S/PALSAR卫星影像头文件中提供的参数。为了 让模拟数据具有真实性，将均值为零、标准偏差为2mm的高斯白噪声加入到InSAR斜距向形 变测量值，结果如图5(d)所示。由于该模拟数据是直接模拟在地理坐标系下，因此在这次 试验中无需地理编码步骤。
[0081] 通过本发明所提出的方法处理，就可以利用上述模拟的含噪InSAR斜距向形变测 量值同时反演出三维地表形变和地下流体的体积变化。为了保证公式（8)在解算时能够有 足够多的冗余观测，本次实验在解算时，将地下流体块源划分到40 X 45的规则格网中。图 6显示的就是本发明反演得到的地下流体体积变化。可以看出，虽然该反演结果在分辨率 上比原始模拟的地下流体体积变化降低了 100倍，但二者的空间变化趋势和幅度都非常一 致。图7(a)-7 (c)显示的分别是本发明反演得到的东西向、南北向和垂直向地表形变，总体 而言该三维形变与原始模拟的三维形变非常吻合。图7(d)-7(f)给出的则是反演和模拟的 三维形变场之间的差值。可以看出，东西向和垂直向上的差值主要呈现为高斯噪声，这主要 是由于这两个方向的地表形变对InSAR斜距测量值较为敏感；而南北向上的差值则主要表 现为格网式的噪声，这主要是因为南北向形变受模型的影响更为显著，因此对格网的降采 样处理更为敏感。为了定量验证本发明的效果，实施例中分别计算东西向、南北向和垂直向 形变的均方根误差，均不到1mm，明显小于InSAR斜距测量值中噪声的标准差，从而说明本 发明是可行的和可靠的。
[0082] 参考文献：
[0083] [l]Massonnet，D.，Rossi，Μ·，Carmona，C.，Adragna，F.，Peltzer，G.，Feigl，Κ· L. , Rabautej T. , 1993. The displacement field of the Landers earthquake mapped by radar interferometry. Nature 364 (6433)，138-142 ;
[0084] [2] Ferretti，A.，Prati，C.，Rocca，F.，2001. Permanent scatterers in SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing 39(1),8-20;
[0085] [3]Berardinoj P. , Fornaroj G. , Lanarij R. , Sansostij E. , 2002. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms.IEEE Transactions on Geoscience