多辐射场源地空瞬变电磁一维反演方法研究

图1 单辐射场源剖分俯视图

Fig.1 The platform of single-radiation source subdivision

1.1.1 各向同性水平层状大地频率域瞬变响应

现有各向同性水平层状大地模型如图2所示,其层数为n,各层电阻率分别为ρ₁,ρ₂,ρ₃,…,ρ_n,对应各层厚度为h₁,h₂,h₃,…,h_n_-1,直角坐标系的原点为O,将长度为ds的电偶极子AB置于地表,其中心位于点O,位于空中的测点M在地表的投影为P,距原点O的距离|OP|=r,线段OP与x轴正方向即电偶极子AB的夹角为φ,点M到地表的距离为z。

图2

图2 层状大地地表电偶极子坐标系示意

Fig.2 The diagram of layered earth surface electric dipole coordinate system

由此可知,在距地表高度为h处放置的电偶极子对于各向同性水平层状大地在空中产生的电矢量位为^[14,15]

(1) $A_{x 0} = \frac{P_{E}}{4 π} \int_{0}^{\infty} [\frac{λ}{u_{0}} e^{- u_{0} | z + h |} + \frac{- \frac{λ}{u_{0}} \frac{u_{1}}{R_{1}} + λ}{u_{0} + \frac{u_{1}}{R_{1}}} e^{- u_{0} h} e^{u_{0} z}] J_{0} (λr) d λ 。$

通过对电矢量位进行一系列变换(见附录),再求取旋度,并做进一步的整理可得地表电偶极子在空中产生的频率域磁感应强度表达式为

(2) $\{\begin{array}{l} B_{x} (ω) = \frac{P_{E}}{2 π} μ_{0} \sin φ \cos φ \int_{0}^{\infty} \frac{λ^{2}}{λ + u_{1} / R_{1}} e^{λz} J_{0} (λr) d λ - \frac{P_{E}}{2 π} \frac{μ_{0} \sin φ \cos φ}{r} \int_{0}^{\infty} \frac{2 λ}{λ + u_{1} / R_{1}} e^{λz} J_{1} (λr) d λ \\ B_{y} (ω) = \frac{P_{E}}{2 π} μ_{0} si n^{2} φ \int_{0}^{\infty} \frac{λ^{2}}{λ + u_{1} / R_{1}} e^{λz} J_{0} (λr) d λ + \frac{P_{E}}{2 π} \frac{μ_{0} (co s^{2} φ - si n^{2} φ}{r} \int_{0}^{\infty} \frac{λ}{λ + u_{1} / R_{1}} e^{λz} J_{1} (λr) d λ \\ B_{z} (ω) = \frac{P_{E}}{2 π} μ_{0} \sin φ \int_{0}^{\infty} \frac{λ^{2}}{λ + u_{1} / R_{1}} e^{λz} J_{1} (λr) d λ \end{array},$

(3) $其中 : \begin{matrix} R_{1} = cth [u_{1} h_{1} + arcth \frac{u_{1}}{u_{2}} cth (u_{2} h_{2} + \dots + arcth \frac{u_{n - 1}}{u_{n}})] 。 \end{matrix}$

u_i= $\sqrt[]{λ^{2} + k_{i}^{2}}$ , $k_{i}^{2}$ =-iωμσ_i-ω²με,在导电介质中,若忽略位移电流 $k_{i}^{2}$ =-iωμσ_i,σ_i是第i层地层的电导率。

在计算单辐射场源地空瞬变电磁响应时,首先要运用剖分的思想,可以把一个电性源分成若干小段,而每一小段又都可以被看成是一个电偶极子。因此,在计算单辐射场源地空瞬变响应时就在电偶极子的基础上进行计算,将每个电偶极子在测点上产生的瞬变响应都计算之后,然后再将剖分成的每个电偶极子产生的瞬变响应均换算到整体坐标系下,最后在整体坐标系下,对所剖分的各个电偶极子产生瞬变响应进行叠加、求和,就可得层状大地单辐射场源在空中产生的频率域磁感应强度表达式为

(4) $\{\begin{array}{l} B_{x} (ω) = \overset{m}{\sum_{j = 1}} [- \frac{P_{Ej}}{2 π} \frac{μ_{0} \sin φ_{j} \cos φ_{j}}{r_{j}} \int_{0}^{\infty} \frac{2 λ}{λ + u_{1} / R_{1}} e^{λz} J_{1} (λ r_{j}) d λ + \frac{P_{Ej}}{2 π} μ_{0} \sin φ_{j} \cos φ_{j} \int_{0}^{\infty} \frac{λ^{2}}{λ + u_{1} / R_{1}} e^{λz} J_{0} (λ r_{j}) d λ] \\ B_{y} (ω) = \overset{m}{\sum_{j = 1}} [\frac{P_{Ej}}{2 π} \frac{μ_{0} (co s^{2} φ_{j} - si n^{2} φ_{j}}{r_{j}} \int_{0}^{\infty} \frac{λ}{λ + u_{1} / R_{1}} e^{λz} J_{1} (λ r_{j}) d λ + \frac{P_{Ej}}{2 π} μ_{0} si n^{2} φ_{j} \int_{0}^{\infty} \frac{λ^{2}}{λ + u_{1} / R_{1}} e^{λz} J_{0} (λ r_{j}) d λ] \\ B_{z} (ω) = \overset{m}{\sum_{j = 1}} \frac{P_{E_{j}}}{2 π} μ_{0} \sin φ_{j} \int_{0}^{\infty} \frac{λ^{2}}{λ + u_{1} / R_{1}} e^{λz} J_{1} (λ r_{j}) d λ \end{array} 。$

1.1.2 频率域瞬变响应与时间域瞬变响应的变换

当发射波形为阶跃波时可知

(5) $I = \{\begin{array}{l} I_{0}, t ＜ 0; \\ 0, t \geq 0 \end{array} 。$

由频谱分析理论,可知谐变场量F(ω)与时间场量f(t)之间的转换关系为^[14,16-17]

(6) $f (t) = \frac{1}{2 π} \int_{-}^{\infty} \frac{F (ω)}{- i ω} e^{- i ωt} d ω,$

又由积分变换可得到频率域与时间域磁感应强度之间的转换关系^[18,19]:

(7) $B (t) = \frac{1}{2 π} \int_{-}^{\infty} \frac{B (ω)}{- i ω} e^{- i ωt} d ω 。$

最后,经过化简就可得

(8) $\{\begin{array}{l} B (t) = B_{0} - \frac{2}{π} \int_{0}^{\infty} \frac{Im B (ω)}{ω} \cos (ωt) d ω \\ B (t) = \frac{2}{π} \int_{0}^{\infty} \frac{Re B (ω)}{ω} \sin (ωt) d ω \end{array} 。$

进而可得对应的时间导数:

(9) $\{\begin{array}{l} \frac{\partial B (t)}{\partial t} = - \frac{2}{π} \int_{0}^{\infty} Im B (ω) \sin (ωt) d ω \\ \frac{\partial B (t)}{\partial t} = - \frac{2}{π} \int_{0}^{\infty} Re B (ω) \cos (ωt) d ω \end{array} 。$

根据式(8)、式(9)所体现的频率域瞬变响应向时间域瞬变响应转换的关系,结合式(4)就可得到单辐射场源情况下各向同性水平层状大地时间域瞬变响应。

1.2 多辐射场源地空瞬变电磁法正演原理

前面已经研究了单辐射场源地空瞬变电磁法的正演问题,现在对于多辐射场源的情况,可以将其看作单辐射场源的叠加。在这里应用到坐标系的转换,也就是说在计算出每个电偶极子源在其自身的坐标系下产生的频率域瞬变响应后,要将其响应均转换到整体坐标系下并叠加,从而得到多辐射场源的频率域瞬变响应,然后再根据时频转换的相关理论进行转换,从而得到多辐射场源时间域瞬变响应。

1.2.1 各向同性水平层状大地频率域瞬变响应

针对多辐射场源的情况,本文中将继续采用剖分叠加的方法,以单辐射场源时的瞬变响应为基础,从而得到多辐射场源情况下的瞬变响应。

图3所示为多辐射场源所在坐标系俯视图,如图所示,本文以其中一个电性源所在的坐标系为整体坐标系,图3中就是以源A₁B₁所在的坐标系为准;而后,将所有源产生的瞬变响应转换到整体坐标系下并叠加,便可得到各向同性水平层状大地多辐射场源在空中产生的频率域瞬变响应如式(10):

图3

图3 多辐射场源所在坐标系俯视图

Fig.3 The platform of multi-radiation source subdivision

(10) $\{\begin{array}{l} B_{x} (ω) = \overset{n}{\sum_{i = 1}} \overset{n_{j}}{\sum_{j = 1}} \cos α_{i} [\begin{array}{l} - \frac{P_{Eij}}{2 π} \frac{μ_{0} \sin φ_{ij} \cos φ_{ij}}{r_{ij}} \int_{0}^{\infty} \frac{2 λ}{λ + u_{1}} e^{λz} J_{1} (λ r_{ij}) d λ \\ + \frac{P_{Eij}}{2 π} μ_{0} \sin φ_{ij} \cos φ_{ij} \int_{0}^{\infty} \frac{λ^{2}}{λ + u_{1}} e^{λz} J_{0} (λ r_{ij}) d λ \end{array}] \\ - \overset{n}{\sum_{i = 1}} \overset{n_{i}}{\sum_{j = 1}} \sin α_{i} [\begin{array}{l} \frac{P_{Eij}}{2 π} \frac{μ_{0} (co s^{2} φ_{ij} - si n^{2} φ_{ij})}{r_{ij}} \int_{0}^{\infty} \frac{λ}{λ + u_{1} / R_{1}} e^{λz} J_{1} (λ r_{ij}) d λ \\ + \frac{P_{Eij}}{2 π} μ_{0} si n^{2} φ_{ij} \int_{0}^{\infty} \frac{λ^{2}}{λ + u_{1} / R_{1}} e^{λz} J_{0} (λ r_{ij}) d λ \end{array}] \\ B_{y} (ω) = \overset{n}{\sum_{i = 1}} \overset{n_{j}}{\sum_{j = 1}} \sin α_{i} [\begin{array}{l} - \frac{P_{Eij}}{2 π} \frac{μ_{0} \sin φ_{ij} \cos φ_{ij}}{r_{ij}} \int_{0}^{\infty} \frac{2 λ}{λ + u_{1} / R_{1}} e^{λz} J_{1} (λ r_{ij}) d λ \\ + \frac{P_{Eij}}{2 π} μ_{0} \sin φ_{ij} \cos φ_{ij} \int_{0}^{\infty} \frac{λ^{2}}{λ + u_{1} / R_{1}} e^{λz} J_{0} (λ r_{ij}) d λ \end{array}] \\ + \overset{n}{\sum_{i = 1}} \overset{n_{i}}{\sum_{j = 1}} \cos α_{i} [\begin{array}{l} \frac{P_{Eij}}{2 π} \frac{μ_{0} (co s^{2} φ_{ij} - si n^{2} φ_{ij})}{r_{ij}} \int_{0}^{\infty} \frac{λ}{λ + u_{1} / R_{1}} e^{λz} J_{1} (λ r_{ij}) d λ \\ + \frac{P_{Eij}}{2 π} μ_{0} si n^{2} φ_{ij} \int_{0}^{\infty} \frac{λ^{2}}{λ + u_{1} / R_{1}} e^{λz} J_{0} (λ r_{ij}) d λ \end{array}] \\ B_{z} (ω) = \overset{n}{\sum_{i = 1}} \overset{n_{i}}{\sum_{j = 1}} \frac{P_{Eij}}{2 π} μ_{0} \sin φ_{ij} \int_{0}^{\infty} \frac{λ^{2}}{λ + u_{1} / R_{1}} e^{λz} J_{1} (λ r_{ij}) d λ \end{array} 。$

1.2.2 频率域瞬变响应与时间域瞬变响应的变换

在单辐射场源地空瞬变响应的计算中,采用正余弦变换的方法对其频率域瞬变响应和时间域瞬变响应进行转换,而对于多辐射场源产生的瞬变响应这种方法同样适用,因此,在进行多辐射场源的瞬变响应的频率域与时间域之间的转换本文仍然采用正余弦变换方法,即式(8)与式(9)。

1.3 多辐射场源地空瞬变电磁法一维Occam反演原理

Occam反演方法其实是一种带平滑约束的最小二乘反演方法,即在一定的拟合误差标准下使模型的粗糙度最小。

模型的粗糙度表示的是模型的光滑程度。模型粗糙度的表示形式为

(11) $R = \int {(\frac{d m}{d z})}^{2} d z,$

其中:z是深度,m(z)是模型每个深度下与之相对应的地电参数,在本文中则为电阻率。为使计算简便,将模型的地电参数表示为

(12) $m (z) = m_{i}, i = 1,2, \dots, N;$

随后,将模型粗糙度即式(11)进行离散化,可得

(13) $R_{1} = \overset{N}{\sum_{i = 2}} (m_{i} - m_{i - 1})^{2} 。$

可以将式(13)由离散形式改写成矩阵形式:

(14) $R_{1} = ‖ \partial m ‖^{2},$

其中$\partial$为粗糙度矩阵,表示为

(15) $\partial = [\begin{array}{l} 0 & 0 \\ - 1 & 1 \\ - 1 & 1 \\ \dots & \dots \\ 0 & - 1 & 1 \end{array}] 。$

对于模型计算数据与实测数据的拟合差X²,可以由

(16) $X^{2} = \overset{M}{\sum_{i = 1}} \frac{(d_{i} - F_{i} {(m))}^{2}}{σ_{i}^{2}}$

得到。其中:d_i(i=1,2,3,…,M)是观测数据,F_i(m)(i=1,2,3,…,M)是正演数据,σ_i(i=1,2,3,…,M)是对应数据的标准差。而Occam反演方法就是对于测定的数据结果及其标准差,找到一个粗糙度尽可能的小而拟合差又能够达到期望能够达到的值,因此引入拉格朗日乘子μ继而建立目标函数:

(17) $X^{2} = \overset{M}{\sum_{i = 1}} \frac{(d_{i} - F_{i} {(m))}^{2}}{σ_{i}^{2}},$

其中:d_i(i=1,2,3,…,M)是观测数据,F_i(m)(i=1,2,3,…,M)是正演数据,σ_i(i=1,2,3,…,M)是对应数据的标准差。

根据Occam反演方法的基本原理,可以建立目标函数^[20-23,29]:

(18) $U = R_{1} + μ^{- 1} {‖ Wd - WF [m] ‖^{2} - X_{\infty}^{2}},$

式中:μ是拉格朗日乘子,d是观测数据向量,F是正演算子, $X_{\infty}^{2}$ 是反演所要求达到的拟合差,W=diag{1/σ₁,1/σ₂,…,σ_j,…,1/σ_M},σ_j是第j个数据的标准差,R₁ 是模型粗糙度。

首先构建初始模型m₁,然后在其附近将该非线性问题进行线性化,则可以得到结果^{[24,25,26,27,28]}:

(19) $F [m_{1} + Δ] = F [m_{1}] + J_{1} Δ 。$

式中:Δ=m₂-m₁是模型修改量,m₂是反演一次迭代的解,J₁是初始模型m₁的偏导数矩阵,其元素就可表示为J_ij=$\partial$F_i[m]/ $\partial$m_j。

将式(19)代入式(20 )并整理就可得到:

(20) $\begin{array}{l} U = ‖ \partial m_{2} ‖^{2} + μ^{- 1} {‖ W (d - F [m_{1}] + J_{1} m_{1}) \\ - W J_{1} m_{2} ‖^{2} - X_{\infty}^{2}} 。 \end{array}$

反演迭代的目的就是使目标函数达到最小,因此通常求取其梯度并令梯度值为零即▽U=0。由此便可得到

(21) $μ^{- 1} (W J_{1})^{T} W J_{1} m_{2} - μ^{- 1} (W J_{1})^{T} W d_{1} + \partial^{T} \partial m_{2} = 0 。$

在Occam反演求μ的方法是在每次迭代中搜索找到一个μ^*,使得数据残差平方和的值最小,即

(22) $X_{k + 1}^{2} (μ^{*}) = \min {‖ Wd - WF [m_{k + 1} (μ^{*})] ‖^{2}} 。$

而在这里,对于μ值的求取方法本文采用的线性搜索的办法来找到能够满足条件的数值。在求解过程中,满足条件的值可能会出现多个,这时选取其中最大的,因为可以使得模型的粗糙度最小。

2 模型计算

针对本章的模型计算,其场源的位置及测点坐标均如图4所示,将两个长度均为1 000 m的电性源平行放置,平行电性源的电流方向如图中箭头所示,二者的电流大小均为100 A,飞行高度为100 m,偏移距为500 m。观测点M在xyz坐标系中对应的坐标为(300 m,400 m,-100 m)。

图4

图4 场源及接收点的位置示意

Fig.4 Field source location and measuring point

2.1 典型三层模型

运用各个模型所得瞬变响应B_p(t)(p=x,y,z)进行反演计算时,以电阻率100 Ω·m为的半空间模型为初始模型,共30层,迭代的终止条件为拟合差小于5%或最大迭代次数为30次。以下模型系统装置如无特殊说明都是如此。

图5b、图2b中可以看出两种模型的反演的计算结果与理论模型拟合情况。图5c、d、e,图6c、d、e分别为反演模型与理论模型关于瞬变响应B_p(t)(p=x,y,z)的拟合曲线。从理论模型和反演模型计算所得的B_p(t)(p=x,y,z)曲线图中可以看出B_p(t)(p=x,y,z)的曲线形态越简单其拟合效果越好,通过对比可以看出在B_p(t)(p=x,y,z)的曲线图中,B_y(t)的曲线拟合的最好,几乎完全重合,而B_p(t)(p=x,z)的曲线在中间拐点处拟合程度较差,误差较大。对比图5b、图6b可以发现,在此情况下运用瞬变响应B_y(t)所得的反演结果与理论模型拟合相对较好。因此,在反演计算的过程中应尽量运用曲线形态相对简单的瞬变响应数据进行反演计算。

图5

图5 A型模型示意图及反演结果对比

Fig.5 The diagram and inversion results comparison chart of A-type

图6

图6 Q型模型示意图及反演结果对比

Fig.6 The diagram and inversion results comparison chart of Q-type

2.2 四层模型计算

设计KH型、HK型模型,其各项参数如图7a、图8a所示。每层厚度为15 m。图7b、图8b中可以看出两种模型的反演的计算结果与理论模型拟合情况。图7c、d、e,图8c、d、e分别为反演模型与理论模型关于瞬变响应B_p(t)(p=x,y,z)的拟合曲线,从理论模型和反演模型计算所得的B_p(t)(p=x,y,z)曲线图中可以看出B_p(t)(p=x,y,z)的曲线形态越简单其拟合效果越好,通过对比可以看出在B_p(t)(p=x,y,z)的曲线图中,B_y(t)的曲线拟合的最好,几乎完全重合,而B_p(t)(p=x,z)的曲线在中间拐点处拟合程度较差,误差较大。对比图7b、图8b可以发现,在此情况下运用瞬变响应B_y(t)所得的反演结果与理论模型拟合相对较好。

图7

图7 KH型模型示意及结果对比

Fig.7 The diagram and inversion results comparison chart of KH-type

图8

图8 HK型模型示意及结果对比

Fig.8 The diagram and inversion results comparison chart of HK-type

2.3 五层模型计算

设计五层模型,其各项参数如图9a所示。每层厚度为15 m。根据上述模型的计算,验证了多辐射场源地空瞬变电磁法一维Occam 反演方法的可靠性,因此应用该方法对五层模型进行反演计算,在计算过程中将多辐射场源地空瞬变电磁法的正演结果作为实测数据继而进行反演计算。计算结束后,将计算得到最终结果与给定的理论模型进行比较。图9b中可以看出该模型的反演的计算结果与理论模型拟合情况。图9c、d、e分别为反演模型与理论模型关于瞬变响应B_p(t)(p=x,y,z)的拟合曲线,从理论模型和反演模型计算所得的B_p(t)(p=x,y,z)曲线图中可以看出B_p(t)(p=x,y,z)的曲线形态越简单其拟合效果越好,通过对比可以看出:在B_p(t)(p=x,y,z)的曲线图中,B_y(t)的曲线拟合的最好,几乎完全重合,而B_p(t)(p=x,z)的曲线在中间拐点处拟合程度较差,误差较大。通过对比可以发现,在此情况下运用瞬变响应B_y(t)所得的反演结果与理论模型拟合相对较好。

图9

DOI:10.3969/j.issn.1001-1552.2003.01.002 URL [本文引用: 1]

图9 五层模型示意及结果对比

Fig.9 The diagram and inversion results comparison chart of five-layer model

Occam反演方法具有的最大特点就是它不依赖初始模型的选择,而且反演出的模型具有很好的圆滑度,符合实际要求,不会出现过度解释的现象。在反演过程中,针对灵敏度矩阵的计算,采用差分形式代替其微分形式,可使运算化简,使得计算更加便捷。本文采用线性搜索拉格朗日乘子的方法,最终求得拟合差达到精度且粗糙度达到最小的拉格朗日乘子,从而保证了计算的准确性。对典型三层模型、四层模型及五层模型进行了反演计算,从理论模型及反演结果的对比可以看出,Occam反演方法可以应用到多辐射场源地空瞬变电磁系统中,为多辐射场源地空瞬变电磁解释增加了新途径。而根据图5~图9可以看出,Occam 反演对模型具有平滑效果,在实际地层中电阻率突变的分界处,Occam反演的电阻率结果在此有一个渐变过渡,能够较好地反映出地下地质体的电性分布特征,并且对于形态简单的瞬变响应数据,其反演结果相对较好,因此应尽量运用曲线形态相对简单的瞬变响应数据进行反演计算。

3 结论

本文针对多辐射场源地空瞬变电磁法的一维反演问题进行了研究与讨论。主要解决了以下两方面的问题:第一,解决多辐射场源地空瞬变电磁法的正演问题;第二,以Occam 反演理论为基础,针对多辐射场源地空瞬变电磁法一维反演进行研究,实现对多辐射场源地空瞬变响应的反演计算;随后,本文采用多辐射场源地空系统对地下层状模型进行计算,而后对所得数据进行反演,从而得到相应的地电模型进行地球物理解释。以Occam 反演理论为基础的多辐射场源地空瞬变电磁一维反演方法还不够完善,仍存在许多问题需要继续研究和改进,本文采用线性搜索的方法来搜索合适的拉格朗日乘子,该方法虽简单,但其运算时间较长。因此可以考虑运用其他方法来解决拉格朗日乘子搜索的问题,使得运用该方法进行计算的速度得到提升,也为实现三维反演做准备。

The authors have declared that no competing interests exist.

作者已声明无竞争性利益关系。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

滕吉文

地球深部物质和能量交换的动力过程与矿产资源的形成

[J]. 大地构造与成矿学, 2003,27(1):3-21.

地球深部物质与能量的交换和深层动力过程这一科学问题的提出,乃是近年来地球科学发展与逐步向量化“进军”的必然,它是探索一系列地学前沿课题的基础。地球表面所见到的一系列地球物理场异常,地质构造格局,地球化学组分变异,无一不受到地球内部物质与能量的交换和深层动力过程的制约。如地球圈层的形成与演化,大陆伸展与裂谷,资源与能源等,均为深部物质运移和物理学、化学效应及地质构造耦合的产物。本文讨论了金属矿产资源与地壳、地幔结构及深部物质运移的动力学响应。通过几个典型矿床和其形成要素分析了成矿作用的深层动力过程。文中主要讨论了五个方面的问题:第一,金属矿产资源成矿和分布与地壳、地幔结构及深部物质运移和成矿带。第二,金的形成与深层过程。第三,金属矿床的形成、演化和分布与深部物质及能量的交换。第四,地幔热柱与成矿作用。第五,地球内部深层动力过程与流体运移和必须深化研究的几个问题。

Teng J

The dynamic process of material and energy exchange and the formation of mineral resources in deep earth

[J]. Geotectonics and Metallogenesis, 2003,27(1):3-21.

DOI:10.3321/j.issn:1002-008X.2006.05.001 URL [本文引用: 1]

[2]

邓军, 杨立强, 葛良胜 , 等.

胶东矿集区形成的构造体制研究进展

[J]. 自然科学进展, 2006,16(5):513-518.

通过系统总结胶东矿集区构造体制及其相关的华北东部中生代构造体制转折研究进展,提出了待深入研究的关键问题及对策.认为胶东金矿尽管矿化类型、产出空间、地质背景各不相同,但都是中生代一次或几次统一的重大地质事件的结果;目前对该地质事件的时间演化和空间迁移序列、尤其是地球动力学机制,仍然缺乏必要研究.定量表述各构造体制的作用方式和时空结构,区分其相对成矿贡献,仍然是胶东矿集区成矿研究中最薄弱的环节之一.认为基于区域构造解析和构造动力学的深入研究及矿化网络的全时空观测-探测,以多学科资料、现象的综合分析为约束,合理构建地质 -物理-化学模型,进行构造-热-流体-化学反应全耦合成矿动力学仿真,是深入剖析胶东矿集区构造体制时空演化轨迹和深刻理解大规模成矿作用本质的有效途径.

Deng

, Yang L

, Ge L

, et al.

Research progress on the tectonic system formed in the Jiaodong ore concentration area

[J]. Advances in Natural Science, 2006,16(5):513-518.

[3]

滕吉文, 杨立强, 刘洪臣 , 等.

岩石圈内部第二深度空间金属矿产资源形成与集聚的深层动力学响应

[J]. 地球物理学报, 2009,52(7):1734-1756.

DOI:10.3969/j.issn.0001-5733.2009.07.008 URL Magsci [本文引用: 1]

中国的快速工业化发展和经济腾飞必须有大量金属矿产资源的支撑，在共享世界资源的同时，其根本出路在于立足本土．因此提出第二深度空间（500～2000 m）金属矿产资源的找矿勘探新理念．通过对国内外金属矿产资源地球物理勘探发展的概况分析和研究提出：①金属矿产资源的集聚和分布受控于地壳内部物质与能量的交换和其深层动力过程；②在地壳内部第二深度空间存在着丰富的矿产资源，包括大型和超大型矿床；③必须充分发挥高精度地球物理探测方法的效能，并进行综合技术集成．

Teng J

, Yang L

, Liu H

, et al.

Deep dynamic response of metal mineral resources formation and accumulation in the second deep space of the lithosphere

[J]. Chinese Journal of Geophysics, 2009,52(7):1734-1756.

[4]

刘光鼎, 郝天珧, 刘伊克 .

中国大地构造宏观格架及其与矿产资源的关系

[J]. 科学通报, 1997,42(2):113-118.

<正>大地构造及其演化是地质学的理论基础,也是地球科学中对地球的结构、成分、运动和历史进行考察所取得的规律性总结,从而,它对于生产实践活动也具有重要的意义. 地质学家长期对岩石露头反复进行细致的观测,编绘出各种比例尺的地质图件,而根据不同时代的构造运动及其作用的认识与判断所编制的大地构造图件,则是从理论上对各种地质现象做出规律性的概括,进而为找矿勘探服务,也是环境保护和灾害防治的参考.

Liu G

, Hao T

, Liu Y

The macroscopic framework of China's geotectonics and its relationship with mineral resources

[J]. Chinese Science Bulletin, 1997,42(2):113-118.

[5]

Nabighian M

Electromagnetic methods in applied geophysics-theory (Volume 1)

[M]. Tulsa OK: Society of Exploration, 1988, 217-231.

[6]

Elliott

New airborne electromagnetic method provides fast deep-target data turnaround

[J].The Leading Edge,1996(4):309-310.

[7]

Smith R

, Annan

, McGowan P

A comparison of data from airborne, semi-airborne, and ground electromagnetic systems

[J]. Geophysics, 2001,66(5):1379-1385.

DOI:10.1190/1.1487084 URL [本文引用: 2]

[8]

张莹莹, 李貅, 姚伟华 , 等.

多辐射场源地空瞬变电磁法多分量全域视电阻率定义

[J]. 地球物理学报, 2015,58(8):2745-2758.

DOI:10.6038/cjg20150811 URL [本文引用: 1]

本文针对多辐射场源地空瞬变电磁法理论,建立了相应的多分量全域视电阻率定义方法.分析了利用磁场强度进行多辐射场源地空系统全域视电阻率定义的优点,针对磁场强度的各个分量,提出了各自的全域视电阻率算法,实现了多分量、全时域、全空域视电阻率计算,并分析了偏移距对全域视电阻率的影响.通过调整源的相对位置及电流方向等参数,多辐射源瞬变电磁地空系统不仅可以加强不同分量信号强度,削弱随机干扰,还可更好地分辨地下异常体的位置.通过对所设计模型的处理,证实了多辐射场源地空系统多分量全域视电阻率算法的有效性,也验证了多辐射场源地空系统的优势.

Zhang Y

, Li

, Yao W

, et al.

Determination of apparent resistivity of multi-component global TEM in multi-radiation field

[J]. Chinese Journal of Geophysics, 2015,58(8):2745-2758.

[9]

李貅, 张莹莹, 卢绪山 , 等.

电性源瞬变电磁地空逆合成孔径成像

[J]. 地球物理学报, 2015,58(1):277-288.

DOI:10.6038/cjg20150125 URL Magsci [本文引用: 1]

电性源地空瞬变电磁法具有工作效率高、勘探深度大、采集信号信噪比高、适用于地形地质条件复杂地区等优点.但是,到目前为止,由于尚未建立起该方法的解释系统,大大制约了该方法的发展.本文旨在建立起完整的地空电磁探测系统,丰富整个探测系统的理论.本文围绕地空瞬变电磁法全域视电阻率定义、瞬变电磁虚拟波场的克希霍夫偏移成像、逆合成孔径成像方法三个科学问题进行了系统研究.提出了用磁场强度定义全域视电阻率的迭代算法,理论模型试验结果表明计算出的视电阻率曲线首支趋于第一层电阻率,尾支趋于最后一层电阻率,实现了全空域、全时域视电阻率的计算;在先前研究的基础上,实现了适合电性源地空装置的瞬变电磁虚拟波场的克希霍夫偏移成像;采用相关迭加技术,实现孔径内多测点数据合成,将传统的单点处理方式发展成为逐点推移多次覆盖的逆合成孔径处理方法.层状模型试验表明:(1)全域视电阻率能够光滑、完整、渐变地反映出模型的电性信息变化;(2)当改变三层模型中间层电阻率时,全域视电阻率曲线随着参数的改变分异明显,对电性层的识别容易且直观;(3)由于在电阻率计算中同时考虑了接收机高度、偏移距、时间等各参数的影响,全域视电阻率可实现全空域、全时域的视电阻率计算.含水采空区的复杂模型算例表明:(1)根据不同测线的全域视电阻率结果可以看出,在靠近采空区的位置,全域视电阻率断面可以清晰地反映出采空区的空间位置,随着测线离采空区越来越远,采空区异常越来越弱直至消失;(2)波场变换和偏移成像的结果显示存在两个电性差异较大的界面,上界面指示地表,由于空气和大地之间的电性差异较大,故该界面波场信号反映强烈,遍布整个区域,下界面异常信号则主要集中在中部,向外逐渐减弱,指示采空区;(3)逆合成孔径成像结果表明地表界面在合成前后没有变化,而采空区异常合成后范围明显变小,且异常边界清晰,指示的采空区位置与模型吻合很好.本文借助于逆合成孔径雷达成像的基本思想,建立了一套电性源瞬变电磁地空逆合成孔径成像方法.基于反函数思想结合迭代算法提出的电性源地空瞬变电磁法的全域视电阻率定义方法,实现了全空域、全时域的视电阻率计算;借鉴瞬变电磁拟地震偏移成像算法,实现了瞬变场的三维成像;借鉴逆合成孔径雷达的思想,提出电性源瞬变电磁地空逆合成孔径算法,进一步提高了成像的分辨率.采空区模型算例表明相关叠加合成确实具有增强有用信号、提高信噪比、提高分辨率的诸多优点,证实了瞬变电磁地空逆合成孔径成像方法的有效性.

, Zhang Y

, Lu X

, et al.

Simulation of transient aperture-inverse synthetic aperture imaging of electrical source transients

[J]. Chinese Journal of Geophysics, 2015,58(1):277-288.

[10]

张莹莹, 李貅, 李佳 , 等.

多辐射场源地空瞬变电磁法快速成像方法研究

[J]. 地球物理学进展, 2016,31(2):869-876.

URL [本文引用: 1]

多辐射场源地空瞬变电磁法目前仍处于研究阶段,迫切需要完善该法的理论和解释体系,以便在生产中推广使用.本文利用地空瞬变电磁法数据量大的特点,结合微分电导与相关叠加合成算法,实现了多辐射场源地空瞬变电磁法快速成像解释.首先基于等效导电平面原理,从电偶极子源出发推导了较高精度的瞬变场快速近似正演计算公式,采用剖分叠加的思想,将其推广到多辐射场源地空瞬变电磁法空间任意一点瞬变场的快速近似正演计算,并用四种典型的层状模型验证了该法的计算精度.以此为基础,进一步推导了多辐射场源地空瞬变电磁法的视纵向电导和视深度计算公式.利用微分电导对电性界面敏感和相关叠加合成算法突出异常、压制干扰的特点,实现了多辐射场源地空瞬变电磁相关叠加快速成像.三维模型试验表明,本文提出的快速成像算法计算速度快,能有效压制噪声,且对地下界面有较好的分辨能力,有助于提高多辐射场源地空瞬变电磁法的解释精度。

Zhang Y

, Li

, et al.

Fast imaging method of ground-to-air transient electromagnetic method for multi-radiation field source

[J]. Progress in Geophysics, 2016,31(2):869-876.

DOI:10.3969/j.issn.0495-5331.2007.05.014 URL [本文引用: 1]

[11]

翁爱华

Occam反演及其在瞬变电磁测深中的应用

[J]. 地质与勘探, 2007,43(5):74-76.

Occam反演技术被应用中心回线瞬变电磁测深数据的反演中。理论模型和实际数据的反演结果表明，利用Occam反演技术可以较为准确地获得地电断面的电阻率分布。一般5次迭代目标函数就可以收敛到5×10^-2。

Weng A

Occam inversion and its application in transient electromagnetic sounding

[J]. Geology and exploration, 2007,43(5):74-76.

[12]

杨文采

非线性地震道的混沌反演——Ⅱ.关于Lyapunov指数和吸引子

[J]. 地球物理学报, 1993,36(3):376-387.

DOI: Magsci

在本文上篇《非线性地震道的混沌反演——1.理论和数值试验》中讨论了非线性反问题逐次线性化方法及迭代过程的系统特性、递推法和广义线性反演法输出波阻抗序列的特征。本篇参照混沌理论中关于Lyapunov指数、相空间和吸引子的概念,深入研究逐次线性化迭代的状态变化的阶段性及其定量描述.当引用相空间描述系统的演化时,对于随迭代变化的阻尼因子可以发现几种吸引子.最后,总结了非线性地震混沌反演的要点,指出非线性反演迭代的复杂阶段性要川混沌理论来描述,而Lyapunov指数等才是控制迭代过程的实质性参数.

Yang W

Chaotic inversion of nonlinear seismic traces—II.on Lyapunov exponents and attractors

[J]. Chinese Journal of Geophysics, 1993,36(3):376-387.

Magsci

[13]

江玉乐, 雷宛 . 地球物理数据处理教程[M]. 北京: 地质出版社, 2006.

Jiang Y

, Lei

Geophysical data processing course[M]. Beijing: Geological Publishing House, 2006.

[14]

朴化荣

. 电磁测深法原理[M]. 北京: 地质出版社, 1990:83- 101, 112-126.

[本文引用: 2]

Piao H

Principles of electromagnetic sounding[M]. Beijing: Geological Publishing House, 1990:83- 101, 112-126.

[本文引用: 2]

[15]

Kaufman A

, Keller G

. 频率域和时间域电磁测深[M]. 王建谋, 译.北京:地质出版社, 1987:153- 194, 257-279.

Kaufman A

, Keller G

Frequency and transient soundings(in Chinaese)[M].Translated by Wang J M. Beijing: Geological Publishing House, 1987:153- 194, 257-279.

[16]

李貅

. 瞬变电测深的理论与应用[M]. 西安: 陕西科学技术出版社, 2002:49- 51, 79-82.

Theory and application of transient electrical sounding[M]. Xi'an: Shaanxi Science and Technology Press, 2002:49- 51, 79-82.

[17]

方文藻, 李予国, 李貅 . 瞬变电磁测深法原理[M]. 西安: 西北大学工业出版社, 1993.

Fang W

, Li Y

, Li

Principle of transient electromagnetic sounding [M]. Xi'an: Northwestern University Industrial Press, 1993.

[18]

何继善

. 广域电磁法和伪随机信号电法[M]. 北京: 高等教育出版社, 2010: 77-88.

He J

Wide-area electromagnetic method and pseudo-random signal electrical method[M]. Beijing: Higher Education Press, 2010: 77-88.

DOI:10.3969/j.issn.1672-7940.2004.04.008 URL [本文引用: 1]

[19]

王华军

正余弦变换的数值滤波算法

[J]. 工程地球物理学报, 2004,1(4):329-335.

基于正余弦变换是汉克尔变换的特例，运用成熟的快速汉克尔变换理论，导出了正余弦变换的数值滤波算法，并给出了各250个滤波系数。以有解析解的李普希兹积分为例，与现有的算法作了试算对比，证明该算法的计算精度高，值得推广。

Wang H

Numerical screening algorithm for several cross-forming transformation

[J]. Journal of engineering geophysics, 2004,1(4):329-335.

DOI:10.3969/j.issn.0001-5733.2010.03.031 URL [本文引用: 1]

[20]

李永兴, 强健科, 汤井田 .

航空瞬变电磁一维正反演研究

[J]. 地球地球物理学报, 2010,53(3):751-759.

First we optimized the ATEM 1-D forward algorithm with new Hankel filters, then theoretically higher precision can be acquired. We developed a new ATEM 1-D inverse algorithm——alternate model adjustment method. The principle was discussed and the algorithm was described. The program was worked out, and applied for several typical models′ forward responses, the good results prove its advantages. Compared with Zohdy′s method, this method can achieve higher precision.

Li Y

, Qian J

, Tang J

Study on one-dimensional forward and inversion of aeronautical transient electromagnetic method

[J]. Chinese Journal of Geophysics, 2010,53(3):751-759.

[21]

毛立峰, 王绪本, 陈斌 .

直升机航空瞬变电磁自适应正则化一维反演方法研究

[J]. 地球物理学进展, 2011,26(1):300-305.

DOI:10.3969/j.issn.1004-2903.2011.01.035 URL Magsci

本文以中心回线式直升机航空瞬变电磁法为例对时间域航空电磁资料的一维反演算法进行理论研究.由于时间域航空电磁反演问题是典型的不适定问题,多解性严重且稳定性能较差.为此,论文使用Tikhonov正则化反演方法进行一维解释,可以得到用最光滑模型、最平坦模型和最小构造模型约束条件下的反演结果,而正则化因子是根据各次迭代的数据目标函数与约束模型目标函数的关系,以自适应的方式给出,避免了反复地搜索正则因子,使反演能够稳定快速进行.论文介绍了自适应正则化反演方法之后,用多层理论模型的垂直感生电动势响应数据进行反演试算,对比分析不同约束模型下的反演结果.结果表明,三种模型约束下,算法对初始模型要求不高,可以使用均匀半空间模型作为反演初始模型,并且均能够在不超过十次的反演迭代计算就能较好地重构理论模型,体现了良好的稳定性和可靠性.

Mao L

, Wang X

, Chen

Study on one-dimensional inversion method of helicopter aviation transient electromagnetic adaptive regularization

[J]. Progress in Geophysics, 2011,26(1):300-305.

Magsci

[22]

Constable B S

, Parker R

, Constable C

Occam’s inversion: A practical algorithm for generating smooth model from electromagnetic sounding data

[C]// Geophysics. 1987.

[23]

吴小平

大地电磁数据的Occam反演改进

[J]. 地球物理学报, 1998,41(4):547-554.

DOI: URL Magsci [本文引用: 1]

Occam 反演以其稳定的收敛性,在大地电磁（MT）数据解释中有广泛的应用。但是其每次迭代均用一维搜索方法求拉格朗日因子μ值,需要许多次正演计算,速度非常慢。在讨论了 Occam 反演中数据拟合差随μ变化的基础上,本文采用了μ值在一定步长下逐次递减的求取方法,每次迭代只需一次正演,极大地提高了计算速度。另外,反演求得光滑模型而非最光滑模型,分辨率更高。理论及实际数据的反演试算均表明,和 Occam 反演相比,反演依然稳定,但速度更快,结果也更真实。

Wu X

Improvement of Occam inversion of magnetotelluric data

[J]. Chinese Journal of Geophysics, 1998,41(4):547-554.

[24]

周道卿, 谭捍东, 王卫平 .

频率域航空电磁资料Occam反演研究

[J]. 物探与化探, 2006,30(2):162-165.

DOI:10.3969/j.issn.1000-8918.2006.02.016 URL Magsci [本文引用: 1]

通过对频率域航空电磁测量常用装置类型电磁响应分析和Occam反演方法介绍,成功地将Occam思想运用到频率域航空电磁资料的反演解释工作中。Occam反演是一种正则化的反演方法,它在追求模拟数据与原始测量曲线最大拟合的同时,要求模型数据最平滑或最圆滑,因而受初始模型影响小,能够达到稳定收敛。经过对理论模型和实测数据试算表明,该方法运算稳定,收敛速度快,是一种有效的航空电磁数据处理解释方法。

Zhou D

, Tan H

, Wang W

Occam inversion of aeronautical electromagnetic data in frequency domain

[J]. Geophysical and Geochemical Exploration, 2006,30(2):162-165.

[25]

何梅兴, 胡祥云, 叶益信 , 等.

2.5维可控源音频大地电磁法Occam反演理论及应用

[J]. 地球物理学进展, 2011,26(6):2163-2170.

DOI:10.3969/j.issn.1004-2903.2011.06.033 URL Magsci [本文引用: 1]

Occam反演是一种正则化的光滑模型反演方法,它在寻找最小拟合差的同时追求最光滑模型.因地下介质的电性通常是连续变化的,Occam反演考虑了地下介质横向和纵向的光滑情况,应用Occam反演避免了模型电性参数被随意间断.通过建立层状和存在两个异常体的模型,利用2.5维可控源音频大地电磁有限元法作正演响应计算,对响应数据基于Occam反演理论拟合反演.结果表明,Occam反演对可控源音频大地电磁法响应数据反演稳定收敛,理论模型与反演结果的视电阻率曲线形态基本一致,测点数据的拟合误差小于3%,反演结果能反映出模型电阻率的光滑效果.

He M

, Hu X

, Ye Y

, et al.

Occam inversion theory and application of 2.5D controllable source audio magnetotelluric method

[J]. Progress in Geophysics, 2011,26(6):2163-2170.

[26]

刘方镝

海洋可控源电磁法一维正演与Occam反演研究

[D]. 成都:成都理工大学, 2012.

Liu F

1D direct problem and Occam inversion of ocean controlled source electromagnetic method

[D]. Chengdu: Chengdu University of Technology, 2012.

DOI:10.3969/j.issn.1001-1749.2015.06.03 URL [本文引用: 1]

[27]

张君涛, 周军, 王绪本 , 等.

一维大地电磁Occam反演拉格朗日乘子的搜索

[J]. 物探化探计算技术,2015(6):687-692.

在大地电磁反演中,Occam法因其在反演稳定性和模型分辨率等方面的优势,得到广泛应用。但由于其每次迭代都需要不断地搜索拉格朗日乘子,因而拉格朗日乘子的搜索效率对Occam法反演的运算速度起着至关重要的作用。为提高拉格朗日乘子的搜索效率,这里提出将拉格朗日乘子的搜索从自然数域中转到以10为底的对数域中进行,同时在区间最小值的搜索中采用二次函数极小值搜索法。通过大量一维大地电磁反演验证,该方法将拉格朗日乘子的搜索效率提高了20%~50%,大大提高了Occam反演的运算速度,方便了其在高维反演中的应用。

Zhang J

, Zhou

, Wang X

, et al.

Search of one-dimensional magnetotelluric Occam inversion of Lagrange multiplier

[J]. Computing Techniques For Geophysical and Geochemical Exploration, 2015,37(6):687-692.

[28]

陈小斌, 赵国泽, 汤吉 , 等.

大地电磁自适应正则化反演算法

[J]. 地球物理学报, 2005,48(4):937-946.

DOI:10.3321/j.issn:0001-5733.2005.04.029 URL Magsci [本文引用: 1]

针对大地电磁正则化反演中正则化因子的选取困难问题提出了自适应正则化反演算法（Adaptive Regularized Inversion Algorithm, ARIA）. 在该算法中， ①提出了一种新的数据方差处理方法：数据方差规范化，使得数据方差的大小只对数据的拟合发生影响，不对数据目标函数和模型约束目标函数的权重产生影响，从而减少了正则化因子取值的影响因素；②提出了粗糙度核矩阵的概念，并给出了由基本结构插值基函数计算粗糙度核矩阵的公式，使得模型目标函数的构建更为简便、直接；③根据数据目标函数、模型约束目标函数和正则化因子之间的关系，提出了两种正则化因子自适应调节方法. 本文详细阐述了最平缓模型约束下的大地电磁一维连续介质反演的ARIA实现，以几个算例的分析比较来说明ARIA的有效性.

Chen X

, Zhao G

, Tang

, et al.

Magneto-electromagnetic adaptive regularization inversion algorithm

[J]. Chinese Journal of Geophysics, 2005,48(4):937-946.