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

∑
N
(
( 27)
m
i = 1 N
ρ ci ( k) 2 Χ= Δλ ]+ α i i= 1 1 λ j= 1 i 2 ( k) (Δλ ( 28) j ) T 由极值必要条件 5 Χ = 0, 求得法方程组 : [ J J ( k) + α = S ( 29) I] Δλ ( k) = [ PT P + α I ]- 1 G 其中: Δλ 1 ρ 1 ρ 1 ρ c1 c1 c1 … ρ ρ ρ T1 λ 1 T1 λ 2 T1 λ m 1 ρ 1 ρ 1 ρ c2 c1 c1 … ρ ρ T2 λ 1 T2 λ 2 T2 λ m J = ρ … … … … cN cN cN 1 ρ 1 ρ 1 ρ … ( k) ρ ρ ρ TN λ 1 TN λ 2 TN λ m λ =λ N 1 ρ c1 ( k) ) (ρ Ti - ρ c1 ) 2( λ 1 Ti i = 1 ρ
∑
1 ( K) Ti - ρ ci ) 2 [ (ρ j= ρ Ti
∑
m
∑
∑ ∑
…
N
S=
i = 1
1 ρ c1 ( k) ) (ρ Ti - ρ c1 ) 2( ρ λ Ti 2
1 ρ c1 ( k) ) (ρ Ti - ρ c1 ) 2( λ m ρ Ti ( k= 0, 1, 2, … ) 如下给出一维反演的流程图 (如图 2): 4 M A T L A B简介以及实现的正反演结果 大地电磁测深数据的正反演计算很早就实现了 计算机程序自动完成 , 但大多数地球物理工作者都 是采用 FO R T RAN , C, C + + 等语言编写源程序 , 然 后编译连接成可执行程序 , 有学者将此类语言称为 第二代计算机语言 ; 而本文采用了 M A T L A B (矩阵 实验室 )环境中的 M 语言编写源程序 , 有学者将 M 语言称为第三代计算机语言 , 因为其运行是翻译源 代码 , 而非编译连接 . 以下是实现的正演曲线和反 演曲线图件 , 从反演结果来看 ,当迭代若干次后和正 演曲线基本重合 , 验证了算法的正确性。 考虑篇幅 ,
ro t A= 5 × A=
36
如下给出一维正演的流程图 :
内蒙古石油化工 2010年第 9期
本文未给出程序具体实现代码 , 但根据流程图可自 己编制出代码来验证 。
图 1 一维正演的程序流程图 3 大地电磁反演理论 阻尼最小二乘法也称之为麦夸特法 , 他将高斯 = G 修改为 [ PT P + αI ] - 牛顿法的法方程组 PT P Δλ = G ,式中 α 称为阻尼因子 (α > 0) , 它用于控制校 Δλ 正方向和步长。 I 是单位矩阵。可以看出 ,当 α = 0时 , 麦夸特法就变为高斯 - 牛顿法 ; 当 α →∞ 时 , 麦夸特 法就退化为梯度法 。 适当的选择阻尼因子 , 可以兼具 高斯 - 牛顿法和最速下降法的优点 , 同时克服二者 的缺点 。 大地电磁测深反演用麦夸特法时 , 目标函数为: ρ Ti - ρ ci 2 ) ρ Ti 逼近目标函数的二次函数取 : Χ=
∑
N
i = 1
5 可扩展的几个方面 5. 1 自 M A T L A B R2009a 版本开始将提供并行运 算 工具箱 , 在原有 M 语言的基础上仅做稍许改动甚 至不做改动 , 便可完成并行程序的翻译 , 最多可支 持 8 核的 C PU , 这对计算量非常大的地球物理数据 处理来说是意义非常 。 5. 2 M A T L A B 支持与外设的串口通信 , 这对于要 求实时处理资料的地球物理方法来说大有用处。 6 总结 M A T L AB 提供了大量成熟可靠的数值计算方 法 , 以为国内外许多科研单位所接受和应用 , 而且 M A T L A B 可以方便 的实 现数据 可视化 和图形 绘 制 ,特别是三维图形绘制 , 如果想脱离 M A T L AB 环 境运 行 , 也有 很 多种 方 法可 以 实现 , 比 如自 带 的 M CC 编译器 、 CO M 等方法 , 这样利用 M FC 或者 V B 等编辑界面 , 利用 M A T L AB 实现正反演计算和绘 图并将之生成 D L L文件或 CO M 组件供 M FC 调用 , 将会大大提高编程效率 , 缩短研发周期 , 所以这也正 是本文的意义所在。 除了大地电磁测深一维正反演 外 , 作者还完成了二维正反演模块 , 限于篇幅 , 考虑 将在另外文章中发表 。 [ 参考文献 ] [1 ] 石应骏 ,等 . 大地电磁测深法教程理论 [ M ]. 北 京: 地质出版社 , 1985. [2 ] 米萨克 N 纳比吉安 . 勘察地球物理电磁法 (第 一卷 )理论 [ M ]. 赵经祥 , 译 . 北京 : 地质出 版 社 , 1992. [ 3] 王长清等 Max well 方程用于电磁脉冲在损耗 介质中的传播问题 [ J ]. 电波科学学报 , 1999, 14( 1) .
2010年第 9期 内蒙古石油化工
35
基于 M A T L AB的大地电磁测深正反演实现
— — 以层状一维介质的正演和阻尼最小二乘法反演为例
李 斌 , 郭嵩巍 , 郑 凯 , 刘 云
( 成都理工大学 地球探测与信息技术重点实验室 , 四川 成都 610059) 摘 要: 本文扼要的叙述了大地电磁的基本原理 , 简要介绍了大地电磁的正反演方法和流程 , 重点 阐述了基于 M A T L AB 软件平台的实现 ,并给出了此种方法可扩展的几个方面。 关键词 : 大地电磁测深 ; 正演 ; 反演 ; M A T L A B; 绘图 ; 接口 ; 并行运算 + 中图分类号 : P631. 3 25 文献标识码: A 文章编号: 1006 — 7981( 2010) 09 — 0035 — 02 1 大地电磁方法简介 1. 1 大地电磁测深理论概述 大地电磁测深是研究地壳和上地幔构造的一种 地球物理探测方法 。 它以天然交变电磁场为场源 , 当 交变电磁场以波的形式在地下介质中传播时 , 由于 电磁感应作用 , 地面电磁场的观测值将包含有地下 介质电阻率的分布信息 。 而且 ,由于电磁场的集肤效 应 , 不同周期的电磁场信号具有不同的穿透深度 [ 1] , 因此 ,在地面上观测大地电磁场 , 可获得地下不同深 度介质电阻率分布的信息。 1. 2 大地电磁测深理论的几点假设和论证 吉洪诺夫院士和卡尼尔研究员提出了假设并论 证了以下几点 : ① 将场源近似地看为平面电磁波垂 直入射大地。 ②引入波阻抗的概念 ( Z = E / H) , 表征 地球电性分布对大地电磁场的响应。 ③利用单点大 地电磁场观测研究地球电性分布是可能的 。 具体可 见参考书目 [ 2 ]。 1. 3 平面电磁波在均匀大地介质中的传播 电磁波在无耗介质中传播严格遵循麦克斯韦方 程 , 但在有耗介质中是否遵循麦克斯韦方程有学者 提出异议 , 但据参考文献 [ 3 ]我们仍然采用麦克斯韦 方程。
→ → → → → → → → →
单色波 (单频 )方程为赫姆霍茨方程 。 5 5
→ 2 2 → → 2
E- k2 E= 0
→
→
( 11) ( 12) ( 13) ( 14)
H- k H= 0 - i ω μ/ ρ
→
k= E = x
→
E y
→
H H = ( 15) x y Ez= 0 ( 16) ( 17) Hz= 0 这是一个二阶常微分方程 , 其一般解为: - kz + Bekz ( 18) x ( z) = Ae E k Hy ( z ) = ( Ae- kz - Bekz ( 19) i ω μ 其中 A 和 B 是由边界条件确定的积分常数 。 边界条件 为: ① B 和 D 在分界面两侧法向 量连 和 H 在分界面两侧切向量连续 。 续 ;② E 2 大地电磁正演理论 (以一维水平层状各向同性介 质为例 ) 引入波阻抗的定义: Z= E / H对于第 m 层的波 阻抗有 - k z k z E x ( z m ) m e m m+ Be m m i ω μ A Z = ( 20) m= - km zm k m zm ) k ( A e Be ) Hy ( z k m m m m - k z k z Ex ( z) i ω μ Am e m m+ 1 + Be m m+ 1 Z = m+ 1 = Hy ( z) km km ( Am e- k mzm+ 1 - Bek m zm+ 1 ) ( 21) 波阻抗在分界面上是连续的 , 任一层底界面的 波阻抗等于其下层相邻介质顶界面的波阻抗 , 因而 可转换为求解同一层顶面和底面波阻抗 之间的关 系。 假定第 m 层顶面深度为 Z m , 底面深度为 Z m+ 1 , 联 立以上两式可得递推公式 - 2k m h m 1- Lm + 1e Z ( 22) m= Z om - 2k m h m 1 + Lm e + 1 Zom - Z m+ 1 L ( 23) m+ 1 = Zom+ Z m+ 1 i ω μ ( 24) n= Z om = Z kn Z ( 25) m+ 1 - Zm = hm ( m= 1, 2,… , n- 1) 1 2 | 1. n| ( 26) T= ρ ω μ Z
→ →
r x Ax
r
i
( 9) y z Ay Az → Ax A Az y + + ( 10) div A= 5 g A= x y z 在大地电磁测深勘探问题中 , 研究的电磁场是 随时间变化很慢 , 对时间的二阶导数与一阶导数相 比可以忽略不计 , 称为场的似稳模型 , 那么似稳场的 收稿日期 : 2009- 12- 04
→ → → →
ro t H= j+ ro t E= div B= 0
→
D / t= c B/ t
( 1) ( 2) ( 3) ( 4)
div D= q ro t H= j= σ E rBiblioteka Baidu t E= i ω μH
→ → → → → →
( 5) ( 6) ( 7) ( 8)
→
div H= 0
→
div E= 0