基于MATLAB的大地电磁测深正_省略_质的正演和阻尼最小二乘法反演为例_李斌

　2010年第9期　　　　　　　　　　　　内蒙古石油化工

35

基于MATLAB的大地电磁测深正反演实现

——以层状一维介质的正演和阻尼最小二乘法反演为例

李　斌,郭嵩巍,郑　凯,刘　云

(成都理工大学地球探测与信息技术重点实验室,四川成都　610059)

　　摘　要:本文扼要的叙述了大地电磁的基本原理,简要介绍了大地电磁的正反演方法和流程,重点阐述了基于MATLAB软件平台的实现,并给出了此种方法可扩展的几个方面。

关键词:大地电磁测深;正演;反演;MATLAB;绘图;接口;并行运算

+

　　中图分类号:P631.325　　文献标识码:A　　文章编号:1006—7981(2010)09—0035—021　大地电磁方法简介1.1　大地电磁测深理论概述

大地电磁测深是研究地壳和上地幔构造的一种地球物理探测方法。它以天然交变电磁场为场源,当交变电磁场以波的形式在地下介质中传播时,由于电磁感应作用,地面电磁场的观测值将包含有地下介质电阻率的分布信息。而且,由于电磁场的集肤效应,不同周期的电磁场信号具有不同的穿透深度[1],因此,在地面上观测大地电磁场,可获得地下不同深度介质电阻率分布的信息。1.2　大地电磁测深理论的几点假设和论证

吉洪诺夫院士和卡尼尔研究员提出了假设并论证了以下几点:①将场源近似地看为平面电磁波垂直入射大地。②引入波阻抗的概念(Z=E/H),表征地球电性分布对大地电磁场的响应。③利用单点大地电磁场观测研究地球电性分布是可能的。具体可见参考书目[2]。1.3　平面电磁波在均匀大地介质中的传播

电磁波在无耗介质中传播严格遵循麦克斯韦方程,但在有耗介质中是否遵循麦克斯韦方程有学者提出异议,但据参考文献[3]我们仍然采用麦克斯韦方程。

D/ t=crotH=j+ rotE=- B/ tdivB=0divD=qrotH=j=σErotE=iωμHdivH=0divE=0

→

→→→

→

→

→

→

→→→

→

→

→

→

→

→

→

单色波(单频)方程为赫姆霍茨方程。

55

→

22→→

E-k2E=0H-k2H=0-iωμ/ρ

→

→

→

(11)(12)(13)(14)

k=

E E

=→

→

H H=(15)Ez=0(16)

(17)Hz=0

这是一个二阶常微分方程,其一般解为:

-kz

+Bekz(18)x(z)=AeE

Hy(z)=-(Ae-kz-Bekz(19)iωμ

其中A和B是由边界条件确定的积分常数。边界条件为:①B和D在分界面两侧法向量连和H在分界面两侧切向量连续。续;②E2　大地电磁正演理论(以一维水平层状各向同性介质为例)

引入波阻抗的定义:Z=E/H对于第m层的波阻抗有

-kzkz

Ex(zm)memm+BemmAZ=-(20)m=-kzkz

Hy(zkmkm(Am)memm-Bemm)E(z)iωμAe-kmzm+1+Bekmzm+1

Z=-m+1=Hy(z)kmkm(Ame-kmzm+1-Bekmzm+1)

(21)

波阻抗在分界面上是连续的,任一层底界面的波阻抗等于其下层相邻介质顶界面的波阻抗,因而可转换为求解同一层顶面和底面波阻抗之间的关系。假定第m层顶面深度为Zm,底面深度为Zm+1,联立以上两式可得递推公式

-2kmhm

1-Lm+1e

Z(22)m=Zom-2kmhm

1+Lm+1eZom-Zm+1

L(23)m+1=Zom+Zm+1

iωμ

(24)n=Zom=-Zkn

Z(25)m+1-Zm=hm(m=1,2,…,n-1)

12|1.n|(26)T=ρωμZ

→

→

→

→

(1)(2)(3)(4)(5)(6)(7)(8)

rri

→→

rotA=5×A=(9)AxAyAz

→

Ax A Azy

++(10)divA=5gA=在大地电磁测深勘探问题中,研究的电磁场是随时间变化很慢,对时间的二阶导数与一阶导数相比可以忽略不计,称为场的似稳模型,那么似稳场的

:

36

　　如下给出一维正演的流程图

:

内蒙古石油化工　　　　　　　　　　2010年第9期　

本文未给出程序具体实现代码,但根据流程图可自己编制出代码来验证。

图1　一维正演的程序流程图

3　大地电磁反演理论

阻尼最小二乘法也称之为麦夸特法,他将高斯

=G修改为[PTP+αI]-牛顿法的法方程组PTPΔλ

=G,式中α称为阻尼因子(α>0),它用于控制校Δλ

正方向和步长。I是单位矩阵。可以看出,当α=0时,麦夸特法就变为高斯-牛顿法;当α→∞时,麦夸特法就退化为梯度法。适当的选择阻尼因子,可以兼具高斯-牛顿法和最速下降法的优点,同时克服二者的缺点。

大地电磁测深反演用麦夸特法时,目标函数为:

ρ-ρ2

)ρi=1Ti

逼近目标函数的二次函数取:Χ=

∑

N

((27)

m

ρci(k)2

Χ=Δλ]+αi

i=11 j=1λi

2(k)(Δλ(28)j)

T

由极值必要条件5Χ=0,求得法方程组:[JJ

(k)

+α(29)I]Δλ=S

(k)T-1

+αI]G其中:Δλ=[PP

1 1 ρ1 ρc1c1c1

…

T11T12T1m1 ρ1 ρ1 ρc2c1c1

…

T21T22T2mJ=

…………ρρρcNcNcN

…(k)λλρρρTN 1TN 2TN λ=λ

N

1 ρc1(k)

)(ρTi-ρc1)2(1Tii=1ρ

∑

N

(K)Ti-ρci)2[(ρj=ρTi

∑

m

∑

∑∑

i=N

S=

1 ρc1(k)

)(ρTi-ρc1)2(Ti21

…

1 ρc1(k)

)

(ρTi-ρc1)(mρTi(k=0,1,2,…)i=1

如下给出一维反演的流程图(如图2):4　MATLAB简介以及实现的正反演结果

大地电磁测深数据的正反演计算很早就实现了计算机程序自动完成,但大多数地球物理工作者都是采用FORTRAN,C,C++等语言编写源程序,然后编译连接成可执行程序,有学者将此类语言称为第二代计算机语言;而本文采用了MATLAB(矩阵实验室)环境中的M语言编写源程序,有学者将M语言称为第三代计算机语言,因为其运行是翻译源代码,而非编译连接.以下是实现的正演曲线和反演曲线图件,从反演结果来看,当迭代若干次后和正,,

∑

N

5　可扩展的几个方面5.1　自MATLABR2009a版本开始将提供并行运算工具箱,在原有M语言的基础上仅做稍许改动甚至不做改动,便可完成并行程序的翻译,最多可支持8核的CPU,这对计算量非常大的地球物理数据处理来说是意义非常。5.2　MATLAB支持与外设的串口通信,这对于要求实时处理资料的地球物理方法来说大有用处。6　总结

MATLAB提供了大量成熟可靠的数值计算方法,以为国内外许多科研单位所接受和应用,而且MATLAB可以方便的实现数据可视化和图形绘制,特别是三维图形绘制,如果想脱离MATLAB环境运行,也有很多种方法可以实现,比如自带的MCC编译器、COM等方法,这样利用MFC或者VB等编辑界面,利用MATLAB实现正反演计算和绘图并将之生成DLL文件或COM组件供MFC调用,将会大大提高编程效率,缩短研发周期,所以这也正是本文的意义所在。除了大地电磁测深一维正反演外,作者还完成了二维正反演模块,限于篇幅,考虑将在另外文章中发表。

[参考文献]

[1]　石应骏,等.大地电磁测深法教程理论[M].北

京:地质出版社,1985.[2]　米萨克N纳比吉安.勘察地球物理电磁法(第

一卷)理论[M].赵经祥,译.北京:地质出版社,1992.

[3]　王长清等Maxwell方程用于电磁脉冲在损耗

介质中的传播问题[J].电波科学学报,1999,