EN中文

世界核地质科学 >

2025 , Vol. 42 >Issue 3: 565 - 581

DOI: https://doi.org/10.3969/j.issn.1672-0636.2025.03.008

研究论文

可控源音频大地电磁法和微动谱比法的横向约束联合反演

  • 王崇龙 , 1, 2 ,
  • 彭淼 1
展开
  • 1 中国地质大学（北京） 地球物理与信息技术学院，北京 100083
  • 2 核工业二四○研究所，辽宁 沈阳 110000

王崇龙，男，1998年生，助理工程师，主要从事铀矿勘探的物探工作。E-mail：

收稿日期: 2025-02-26

  修回日期: 2025-04-11

  网络出版日期: 2025-11-06

基金资助

国家自然科学基金项目(42174086)

收起

Laterally constrained joint inversion of CSAMT and microtremor spectrum method

  • WANG Chonglong , 1, 2 ,
  • PENG Miao 1
Expand
  • 1 China University of Geosciences (Beijing),School of Geophysics and Information Technology，Beijing 100083，China
  • 2 Research Institute No.240，CNNC，Shenyang 110000，China

WANG Chonglong，male，born in 1998，assistant engineer，focusing on geophysical exploration for uranium mines. E-mail：

Received date: 2025-02-26

  Revised date: 2025-04-11

  Online published: 2025-11-06

Supported by

National Natural Science Foundation of China(42174086)

Fold

摘要

为使相邻测点间的横向物性和层参数分布更加光滑连续，降低单一地球物理反演方法的局限性，开展了可控源音频大地电磁法（CSAMT）和微动谱比法的拟二维横向约束联合反演研究。微动数据采用谱比法进行数值模拟，联合CSAMT基于有限内存拟牛顿法（L-BFGS）反演算法，引入横向约束理论，加入交叉梯度函数实现两种不同物性参数的相互耦合，开发出一套拟二维横向约束联合反演算法，并通过两组理论模型验证了算法的准确性和有效性。同时，利用反演算法对北京延庆地区的实测数据进行反演，结果表明：电阻率和横波速度的突变界面形态具有较好的对应关系，证明横向约束联合反演算法具备一定实用价值。

关键词： 可控源音频大地电磁法; 微动; 横向约束; 联合反演

本文引用格式

王崇龙 , 彭淼 . 可控源音频大地电磁法和微动谱比法的横向约束联合反演[J]. 世界核地质科学, 2025 , 42(3) : 565 -581 . DOI: 10.3969/j.issn.1672-0636.2025.03.008

Abstract

In order to make the distribution of lateral physical properties and layer parameters between adjacent measuring points smoother and more continuous and reduce the limitations of a single geophysical inversion method，a pseudo-two-dimensional lateral constrained joint inversion study of controlled source audio frequency magnetotelluric method (CSAMT) and micro-motion spectrum ratio method was carried out. The microtremor data is numerically simulated using the spectral ratio method，combined with the CSAMT-based limited memory BFGS (L-BFGS) inversion algorithm，introducing the lateral constraint theory，and adding the cross-gradient function to achieve the mutual coupling of two different physical parameters. A set of quasi-two-dimensional lateral constraint joint inversion algorithms was developed，and the accuracy and effectiveness of the algorithms were verified through two sets of theoretical models. Meanwhile，the inversion algorithm is used to invert the measured data in Yanqing，Beijing. The results show that there is a good correspondence between the abrupt interface morphology of resistivity and shear wave velocity，which proves the practical value of the laterally constrained joint inversion algorithm.

Key words： controlled source audio magnetotelluric method; microtremor; laterally constraint; joint inversion

在传统单点一维反演中，反演的物性参数会出现不连续现象，这使得反演结果中出现很多变异异常，给后续解释工作增加了难度。近些年的研究结果表明，横向约束方法可以有效压制上述突变异常，提高反演结果的横向连续性，进一步提升反演解释的精度[1]。蔡晶等[2]以横向约束反演理论作为基础，引入了加权因子的概念，提出加权横向约束反演理论, 并应用于频率域航空电磁数据的拟二维反演中。王若等[3]第一次在可控源音频大地电磁法的反演中加入了横向约束反演，并用最小二乘法对实测数据进行了反演，反演效果明显提高。殷长春等[4]采用阻尼最小二乘法，将横向约束反演应用于时间域航空电磁法拟二维反演中。
电法勘探和微动勘探联合反演属于不同物性参数之间的联合反演，而对于不同物性参数的联合反演，需要建立两种形式的约束，即岩石物性约束和结构耦合约束。岩石物性约束是对岩石物性的耦合，这种约束依赖于不同物性之间的经验关系式，而经验关系式的准确程度又会影响联合反演的精度，所以基于岩石物性的约束仍具有局限性。结构耦合约束最常用的方法就是交叉梯度约束。这种约束方法首先需要假设在空间上不同物性的分布具有相似性，不对其进行梯度标准化，这样不管梯度的规模是大还是小，都能检测到变化的差异。这个假设是成立的，因为在地下介质中，如果一种物性参数发生变化，其他的物性参数往往也会随之变化。Dobroka等[5]将最大频率值方法应用于地震与电法的联合反演中，即使实测数据中存在噪声，也能以足够的精度得到反演结果。杨振武等[6]利用大地电磁和地震的联合反演改善了反演结果的稳定性。过仲阳等[7]在大地电磁和地震联合反演中运用遗传算法，提高了反演结果的精度，具有较强的抗干扰能力，但时间成本较高。杨辉等[8]在仿真退火算法的基础上，实现了一维及二维起伏地形条件下大地电磁与地震资料多参量信息的联合反演。Heincke等[9]实现重力、地震层析成像和大地电磁的联合反演，可以得到将单方法无法单独恢复出的地下结构。
关于电阻率和速度模型的联合反演中目前已经取得诸多成果，但关于可控源音频大地电磁法（CSAMT）与微动谱比法的联合反演研究较少。两种方法的应用场景均适用于干扰较大的地区，同时CSAMT勘查效率高，探测深度大[10]，又具有较高的纵向分辨率以及不受高阻屏蔽的影响等优点；微动谱比法则有着较好的横向分辨率等优点，对层状介质恢复效果较好。因此考虑综合两种方法的特点，并根据各自方法的优势进行互补，完成两种地球物理方法的联合反演算法，提高反演结果准确性。
在前人提出的横向约束和交叉梯度的基础上，建立基于交叉梯度约束的CSAMT和微动数据横向约束反演的目标函数，实现电阻率和地震波速度两种物性参数之间相互约束的联合反演算法，通过模型算例的测试，讨论了横向约束联合反演算法的有效性和适用性。

1 CSAMT和微动拟二维横向约束反演方法

1.1 CSAMT拟二维横向约束反演方法

目前地球物理反演中广泛采用正则化反演，通过对目标函数求解来获得反演结果。创建如下的正则化目标函数：
φ m = φ d + λ φ m
式（1）中：
φ d = d o b s - f m T C d - 1 d o b s - f m
式（2）中： φ d—数据拟合差项，表示反演结果与观测数据的拟合程度； φ m—模型圆滑项； λ—权重因子； d o b s—观测数据向量，向量个数与观测频点数与测点数具有对应关系； f m—正演响应（所用正演计算，均采用CSAMT一维正演的形式，即某一测点处的正演结果仅与此测点下电阻率与层厚度分布有关，与周围测点电阻率与层厚度分布无关）。
对于频点数为 n，测点数为 a的所有测点的观测数据可以表示为：
d o b s = d 1 , d 2 , … , d a a n × 1 T  
式（3）中： d i—第i个测点处的观测数据； m—模型参数，每一测点下地下为 b层的模型由每一测点下的电阻率与层厚度组成，可以表示为：
m i = l n ρ 1 , l n ρ 2 , … , l n ρ b , l n h 1 , l n h 2 , … , l n h b - 1
测点数为 a的所有测点的模型参数可以表示为：
m = m 1 , m 2 , … , m a a * 2 b - 1 × 1 T
式（5）中： C d - 1—数据加权矩阵，具有如下形式：
C d - 1 = 1 σ 1           ⋱           1 σ j           ⋱           1 σ m n a × n a
式（6）中： σ j—第 j个数据的误差。如果去掉数据加权矩阵，则数据目标函数会变为最小二乘的形式。数据加权矩阵的作用有两个：一个是可以控制数据权重，使观测误差小的数据的权重升高或者使观测误差大的数据的权重降低；另一个是对数量级相差较大的两种数据进行归一化。
φ m = R m T R m
式（7）中： R—模型约束矩阵。
采用的是模型横向约束（Laterally Constrained Inversion，简称LCI）矩阵。LCI最早由奥胡斯大学的Esben Auken学者提出，并将其应用于电阻率法反演中[11-12]。横向约束算法的作用是传递相邻测点间的物性参数，用相对较好的反演结果来约束相对较差的反演结果，进而不断提升相邻测点之间物性参数的相似性，使反演结果变得更加连续。以三个测点的测线为例，横向约束理论可以表示为图1的形式：
图1 横向约束示意图[13]

Fig. 1 Laterally constrained scheme[13]

图1中， ρ i , j、 h i , j分别表示测点 i，第 j层的电阻率与层厚度。横向约束是建立在沿一条测线方向上的相邻测点间的约束，即每一个测点下的每一层的层厚度（数量为 b - 1）与每一层的电阻率（数量为 b）受相邻两个或一个（对于内部测点为两个，对于外部测点为一个）测点的约束，每一个测点下的每一层底界面的埋藏深度受相邻测点对应深度的约束。
用数学的形式表述这种关系为：
R p = 1 0 ⋯ 0 - 1 0 ⋯ 0 0 0 0 1 0 ⋯ 0 - 1 0 ⋯ 0 0 ⋮         ⋮       ⋮ 0 0 0 ⋯ 0 1 0 ⋯ 0 - 1 a - 1 * 2 b - 1 × a * 2 b - 1
R h = ⋯ 0 0 1 0 0 ⋯ 0       ⋯ 0 0 t k , 1 h k , 2 t k , 2 h k , 2 0 ⋯ 0         ⋮ ⋮ ⋮ ⋮     ⋮       ⋯ 0 0 t k , 1 h k , n t k , 2 h k , n t k , 3 h k , n ⋯ t k , n h k , n           ⋯ 0 ⋯ 0 - 1 0 ⋯ 0 ⋯     ⋯ 0 ⋯ 0 - t k + 1,1 h k + 1,2 - t k + 1,2 h k + 1,2 ⋯ 0 ⋯       ⋮   ⋮ ⋮     ⋮       ⋯ 0 ⋯ 0 - t k + 1,1 h k + 1 , n - t k + 1,2 h k + 1 , n ⋯ - t k + 1 , n h k + 1 , n ⋯ a - 1 * b - 1 × a * 2 b - 1
式（9）中： R p—电阻率与层厚度的横向约束矩阵，矩阵的列数等于反演参数总数（一维反演参数个数乘以测点数），矩阵的行数等于横向约束连接数（一维反演参数个数乘以测点数减一）； R h—层深度的横向约束矩阵，矩阵的列数等于反演参数总数，矩阵的行数等于深度横向约束连接数（一维反演深度数乘以测点数减一）； t i , j—表示第 i个测点第 j层下界面的深度。
采用的是正则化反演，模型约束项同时受电阻率、层厚度的约束与深度约束的影响，横向约束项和模型约束项的权重相差不大，考虑到两者的不同权重问题，拟二维模型横向约束矩阵 R直接表示为两个约束矩阵的相加：
[ R ] = [ R p ] + [ R h ]
对正则化目标函数求导可得：
∇ φ m = ∇ φ d + λ ∇ φ m
∇ φ d = - 2 J T C d - 1 d - f m
∇ φ m = 2 R T R m
式（12）中： J—灵敏度矩阵或雅可比矩阵。
L-BFGS算法为有限内存拟牛顿法，该方法计算建立在牛顿法的基础上，无需直接求解海森矩阵，具有高效快速的特点。反演算法均是基于L-BFGS算法开发的。下面对基本原理进行简要介绍。
首先对目标函数进行泰勒展开：
φ m k + 1 = φ m k + φ ' m k ∙ m k + 1 - m k + 1 2 φ ' ' m k ∙ m k + 1 - m k 2 + O m k + 1 - m k 2
其次忽略二阶及以上的高阶项，使目标函数趋于零，可以得到：
φ ' m k ∙ m k + 1 - m k = 0
式（15）对 m k求导，可以得到目标函数在二阶泰勒展开下，下降最快的方向：
φ ' ' m k ∙ m k + 1 - m k - φ ' m k = 0
m k + 1 = m k - φ ' m k φ ' ' m k  
最后表示为一种更为常见的形式为：
m k + 1 = m k - H k - 1 ∇ φ m k  
式（17）~（18）中： H k—海森矩阵，是目标函数的二阶导数； ∇ φ m k—目标函数的梯度。对于大规模数据与大量模型参数而言，求解海森矩阵是十分困难的。为此，一些学者提出了L-BFGS法，采用前几次迭代的模型参数、目标函数梯度，来贴近当前迭代所需的海森矩阵的逆，这样做大幅减少了对内存的要求，提高了计算效率。
海森矩阵逼近算法为：
s i = m i + 1 - m i
y i = ∇ φ m i + 1 - ∇ φ m i
ρ k = 1 y k T s k
$\begin{aligned} V_{k}= & \left(V_{k-1}^{T}, \cdots, V_{k-n}^{T}\right) H_{k}^{0}\left(V_{k-n}, \cdots, V_{k-1}\right)+\rho_{k-n}\left(V_{k-1}^{T}, \cdots, V_{k-n+1}^{T}\right) s_{k-n} s_{k-n}^{T}\left(V_{k-n+1}, \cdots, V_{k-1}\right) \\ & +\rho_{k-n+1}\left(V_{k-1}^{T}, \cdots, V_{k-n+2}^{T}\right) s_{k-n+1} s_{k-n+1}^{T}\left(V_{k-n+2}, \cdots, V_{k-1}\right)+\cdots+\rho_{k-1} s_{k-1} s_{k-1}^{T} \\ i= & k-n, \cdots, k-1 \end{aligned}$
结合上述信息，简要介绍L-BFGS反演流程如下：
1）设定初始模型 m 0及修正向量的组数 n = 3 ~ 20。
2）计算目标函数梯度 ∇ φ m k，选择初始的海森矩阵 H k 0，一般选择单位矩阵。根据L-BFGS理论通过修正向量计算第 k次的海森矩阵的逆 H k - 1，进而获得当前迭代所需的拟牛顿方向，即：
p k = - H k - ∇ φ m k
3）计算修正步长：
m k + 1 = m k + α k p k
初始步长 α k为1，线性搜索得到最佳的 α k，然后得到 m k + 1。
最佳步长的计算需满足Wolf准则：
φ m k + α k p k ≤ φ m k + β ' α k ∇ φ m k p k
4）计算当前模型拟合差是否满足拟合要求，若小于则停止迭代，同时判断是否达到最大迭代次数，若两者都未满足，则跳到第二步，继续迭代。
本节公式的具体推导过程见参考文献[14]。

1.2 微动拟二维横向约束反演方法

与CSAMT拟二维横向约束反演方法相似，在微动一维正演理论的基础上，加入横向约束理论，基于L-BFGS算法实现微动拟二维横向约束算法。
微动数据中一般用 H来表示水平向分量， V表示垂直向分量，它们都是关于频率 f的函数，水平向分量与垂直向分量的傅里叶谱分别表示为 H f和 V f，两者之比即为微动H/V谱比:
H V f = H f V f
在图2所示的沉积盆地典型地质结构模型下（其中 H f和 V f—沉积盆地地表位移的水平分量和垂直分量的频谱振幅， H b和 V b—基岩的水平分量和垂直分量的频谱振幅， H r和 V r—基岩露头位移水平分量和垂直分量的频谱振幅），在沉积盆地地表的水平和垂直频谱可表示为:
图2 典型的沉积盆地地质构造[15]

Fig. 2 Typical sedimentary basin geological structure[15]

H s f = A h f H b b f + H s R f
V s f = A v f V b b f + V s R f
式（27）~（28）中： H s ( f )和 V s ( f )—沉积层表面面波的水平和垂直分量； H b ( f )和 V b ( f )—基岩位置的位移水平和垂直分量；下标 s—沉积层；下标 b—基岩；上标 b—体波；上标 R—瑞利波； f—频率。 A h f， A v ( f )—垂直入射体波的水平和垂直分量的放大因子。由式（27）和（28）得H/V谱比为：
H V f = H s f V s f = H b b f V b b f ∙ A h f + H s R f H b b f A v f + V s R f V b b f
在没有瑞利波的影响下，H/V谱比可以近似地表示为： H / V ( f ) = A h ( f ) / A v ( f )，如果瑞利波的影响越大， H s R ( f ) H b b ( f )—瑞利波的 H / V谱比，也叫瑞利波的椭率，如果瑞利波的能量占比非常大，式（29）近似为瑞利波的椭率。
本节公式的具体推导过程见参考文献[16-17]。
基于L-BFGS的横向约束反演算法流程如下：首先，输入速度的初始模型，生成横向约束矩阵，计算目标函数梯度以及雅可比矩阵，计算当前模型响应的数据拟合差，并判断是否满足拟合要求，若满足要求则停止迭代，反之，继续迭代直至满足要求。反演流程图为：
图3 反演流程图

Fig. 3 The flowchart of inversion

2 基于交叉梯度的拟二维横向约束联合反演

2.1 二维交叉梯度函数

交叉梯度函数表示为两种模型单元参量的梯度的叉积。其表达式为：
t x , y , z = ∇ m s x , y , z × ∇ m r x , y , zW
式（30）中： ∇ m s—速度的梯度； ∇ m r—电阻率的梯度[18]。
交叉梯度函数表示为两种模型单元参量的梯度叉积。在联合反演中直接加入目标函数之中，通过这种方式可以将不同物性参数的电阻率模型和速度模型连接起来。交叉梯度要求不同方法的模型值只能朝着相同方向发生改变，或者其中一种方法的模型值改变而另一种保持不变，是一种松约束的结构性参数耦合法。
拟二维反演算法本质上是一维的，该反演算法是将一条测线的测点排列在一起，在各个测点下方网格之间加入横向约束进行反演，因此将一条测线中所有测点下方的物性参数和层厚度组合到一起，就可以形成二维的网格，二维交叉梯度就是在这样的网格基础上进行计算的。

2.2 联合反演目标函数

根据前面的CSAMT和微动谱比法的L-BFGS反演思路，可以在单方法的目标函数后面加入带权重因子的交叉梯度约束实现联合反演。
因此联合反演目标函数可以写成：
φ H V m s = φ d H V + λ φ m s + α t m s , m r
φ R m r = φ d C S + λ φ m r + β t m r , m s
式（31）~（32）中：
φ d H V = d o b s H V - f H V m s T C d - 1 d o b s H V - f H V m s
φ d C S = d o b s C S - f C S m r T C d - 1 d o b s C S - f C S m r
式（33）~（34）中： φ d H V—微动联合反演目标函数的数据拟合差项； φ d C S—CSAMT联合反演目标函数的数据拟合差项； m r—电阻率的模型向量； m s—横波速度的模型向量； C d - 1—数据协方差矩阵，在公式（33）和（34）中具有不同的 C d - 1； d o b s H V—微动谱比法的观测数据向量； d o b s C S—CSAMT的观测数据向量； α， β和 λ均为权重因子，不同方法的权重因子的选取有一定的区别； t—交叉梯度算子。具体推导过程见参考文献[18]。
联合反演采用的是L-BFGS方法，因此，后续对目标函数的求导以及雅可比矩阵的计算都是与单方法反演一致的，这里不再赘述。

2.3 网格匹配

考虑到不同方法的反演过程中深度参数的变化没有相关性，在联合反演中，为了对两种模型的物性参数进行交叉梯度的耦合，对不同方法的深度网格进行了匹配处理。
最直接的办法就是将原本的单方法反演算法改写成只反演电阻率或速度，不反演地层厚度（即固定层厚）的反演算法。固定层厚的这种反演策略在只进行单一方法的反演时是没有问题的，然而在进行联合反演时，由于不同反演方法需求的层厚的固定值各不相同，反演结果往往只能满足其中一种方法的要求，无法做到同时满足要求。
为了能够利用交叉梯度对两种模型的物性参数进行约束，在单方法反演的基础上，依旧反演物性参数和地层厚度，让不同反演方法的地层层数保持一致，便于交叉梯度的网格匹配。对不同模型的深度网格进行插值的示意图如图4所示，在加入交叉梯度之前，先对电阻率模型的深度网格进行插值，得到新的电阻率模型深度网格。这样新的电阻率模型深度网格就能够与速度网格匹配上，可以进行交叉梯度的计算，从而可以得到公式（31）中的交叉梯度项。同理，只对速度网格进行插值，可得到公式（32）中的交叉梯度项。这样就完成了两种物性参数的模型网格匹配。
图4 深度网格插值示意图

Fig. 4 Depth grid interpolation schemes

2.4 联合反演流程

联合反演算法采取的是同步迭代的反演策略，即电阻率模型和速度模型分别进行迭代更新，达到一定次数后进行交叉梯度的计算，分别对两个模型进行约束，然后再继续进行迭代更新。当迭代次数超过最大迭代次数时，判断目标函数和拟合差是否满足条件，满足条件则结束联合反演，反之，继续进行迭代。
联合反演流程如下：
图5 联合反演流程示意图

Fig. 5 The flowchart of joint inversion

3 算例分析

3.1 背斜体模型

建立如图6所示的背斜体模型。电阻率背斜体理论模型从浅到深三层电阻率分别为 100、 10和 1   000   Ω ∙ m。速度背斜体理论模型从浅到深三层横波速度分别为 900、 1   500和 2   500   m ∙ s - 1。两个模型的背斜体底界面所在深度均为310 m。
图6 背斜体理论模型

a—电阻率理论模型；b—速度理论模型；c—测点布设。

Fig. 6 Theoretical model of anticline

a-Theoretical model of resistivity；b-Theoretical model of velocity；c-Layout of measurement points.

微动信号的接收器沿测线方向每隔40 m放置一个，共计21个（4 780~5 620 m）；水平电偶极源放置在距离测线5 km处，长度为500 m，沿测线方向设置21个测点，间隔40 m（4 780~5 620 m），发射频率设置为13个，分别为1、2、4、8、16、32、64、128、256、512、1 024、2 048和 4 096 Hz。
根据趋肤深度公式[14]，可以知道：
H = 356 ρ 1 f
式（35）中： H—趋肤深度，m； ρ 1—表层电阻率，Ω∙m； f—频率，Hz。
对于表层电阻率为100 Ω ∙ m的地下介质，可观测到的最大趋肤深度为：
H = 356 ρ 1 f = 356 100 1 = 3   560
对于选取频点的穿透深度从地表到地下深度依次为：55、78、111、157、222、314、445、629、890、1 258、1 780、2 517和3 560 m。
根据此深度规律，可以推断，对于异常体规模较小比如存在薄层、异常体埋藏深度超过3 560 m的测区，本观测装置就难以分辨。而对于深度在500 m之内且具有一定规模的异常体，本观测装置采集的数据具有很好的分辨效果。
在理论模型正演响应数据中加入5 %高斯随机噪声作为反演数据。取发射频率，即13个频率进行反演。初始电阻率模型设置为17层的均匀半空间模型，电阻率为 100   Ω ∙ m，层厚度设置16个，分别为20、21、22、23、24、25、26、27、28、29、30、31、32、33、34和35 m。初始速度模型设置为17层的均匀半空间模型，纵波速度为 1   700   m ∙ s - 1，横波速度为 1   000   m ∙ s - 1，层厚度均为20 m，设计不同的横向约束权重（ λ）对反演效果进行对比验证。横向约束权重分别为10、100、1 000和10 000，对比不同横向约束权重的反演结果，从中选出最为合适的反演结果，或者增加更多横向约束权重值，进行更加细致的对比。CSAMT和微动谱比法的反演结果见图7和8。
图7 未加横向约束时背斜体模型反演结果对比

a—CSAMT单独反演结果；b—微动谱比法单独反演结果；c—CSAMT联合反演结果；d—微动谱比法联合反演结果。

Fig.7 Comparison of inversion results of anticlinal model without laterally constraint

a-Separate result of CSAMT；b-Separate result of microtremor spectrum ratio method；c-Joint inversion result of CSAMT；d-Joint inversion result of microtremor spectrum ratio method.

图8 加入横向约束后背斜体模型反演结果对比

a—CSAMT单独反演结果；b—微动谱比法单独反演结果；c—CSAMT联合反演结果；d—微动谱比法联合反演结果。

Fig. 8 Comparison of inversion results of anticlinal model with laterally constraint

a-Separate result of CSAMT；b-Separate result of microtremor spectrum ratio method；c-Joint inversion result of CSAMT；d-Joint inversion result of microtremor spectrum ratio method.

如图7和8所示，CSAMT单独反演的结果只能恢复异常体的基本位置，背斜形态并不明显，且深部地层的恢复效果较于真实模型相差甚远；而联合反演的结果明显改善了异常体的背斜形态以及深部地层的恢复情况。微动谱比法单方法反演对异常体的速度和边界的恢复效果较好，且速度联合反演结果相比于单方法对模型边界的刻画更加准确，改善了反演结果中少许的低速假异常。横向约束的加入使反演结果的横向连续性均得到有效的改善，其中横向约束和联合反演的共同作用使得反演结果更接近真实模型。

3.2 菱形体模型

建立如图9所示两端对称的菱形体模型。菱形体的电阻率为 10   Ω ∙ m，背景电阻率为 100   Ω ∙ m。速度菱形体理论模型从浅到深3层横波速度分别为 900、 1   500和 2   500   m ∙ s - 1。横向约束权重分别为10、100、1 000和10 000，对比不同横向约束权重的反演结果，从中选出最为合适的反演结果或者增加更多横向约束权重值，进行更加细致的对比。
图9 菱形体理论模型

a—电阻率理论模型；b—速度理论模型；c—测点布设。

Fig. 9 Theoretical model of rhombus

a-Theoretical model of resistivity；b-Theoretical model of velocity；c-Layout of measurement points.

CSAMT与微动的观测系统和反演初始模型的设置与第3.1节相同。
如图10和11所示，CSAMT单独反演的结果对异常体的恢复效果极差；只进行联合反演的结果能够反映出异常体大概位置和两端对称的形态，但边界刻画不够清晰；只加入横向约束的反演结果较好的恢复出异常体的形态，清楚的刻画出异常体的顶界面，但底界面的刻画不够准确，异常体下方的电阻率变化很小，说明单方法反演有所局限；加入横向约束后的联合反演结果完美的恢复出异常体的位置与形态，对于深部地层的恢复也有所改善，但受到单方法的影响较大，仍然存在深部地层电阻率变化很小的特征。
图10 未加横向约束时菱形体模型反演结果对比

a—CSAMT单独反演结果；b—微动谱比法单独反演结果；c—CSAMT联合反演结果；d—微动谱比法联合反演结果。

Fig. 10 Comparison of inversion results of rhomboid model without laterally constraint

a-Separate result of CSAMT；b-Separate result of microtremor spectrum ratio method；c-Joint inversion result of CSAMT；d-Joint inversion result of microtremor spectrum ratio method.

图11 加入横向约束后菱形体模型反演结果对比

a—CSAMT单独反演结果；b—微动谱比法单独反演结果；c—CSAMT联合反演结果；d—微动谱比法联合反演结果。

Fig. 11 Comparison of inversion results of rhomboid model with laterally constraint

a-Separate result of CSAMT；b-Separate result of microtremor spectrum ratio method；c-Joint inversion result of CSAMT；d-Joint inversion result of microtremor spectrum ratio method.

微动谱比法单独反演的结果对异常体的恢复效果很差，仅能刻画出菱形体的底界面；只加入横向约束或只进行联合反演均能有效的压制单独反演结果中大量的假异常，但对于菱形体边界的刻画有所偏移；加入横向约束的联合反演结果对真实模型的速度和边界恢复的较好，但相比于只加入横向约束或只进行联合反演的反演结果没有明显的改善。
总体来看，不管是背斜体模型还是菱形体模型，开发的横向约束联合反演算法均能恢复出地下真实模型的形态与真值，且均能对单方法的反演结果有一定的改善效果，说明在层状模型反演中加入交叉梯度是有效的。

4 实测数据应用

4.1 测区概况

测区位于北京市延庆区西北部延庆农场附近。查阅相关地质资料可知，该测区北部附近深部地层由新到老为：第四系、下白垩系东岭台组、上侏罗系髫髻山组和蓟县系雾迷山组[20]。
各地层岩性及电阻率分布为：
1）第四系：该组岩性主要以黏性土、砾石层与粗中砂层为主，厚度介于500~1 000 m之间。电阻率范围介于 20~ 80   Ω ∙ m之间。
2）下白垩系东岭台组：该组在延庆县城北部山区有出露，以流纹岩、砂岩、砂砾岩和火山岩为主。电阻率范围介于71~ 174   Ω ∙ m之间。
3）上侏罗系髫髻山组：以砂岩、白云岩、砾岩和火山岩为主。电阻率范围介于71~ 174   Ω ∙ m之间。
4）蓟县系雾迷山组：以砂泥质白云岩为主。电阻率范围介于540~ 1   000   Ω ∙ m之间。工作区内不同地层间电性参数有明显差异，这是电磁法勘探适用的基础。

4.2 野外数据采集

CSAMT的测线布设如图12所示。本测线采用接地导线场源进行观测，实际观测收发距为10 km，发射偶极长度为1 km，共计34个测点，测点间距40 m。观测仪器为GDP32-II，观测频点共计17个（ 2 - 3~ 2 10   H z），观测电磁场分量为 E x和 H y。
图12 延庆地区地质简图[19]

1—第四系；2—白垩系；3—侏罗系；4—蓟县系雾迷山组；5—长城系；6—太古界；7—燕山期侵入岩；8—正断层；9—逆断层；10—山区断层；11—推测断裂；12—山区界；13—重力计算区；14—热流计算孔。

Fig. 12 Geological map of Yanqing area

1-Quaternary system；2-Cretaceous system；3-Jurassic system；4-Yijinshan formation of Jixian system；5-Changcheng System；6-Archean system；7-Intrusive rocks of Yanshan period；8-Normal fault；9-Reverse fault；10-Mountainous fault；11-Speculative fault；12-Mountainous boundary；13-Gravity calculation area；14-Heat flow calculation hole.

微动探测的测线布设如图13所示。野外采集所用仪器为Smart Solo，依据微动探测空间自相关法原理，采用了环形台站的布站方式进行数据的采集，一个点布设13台台站，采集时长为1 h15 min，观测数据为微动的三分量数据。用来联合反演的实测数据采用的是距离CSAMT各个测点位置最近的台站数据，共计34个台站。两条测线均沿南北方向走线。
图13 实测数据测点分布图

Fig. 13 Distribution of measured data points

横向约束权重分别为10、100、1 000和 10 000，对比不同横向约束权重的反演结果，从中选出最为合适的反演结果，或者增加更多横向约束权重值，进行更加细致的对比。
直接利用采集到的实测数据可以得到拟断面图，CSAMT与微动的拟断面图表明：在测线下方存在层状连续的地层，适合利用开发的算法进行反演。
图14 实测数据视电阻率拟断面图

Fig. 14 Pseudo-section measured apparent resistivity

图15 实测数据阻抗相位拟断面图

Fig. 15 Pseudo-section measured impendence phase

4.3 反演结果

未加横向约束时，电阻率的联合反演结果RMS降至3.3，加入横向约束后，电阻率的联合反演结果RMS降至1.8。图16的反演结果表明：未加横向约束时，电阻率反演结果整体上分布比较凌乱，无法体现出测线下方地层的分层性；加入横向约束后，电阻率反演结果的横向分布变得连续，地层之间的分界线刻画得非常清晰，可以大致看出测线下方的地层分为3层，由浅至深电阻率呈现逐渐增大的趋势。
图16 实测数据HVSR值拟断面图

Fig. 16 HVSR-profiling of measured data

未加横向约束时，速度的联合反演结果RMS降到1.1，加入横向约束后，速度联合反演结果RMS降到1.0。图17的反演结果表明，未加横向约束时，速度联合反演结果的整体分布较为杂乱，虽然也能看出测线下方的浅部地层是低速层，深部地层是高速层，但是地层之间的分界并不明显；加入横向约束后，速度联合反演结果的横向连续性得到明显的改善，浅部低速层和深部高速层的特征依旧存在，且不同速度的地层之间界限更加清晰。相比于电阻率的横向约束联合反演结果，速度的联合反演结果分层不够明显，原因在于未加横向约束时，速度的反演结果RMS已经十分接近1.0，这时横向约束的加入对速度反演结果的优化效果就不如电阻率反演结果明显。
图17 实测数据联合反演电阻率结果

a—未加横向约束；b—加入横向约束。

Fig. 17 Joint inversion resistivity results of measured data

a-Without laterally constraint；b-With laterally constraint.

图18 实测数据联合反演速度结果

a—未加横向约束；b—加入横向约束。

Fig. 18 Joint inversion velocity results of measured data

a-without laterally constraint；b-with laterally constraint.

实测数据的联合反演结果表明：在测线下方的地层存在着连续性；在联合反演中加入横向约束可以进一步对反演结果进行优化。

5 结论

1）引入横向约束理论，分别实现了CSAMT和微动谱比法的拟二维横向约束反演；加入交叉梯度函数，实现了CSAMT和微动谱比法的拟二维横向约束联合反演，并验证了算法的可行性和适用性。
2）设计多个理论模型合成数据，进行联合反演计算，结果表明：联合反演可以有效的改善和提升电阻率的反演结果，一定权重的横向约束可以进一步改善两个方法的反演结果。
3）实测数据的反演结果表明，联合反演可以结合CSAMT和微动方法各自的特点，使反演结果获得更好的一致性；横向约束的加入能够使反演结果更加连续，分层性更好，验证了联合反演算法的实用性。
CSAMT反演算法只对电阻率数据进行了处理和反演，忽略了相位数据，可能是CSAMT反演效果不好的原因，后续会考虑将相位数据加入CSAMT反演中，进一步提升反演算法的可靠性。

参考文献

原文顺序 | 文献年度倒序 | 文中引用次数倒序
1
路俊涛, 雷达, 任浩. 含激电效应的航空瞬变电磁数据横向约束反演[J]. 地球物理学进展, 2022, 37(1):1-14.

LU Juntao, LEI Da, REN Hao. Laterally constrained inversion of airborne transient electromagnetic data with induced-polarization effects[J]. Progress in Geophysics, 2022, 37(1):220-229 (in Chinese).

2
蔡晶, 齐彦福, 殷长春. 频率域航空电磁数据的加权横向约束反演[J]. 地球物理学报, 2014, 57(3):953-960.

CAI Jing, QI Yanfu, YIN Changchun. Weighted Laterally-constrained inversion of frequency-domain airborne EM data[J]. Chinese Journal of Geophysics, 2014, 57(3):953-960 (in Chinese).

3
WANG R, YIN C, WANG M, et al. Laterally constrained inversion for CSAMT data interpretation[J]. Journal of Applied Geophysics, 2015,(121):63-70.

4
殷长春, 邱长凯, 刘云鹤, 等. 时间域航空电磁数据加权横向约束反演[J]. 吉林大学学报(地球科学版), 2016, 46(1):254-261.

YIN Changchun, QIU Changkai, LIU Yunhe, et al. Weighted laterally-constrained inversion of time-domain airborne electromagnetic data[J]. Journal of Jilin University(Earth Science Edition), 2016, 46(1):254-261 (in Chinese).

5
Dobroka M, Gyulai A, Ormos T, et al. Joint inversion of seismic and geoelectric data recorded in an underground coal mine[J]. Geophysical Prospecting, 1991, 39(5):643-665.

6
杨振武, 王家映, 张胜业, 等. 一维大地电磁和地震数据联合反演方法研究[J]. 石油地球物理勘探, 1998, 33(1):78-88.

YANG Zhenwu, WANG Jiaying, ZHANG Shengye, et al. Research on joint inversion of 1-D magnetoteliuric and seismic data[J]. Oil Geophysical Prospecting, 1998, 33(1):78-88 (in Chinese).

7
过仲阳, 王家林, 吴健生. 应用遗传算法联合反演地震—大地电磁测深数据[J]. 石油物探, 1999, 38(1):102-108.

GUO Zhongyang, WANG Jialin, WU Jiansheng. Joint inversion of seismic and magnetotelluric sounding data using genetic algorithm[J]. Geophysical Prospecting for Petroleum, 1999, 38(1):102-108 (in Chinese).

8
杨辉, 王家林, 吴健生, 等. 大地电磁与地震资料仿真退火约束联合反演[J]. 地球物理学报, 2002, 45(5):723-734.

YANG Hui, WANG Jialin, WU Jiansheng, et al. Constrained joint inversion of magnetotelluric and seismic data using simulated annealing algorithm[J]. Chinese Journal of Geophysics, 2002, 45(5):723-734 (in Chinese).

9
Heincke B, Hobbs R. Joint inversion of MT,gravity and seismic data applied to sub-basalt imaging[C]// SEG2006 Expanded Abstract,SEG, 2006,25:784-787.

10
胡英才, 程纪星. 巴楚地区可控源音频大地电磁法与电性源瞬变电磁法正演模拟对比研究[J]. 世界核地质科学, 2024, 41(6):1141-1155.

HU Yingcai, CHENG Jixing. Comparative study on forward simulation of controllable source audio magnetotelluric method and electrical source transient electromagnetic method in Bachu area[J]. 世界核地质科学, 2024, 41(6):1141-1155 (in Chinese).

11
Auken E, Christiansen A V. Layered and Laterally Constrained 2D Inversion of Resistivity Data[J]. Geophysics, 2004, 69(3):752-761.

12
Auken E, Foged N, Sorensen K I. Model recognition by 1-D laterally constrained inversion of resistivity data[C]// 9th Meeting,Environmental and Engineering Geophysical Society-European Section,Proceedings, 2002,241-244.

13
Christiansen A V, Auken E, Foged N, et al. Mutually and laterally constrained inversion of CVES and TEM data:A case study[J]. Near Surface Geophysics, 2007, 5(2):115-123.

14
王帅. 基于横向约束的可控源音频大地电磁法拟三维反演研究[D]. 北京: 中国地质大学(北京), 2022.

WANG Shuai. Research on CSAMT quasi-three-dimensional inversion based on laterally constrained inversion[D]. Beijing: China University of Geosciences(Beijing), 2022 (in Chinese).

15
Nakamura Y. Clear identification of fundamental idea of Nakamura's technique and its applications[C]// Proccedings of the 12th World Conference on Earth-quake Engineering, 2000.

16
张一梵. 微动勘探法在浅层探测中的研究与应用[D]. 北京: 中国地质大学(北京), 2019.

ZHANG Yifan. Research and application of microtremor survey method in shallow exploration[D]. Beijing: China University of Geosciences(Beijing), 2019 (in Chinese).

17
秦彤威, 王少曈, 冯宣政, 等. 微动H/V谱比方法[J]. 地球与行星物理论评, 2021, 52(6):587-622.

QIN Tongwei, WANG Shaotong, FENG Xuanzheng, et al. A review on microtremor H/V spectral ratio method[J]. Reviews of Geophysics and Planetary Physics, 2021, 52(6):587-622 (in Chinese).

18
彭淼, 谭捍东, 姜枚, 等. 基于交叉梯度耦合的大地电磁与地震走时资料三维联合反演[J]. 地球物理学报, 2013, 56(8):2728-2738.

PENG Miao, TAN Handong, JIANG Mei, et al. Three-dimensional joint inversion of magnetotelluric and seismic travel time data with cross-gradient constraints[J]. Chinese Journal of Geophysics, 2013, 56(8):2728-2738 (in Chinese).

19
雷晓东, 李晨, 王立发, 等. 延庆盆地大地热流异常及其构造背景[J]. 地球物理学报, 2022, 65(9):3405-3418.

LEI Xiaodong, LI Chen, WANG Lifa, et al. Genesis of the heat flow anomaly in Yanqing basin,NW Beijing[J]. Chinese Journal of Geophysics, 2022, 65(9):3405-3418 (in Chinese).

20
张进平. 浅析北京市延庆区西北部地热地质特征[J]. 城市地质, 2019, 14(1):26-33.

ZHANG Jinping. Analysis on geothermal geology characteristics in the northwest of Yanqing district,Beijing[J]. Urban Geology, 2019, 14(1):26-33 (in Chinese).

Options
文章导航

/

〈 〉