起伏地形大地电磁二维反演

㊃ ５８８㊃
物㊀ 探㊀ 与㊀ 化㊀ 探
４０ 卷 ㊀
果，甚至不能收敛 ［１７］ ㊂ 大地电磁正则化反演的传统 做法是用光滑函数表示模型参数， 通过求解最小范 但对于地质体界面分辨率较低㊂ 基于此，Ｐｏｒｔｎｉａｇｕ⁃ ｉｎｅ 和 Ｚｈｄａｎｏｖ
［１８］
百度文库
数或者光滑稳定的函数，即可得到光滑的反演模型， 于重力反演中； 随后，Ｍｅｈａｎｅｅ 和 Ｚｈｄａｎｏｖ 提出了聚焦反演方法，并成功应用
（ １． 桂林理工大学 地球科学学院， 广西 桂林 ㊀ ５４１００６； ２． 美国 犹他大学 矿业与地球科学学院， 犹 他州 盐湖城㊀ ８４１１２）
摘 要： 为了适应地形起伏的实际地质情况，开展带地形的最小二乘二维反演研究㊂ 鉴于大地电磁 （ ＭＴ） 反演的不 适定问题，引入 Ｔｉｋｈｏｎｏｖ 的正则化方法，从而获得关于总目标函数的方程，利用光滑约束最小二乘法求解总目标函 数方程㊂ 由于正则化因子值与反演精度以及稳定性相关， 采用主动约束平衡方法获取最优化的正则化因子， 以确 保反演精度和稳定性都达到最佳㊂ 与此同时，利用电磁场互易定理以节省反演迭代过程求解雅可比矩阵的计算时 间㊂ 构建了若干地质构造模型进行试算，分别讨论 ＴＥ㊁ＴＭ 模式以及二者联合模式的反演结果，并与前人研究工作 对比以说明本文方法的反演效果㊂ 关键词： 起伏地形；大地电磁；正则化；最小二乘；反演 中图分类号： Ｐ６３１㊀ ㊀ ㊀ 文献标识码： Ａ㊀ ㊀ ㊀ 文章编号： １０００－８９１８（ ２０１６） ０３－０５８７－０７
Ｘｉｏｎｇ Ｂ，Ｌｕｏ Ｔ Ｙ，Ｃａｉ Ｈ Ｚ，ｅｔ ａｌ． Ｔｗｏ⁃ｄｉｍｅｎｓｉｏｎａｌ ｍａｇｎｅｔｏｔｅｌｌｕｒｉｃ ｉｎｖｅｒｓｉｏｎ ｏｆ ｔｏｐｏｇｒａｐｈｙ［ Ｊ］ ．Ｇｅｏｐｈｙｓｉｃａｌ ａｎｄ Ｇｅｏｃｈｅｍｉｃａｌ Ｅｘｐｌｏｒａｔｉｏｎ，２０１６，４０（ ３） ：
起伏地形大地电磁二维反演
ʏ Ｌ（ Ｖ，������ Ｖ，������ Ｖ） ｄｘｄｚ，
Ω ｘ ｚ
（１）
针对大地电磁反演的不适定性， 引入 Ｔｉｋｈｏｎｏｖ 的正则化方法 Ｐ α（ ｍ） ＝ φ（ ｍ） ＋ αｓ（ ｍ） ， （８） ａ 其中，Ｐ （ ｍ） 为总目标函数； φ （ ｍ ） 为观测数据与预 测数据之差的平方和， α 为正则化因子， ｓ （ ｍ ） 为模 型约束目标函数㊂ 根据不同的先验约束条件， 模型 约束目标函数的构建方法有多种， 常用的 ３ 种是最 小模型约束㊁最平缓模型约束㊁ 最光滑模型约束 ［２５］ ， 文中采用基于先验模型的最光滑模型约束㊂ 因此， 反演问题的总目标函数可表示为 ［２６］ Ｐ α（ ｍ） ＝ Ｗ ｄ ［ ｄ ｏｂｓ － Ｆ （ ｍ） ］ ２ ＋ α Ｗ ｍ（ ｍ － ｍ ｒｅｆ ） ２ ， （９） 式中，Ｗ ｄ 为观测数据权系数矩阵， ｄ ｏｂｓ 为观测数据， Ｆ 为正演响应函数， Ｗ ｍ 为光滑度矩阵， ｍ 为模型参 数，ｍ ｒｅｆ 为先验模型㊂ 经过适当的运算， 可得到模型 修改量 Δｍ 的表达式 Ｔ Ｔ Ｔ Δｍ ＝ （ Ｊ ｋ Ｗ ｄ Ｗ ｄ Ｊ ｋ ＋ αＷ ｍ Ｗ ｍ ） － １ ㊃
， ＶＮ ） Ｔ ， Ｖｎ
ｓ ｓ
ｘ
＝
１㊀ 有限元正演
延伸，假设电㊁磁场 Ｖ（ ｘ， ｚ ） 在求解区域 Ω 内满足微 分方程 ＤＶ ＝ ｆ， 在 Ω 边界 ������Ω 上满足给定的边界条 件 Ｖ ＝ ｆ ０ ；Ｄ 为微分算符， ｆ 为关于场源的函数， ｆ ０ 为 已知函数，根据变分法可知，求解偏微分方程相当于 数 ｌ 表示为 寻找使得拉格朗日函数 ｌ 最小的场值 Ｖ，拉格朗日函 ｌ＝ 方程
熊彬１ ，罗天涯１ ，蔡红柱２ ，刘云龙１ ，吴延强１ ，郭胜男１
㊀ ㊀ 大地电磁测深法以天然交变电磁场为场源， 当 交变电磁场以波的形式在地下介质中传播时， 由于 电磁感应作用，地面电磁场的观测值将包含地下介 地对天然电磁场的频率响应， 可获得地下不同深度 介质电阻率信息
［１］
（ ｒａｐｉｄ ｒｅｌａｘａｔｉｏｎ ｉｎｖｅｒｓｉｏｎ，ＲＲＩ） 等㊂ ＯＣＣＡＭ 二维反 演算法是 ｄｅＧｒｏｏｔ⁃Ｈｅｄｌｉｎ 等 ［１３］ 由 ＯＣＣＡＭ 一维反演 在反演计算中，每一次迭代都将求解模型参数的偏 网格剖分的加密而导致计算时间变长；ＲＥＢＯＣＣ 是 Ｓｉｒｉｐｕｎｖａｒａｐｏｒｎ 等 ［１４］ 在 ＯＣＣＡＭ 二维反演方法的基 础上改进得来，采用模型参数的协方差矩阵将反演 问题从 Ｍ 维模型参数空间转变为 Ｎ 维数据空间，由 于 Ｎ≪Ｍ， 使得反演计算时间大大缩减； ＳＢＩ （ ｓｈａｒｐ ｂｏｕｎｄａｒｙ ｉｎｖｅｒｓｉｏｎ） 反演方法是 ｄｅＧｒｏｏｔ⁃Ｈｅｄｌｉｎ 等 ［１５］ 在原有的 ＯＣＣＡＭ 反演算法的基础上， 构建新的模 型粗糙度矩阵，提出一种针对陡边界条件的二维反 演算法；为了避免 ＯＣＣＡＭ 法直接线性搜索， Ｓｍｉｔｈ 和 Ｂｏｏｋｅｒ ［１６］ 提出 ＲＲＩ 方法，通过分析每个测深点之 下的电阻率变化，将二维反演问题降为一维反演问 题，ＲＲＩ 法不直接求雅可比矩阵，而是对正演求得的 电场值做积分运算， 获得视电阻率对模型参数的偏 但在反演时若某些参数控制不好， 便得不到理想结 导数，每次迭代只要做一次正演， 提高了计算速度， 方法扩展而来的，其中正演模拟方法为有限单元法， 导数矩阵，因此，反演过程会随着测点数的增加以及
Ｔ Ｔ Ｔ ［ Ｊ ｋ Ｗ ｄ Ｗ ｄ（ ｄ ｏｂｓ － ｄ ｋ ） ＋ αＷ ｍ Ｗ ｍ（ ｍ ｒｅｆ － ｍ ｋ ） ］ ，
式中，η 和 γ 的含义见文献 ［ ２１］ ㊂ 将式 （ ３） 代入式 （２） 即可获得电㊁ 磁场满足的具体微分表达式㊂ 加 入给定的边界条件获得边值问题以及用有限单元法 求解边值问题的详细过程参见文献［２１ － ２４］ ㊂ 在获得 Ｖ 之后，为了求解 ＴＭ㊁ＴＥ 模式的阻抗
的正则化方法应用到最小二乘优化方法中；同时，用 电磁场互易定理来计算雅可比矩阵，最后，对若干典
������Ｖ Ｎ ö æ ������Ｖ１ ������Ｖ２ ｓ ç ， ， ， ｘ ÷ ，Ｋ 是包含地表节点单元的系 ������ｎ ������ｎ ������ｎ è ø ｓ 数矩阵，Ｖ ｓ 是地表节点的场值， Ｖ ｎ 是 Ｖ ｓ 的法向导 数，Ｎ ｘ 是 ｘ 方向上地表节点数， 从而可以求出场值 的法向导数 ������Ｖ ／ ������ｎ㊂ 切向导数 ������Ｖ ／ ������ｔ 可利用地表沿 ｘ 方向相邻节点的差商求取㊂ 将 ������Ｖ ／ ������ｎ 和 ������Ｖ ／ ������ｔ 代入 式（６） 即可解出 ������Ｖ ／ ������ｚ 并代入式 （ ５） 可求出辅助场 Ｉ，进而可通过式（４） 获得地下阻抗信息㊂
［１９］
Ｈｙ 还需求出辅助场 Ｉ
ＴＭ
Ｚ ｘｙ ＝
Ｅｘ
＝
Ｅｙ Ｉ Ｖ ＝－ ， ，㊀ Ｚ ｙｘ ＝ Ｖ Ｈｘ Ｉ ������ ｚ Ｖ ＝ － τＩ，
ＴＥ
（４） ， （５）
反演方法引入大地电磁二维反演问题的求解， 且取 得很高的反演精度；ＮＬＣＧ 反演算法是 Ｒｏｄｉ 等
［２０］
将聚焦 提
式中，Ｉ ＝ Ｅ ｘ 或 Ｉ ＝ － Ｈ ｘ ㊂ 沿用 Ｌｅｅ 等
质电阻率分布信息，根据电磁场的集肤效应，研究大 是不可或缺的一部分㊂ 随着正演理论方法的发展， 反演也不断取得新的进展， 经历吉洪诺夫和卡尼亚 演
［２ ５］
㊂ 在大地电磁数据解释中， 反演
时期的一维反演， 到 ２０ 世纪 ７０ 年代开始的二维反 主流㊂ 然而三维反演非常耗时， 对计算机内存要求 比较高，此外，大地电磁野外数据采集周期往往比较 ，至今， 三维地质结构的反演
２㊀ 计算辅助场
（１０）
㊀ ３期
熊彬等：起伏地形大地电磁二维反演
㊃ ５８９㊃
３． ２㊀ 雅可比矩阵的计算
式中，Ｊ 为雅克比矩阵， ｍ ｋ 为模型的第 ｋ 次迭代值， ｄ ｋ ＝ Ｆ（ ｍ ｋ ） ㊂ 题，文中利用电磁场的互换定理 ［２７］ 来求解雅可比矩 阵㊂ 此外，反演过程中，针对电阻率和视电阻率很大 模型参数表示为对数形式 ［２４］ ： ｍ ＝ （ ｌｏｇρ １ ，ｌｏｇρ ２ ， ，ｌｏｇρ Ｍ ） Ｔ ， （１１） 的动态范围和便于在目标函数中采用均方误差， 将 同时也将视电阻率转换成其对数形式， 在二维地质 考虑到雅可比矩阵求解需要耗费大量时间的问
在直角坐标系中， 令地质构造走向沿 ｙ 轴无限
３㊀ 最小二乘光滑约束反演
３． １㊀ 基本公式
Ｌ 为拉格朗日密度㊂ 当场值 Ｖ 满足 Ｅｕｌｅｒ⁃Ｌａｇｒａｎｇｅ ������Ｌ é ������Ｌ ù é ������Ｌ ù ������ ｚ ê （２） ê ������（ ������ Ｖ） ú ú ＋ ������ ｘ ê ê ������（ ������ Ｖ） ú ú ＝ ������Ｖ ë û ë û ｚ ｘ 时，拉格朗日函数取得最小值，拉格朗日密度 Ｌ 可表 示为 Ｌ＝ １ １ γ （ ������ ｚ Ｖ） ２ ＋ （ ������ ｘ Ｖ） ２ ＋ Ｖ２ ， ２η ２η ２ （３）
［６ １２］
成为研究的
长，勘探成本很高， 因而其在实际应用中并未普及㊂ 而对一些地质条件，二维结构的假设也是合理的，且 二维反演具有计算速度快㊁数据存储量小等优势，因 此，深入研究 ＭＴ 二维正反演仍不失为减小勘探成 本的有效手段㊂ 演㊁ ＲＥＢＯＣＣ （ Ｒｅｄｕｃｅｄ ｂａｓｉｓ ＯＣＣＡＭ ） 反 演㊁ ＮＬＣＧ
ｓ ｓ ｓ ｓ ＝ （ Ｖ１ ， Ｖ２ ， ｓ Ｔ
的算法， 地
出的一种快速㊁稳定㊁ 收敛的二维反演方法， 其正演 模拟方法是基于有限差分法， 反演结果具有较高的 分辨率， 计算所需内存小， 另外，ＮＬＣＧ 反演方法可 以进行联合反演，如 ＴＥ 和 ＴＭ 视电阻率和阻抗相位 的联合㊂ 然而，ＮＬＣＧ 对初始模型的选择存在依赖 性，正则化参数选择需要凭经验输入㊂ 程稳定性以及偏导数计算等问题， 将求解病态问题 型模型进行计算分析㊂ 在前人研究的基础上，针对反演分辨率㊁ 反演过
㊀ ㊀ 收稿日期： ２０１５－１１－０２
常用的大 地 电 磁 二 维 反 演 方 法 有 ＯＣＣＡＭ 反
（ ｎｏｎｌｉｎｅａｒ ｃｏｎｊｕｇａｔｅ ｇｒａｄｉｅｎｔｓ） 反演㊁ 快速松弛反演
㊀ ㊀ 基金项目： 国家自然科学基金项目（ ４０９７４０７７㊁４１１６４００４） ； 广西自然科学基金项目 （ ２０１１ＧＸＮＳＦＡ０１８００３㊁ ２０１３ＧＸＮＳＦＡＡ０１９２７７） ； 桂林市 漓江学者 专项资助；桂林理工大学研究生创新项目（ ＢＳ２０１６０１）
㊀ 第 ４０ 卷第 ３ 期 ㊀ ２０１６ 年 ６ 月
ＧＥＯＰＨＹＳＩＣＡＬ ＆ ＧＥＯＣＨＥＭＩＣＡＬ ＥＸＰＬＯＲＡＴＩＯＮ
物㊀ 探㊀ 与㊀ 化㊀ 探
Ｖｏｌ．４０，Ｎｏ．３㊀ Ｊｕｎ．，２０１６㊀
ｄｏｉ： １０．１１７２０ ／ ｗｔｙｈｔ．２０１６．３．２２
熊彬，罗天涯，蔡红柱，等． 起伏地形大地电磁二维反演［ Ｊ］ ． 物探与化探，２０１６，４０（３） ：５８７－５９３．ｈｔｔｐ： ／ ／ ｄｏｉ．ｏｒｇ ／ １０．１１７２０ ／ ｗｔｙｈｔ．２０１６．３．２２ ５８７－５９３．ｈｔｔｐ： ／ ／ ｄｏｉ．ｏｒｇ ／ １０．１１７２０ ／ ｗｔｙｈｔ．２０１６．３．２２
［２４］
表场值沿 ｘ 和 ｚ 轴的导数可表示为其法向和切向分 量导数与地表倾斜角度 θ 的正㊁余弦函数的组合 －１ æ ������Ｖ ö æ ������ｘ ������ｚ ö æ ������Ｖ ö æ ������Ｖ ö ç ������ｘ ÷ ç ������ｎ ������ｎ ÷ ç ������ｎ ÷ ｓｉｎθ － ｃｏｓθ ö － １ ç ������ｎ ÷ æ ç ÷＝ç ÷ ç ÷＝ç ç ÷， ÷ ｓｉｎθ ø ç ������Ｖ ÷ ç ������Ｖ ÷ ç ������ｘ ������ｚ ÷ ç ������Ｖ ÷ è ｃｏｓθ ç ÷ ç ÷ ç ÷ ç ÷ è ������ｚ ø è ������ｔ ������ｔ ø è ������ｔ ø è ������ｔ ø （６） 利用含有地表节点的单元构造方程 ｓ Ｋｓ Ｖｓ ＝ Ｖｎ ， （７） 式 中， Ｖ ｓ