起伏地形条件下的磁异常二维聚焦反演及应用

引用本文

匡星涛, 吴健生, 杨海, 郑广如, 朱晓颖. 起伏地形条件下的磁异常二维聚焦反演及应用[J]. 物探与化探, 2016,40(4): 788-797.
KUANG Xing-Tao, WU Jian-Sheng, YANG Hai, ZHENG Guang-Ru, ZHU Xiao-Ying. Two dimensional inversion of magnetic anomaly and its application to undulating terrain[J]. Geophysical and Geochemical Exploration, 2016,40(4): 788-797.
Doi:10.11720/wtyht.2016.4.24 复制到剪切板

Permissions

《物探与化探》编辑部所有

起伏地形条件下的磁异常二维聚焦反演及应用

匡星涛¹, 吴健生², 杨海¹, 郑广如¹, 朱晓颖¹

1.中国国土资源航空物探遥感中心,北京 100083

2.同济大学海洋地质国家重点实验室,上海 200092

作者简介: 匡星涛(1988-),男,2014年获固体地球物理专业硕士学位。主要研究方向为重、磁数据处理与反演。

收稿日期: 2015-12-15

基金: 中国地质调查局“塔里木东北部地区航空重磁综合调查”项目(1212031513024)

摘要

地形对磁异常的反演和解释有很大的影响,特别是磁异常体贴近地表分布时更是如此,为了降低这种影响,笔者推导出了基于有限延深的二度厚板状体的带地形的磁异常正演公式,并提出了一种基于起伏地形的地下网格模型剖分的新方法。同时,聚焦反演对物性体有很强的聚焦效果,如果能够预先判断磁异常主要由地下浅层的强磁性体所引起,则可以考虑采用聚焦反演的方法得到磁性体的位置,显然这种方法非常适合铁矿体的反演。笔者将带地形的二维磁异常正演与聚焦反演结合起来,进行了相应的模型试验,并与水平地形下的聚焦反演及起伏地形下的光滑反演进行对比,验证了这种结合的有效性。最后,将该方法应用于新疆特克斯的磁异常剖面,对该剖面的磁性体分布做出了合理的解释。

关键词: 起伏地形; 二度体正演; 聚焦反演; 特克斯磁异常

中图分类号:P631 文献标志码:A 文章编号:1000-8918(2016)04-0788-10 doi: 10.11720/wtyht.2016.4.24

Two dimensional inversion of magnetic anomaly and its application to undulating terrain

KUANG Xing-Tao¹, WU Jian-Sheng², YANG Hai¹, ZHENG Guang-Ru¹, ZHU Xiao-Ying¹

1.China Aero Geophysical Survey and Remote Sensing Center for Land and Resources,Beijing 100083,China

2.State Key Laboratory of Marine Geology,Tongji University,Shanghai 200092,China

Abstract

Topography has a great influence on the inversion and interpretation of magnetic anomaly,especially where the abnormal body is distributed close to the ground.In order to reduce this influence,the authors deduced limited-deep two-dimensional thick plate forward formula under the rugged topography,and put forward a new method to divide the underground objects under the rugged topography.In addition,focusing inversion has good focusing effect for physical property;if we can predict that the magnetic anomaly is caused by ferromagnetic body,then focusing inversion can be applied to obtain better situation.It is obvious that this method is suitable for the inversion of ferromagnetic body.The authors combined two-dimensional magnetic anomaly forwarding and focusing inversion under the rugged topography,and carried out corresponding model test,then compared this method with focusing inversion under the flat topography with max smoothness inversion under the rugged topography,and the results show the effectiveness of the method put forward by the authors.Finally,the authors applied the method to magnetic anomaly profile obtained from Tex,Xinjiang,and made a reasonable explanation for the ferromagnetic body.

Keyword: rugged topography; 2D magnetic forwarding; focusing inversion; Tex magnetic anomaly

Show Figures

在地质勘探领域, 地形对地球物理数据的处理、反演与解释都会产生一些影响。尤其在磁法勘探领域, 观测磁异常中包含地形和地表以下各磁性体的磁异常^[1], 由于地表可能含有大量火成岩或者磁性矿体, 如果直接对数据进行地形改正(曲化平), 势必产生很大的误差, 这种误差比重力数据的地形改正大得多。另外, 曲化平的有限项调和级数法、有限元法、格林函数法、等效源法等方法计算量大, 而且计算方程常出现病态, 迭代不易收敛。所以, 将起伏地形条件下的磁异常曲化平之后再进行其他处理与解释存在一定的问题。眭素文等人采用二度半体磁异常正演公式, 研制了最优化反演方法在原起伏地形测量面上直接解释的软件, 该方法简单实用, 可视化强, 但同时具有较强的人为性^[2]。安玉林等研究了起伏地形上规则二度体复数域的重、磁场正演和直接反演方法, 该方法给出了规则二度体复数域正演公式的简单形式, 采用线性化方法, 建立线性方程组, 再结合复高阶方程解法, 比较成功地解决了起伏地形二度体重、磁场的直接反演问题^[3]。由于该方法需要分步求解多个线性方程组, 方程的病态性没有很好地解决。文中指出“ 压制小奇异值” 的“ 奇异值解法” 效果仍然不能令人满意。

每一种地球物理方法都有局限性, 针对不同的地质任务和目标应选用合适的反演和解释方法。地形起伏较大的山区, 如果观测到水平梯度较大的强磁异常, 可以设想这种磁异常是埋深较浅的强磁性体所引起的。此时显然不适合将数据曲化平再进行反演解释, 而采用最优化解释软件其人为性较大。如果能研究出一种适合于山区浅表强磁性地质情况的正、反演方法, 将其结果与欧拉反演和最优化软件反演结果进行对比, 则可以给出更可信的推断。笔者提出基于地形起伏的地下模型的网格剖分方法, 使带地形的磁异常正、反演方法和水平地面的正、反演方法一致。Last和Kubik^[4]首先提出最小化场源体积的紧致反演方法, 有效地解决了以往最小模型、最光滑模型以及最平缓模型约束所不能解决的地质体边界的问题。Silva和Barbosa^[5]也对最小体积反演方法进行了研究, 提出了解释复杂场源分布的重、磁交互反演方法。Pilkingtong^[6]提出预条件的共轭梯度法反演磁异常数据, 提高了反演效率。1998年Li Yaoguo和Oldenburg^[7]采用深度加权函数, 构建了三维重力反演框架。Rudin^[10]提出了总差分泛函稳定器, 相对于最光滑模型具有相对聚焦作用。1994年Acar和Vogel^[11]改进了Rudin的方法, 使模型参数不需要连续而只需分段光滑。1999年Portniaguine和Zhdanov^[12]总结前人结果, 提出基于Tikonov正则化方法的聚焦反演成像法, 与传统稳定泛函相比, 该方法对地质体边界识别效果较好。2002年Portniaguine和Zhdanov^[13]使用聚焦反演方法进行了磁异常三维反演模型试验, 并成功应用于加拿大北部实际航磁数据中。2004年Zhdanov和Mukherjee^[14]提出位场数据的上下延拓能有效改善聚焦反演的分辨率。

笔者基于Zhdanov聚焦反演理论, 使用最小支撑泛函来突出磁异常体边界, 并对模型进行基于总灵敏度矩阵的深度加权, 将此方法与带地形的地下网格剖分方法进行结合, 实现了带地形的聚焦反演模型试验, 并对特克斯地区的实测磁数据进行反演, 取得了较好的效果。

1 起伏地形条件下的地下模型剖分方案

为了将水平地形的聚焦反演方法应用到起伏地形情况中, 对地形进行分段近似, 用一系列首尾相连的线段代替实际地形起伏。图1显示了一个实际的地形情况。通过不等距采样, 原始的48个测量点所刻画的高程曲线可以由15个采样点近似表达。这样, 原始不规则的高程曲线就可以由分段直线段所描述。

	Figure Option View Download New Window
	图1 起伏地形不等距采样近似方法示意

地下模型剖分是根据地形起伏, 将地下空间剖分成分段的平行于地表的一系列平行四边形(图2), 设不等距采样相邻两采样点的坐标为(x_i, z_i)、(x_i+₁, z_i+₁), 如图3所示, 根据实际地质情况, 设定横向、纵向剖分间距为dx_i、dz_i, 最深点埋深为D_i, 则两采样点之间横向、纵向可剖分的个数为:

$\begin{matrix} \begin{matrix} N_{xi} = int [\frac{x_{i + 1} - x_{i}}{d x_{i}}], \\ N_{zi} = int [\frac{D_{i}}{d x_{i}}], \end{matrix} \end{matrix}$

相应剖分节点的横、纵坐标为:

$\begin{matrix} \begin{matrix} x_{i} (j, k) = x_{i} + (k - 1) d x_{i}, \\ z_{i} (j, k) = z_{i} - (j - 1) d z_{i}, \end{matrix} \end{matrix}$

其中, x(j, k)_i和z(j, k)_i分别表示第i个和第i+1个采样点地下网格剖分节点的横坐标与纵坐标。

	Figure Option View Download New Window
	图2 起伏地形地下模型剖分方案示意

	Figure Option View Download New Window
	图3 相邻两点之间地下模型剖分方法

2 地形起伏条件下的磁异常正演

通过对比分析, 选择下延有限的厚板状体为研究对象, 以基本积分方程作为出发点, 推导出适用于任意起伏地形的二度体磁异常正演公式。

假设磁荷面与观测面平行(图4), 则三维磁荷面积分公式变为二维, 即:

$\begin{matrix} \begin{matrix} H_{ax} = \frac{μ_{0}}{4 π} \int \int \frac{M_{n} (x - ξ)}{[{(x - ξ)}^{2} + {(y - η)}^{2} + {(z - ζ)}^{2}]^{\frac{3}{2}}} dξdη, \\ H_{ay} = \frac{μ_{0}}{4 π} \int \int \frac{M_{n} (y - η)}{[{(x - ξ)}^{2} + {(y - η)}^{2} + {(z - ζ)}^{2}]^{\frac{3}{2}}} dξdη, \end{matrix} \end{matrix}$

$\begin{matrix} Z_{a} = \frac{μ_{0}}{4 π} \int \int \frac{M_{n} (z - ζ)}{[{(x - ξ)}^{2} + {(y - η)}^{2} + {(z - ζ)}^{2}]^{\frac{3}{2}}} d ξdη 。 (1) \end{matrix}$

若磁荷面沿走向(y轴)延伸无限长, 则式(1)可写为:

$\begin{matrix} \begin{matrix} H_{ax} = \frac{μ_{0}}{4 π} \int_{x_{Q}}^{x_{R}} M_{n} (x - ξ) dξ \int_{-}^{\infty} \frac{dη}{[{(x - ξ)}^{2} + {(y - η)}^{2} + {(z - ζ)}^{2}]^{\frac{3}{2}}}, \\ H_{ay} = \frac{μ_{0}}{4 π} \int_{x_{Q}}^{x_{R}} M_{n} dξ \int_{-}^{\infty} \frac{(y - η) dη}{[{(x - ξ)}^{2} + {(y - η)}^{2} + {(z - ζ)}^{2}]^{\frac{3}{2}}}, \\ Z_{a} = \frac{μ_{0}}{4 π} \int_{x_{Q}}^{x_{R}} M_{n} (z - ζ) dξ \int_{-}^{\infty} \frac{dη}{[{(x - ξ)}^{2} + {(y - η)}^{2} + {(z - ζ)}^{2}]^{\frac{3}{2}}} 。 (2) \end{matrix} \end{matrix}$

对式(2)作里层的广义积分, 令z=0, ζ =h:

$\begin{matrix} \begin{matrix} H_{ax} = \frac{μ_{0}}{4 π} \int_{x_{Q}}^{x_{R}} \frac{2 M_{n} (x - ξ)}{{(x - ξ)}^{2} + h^{2}} dξ, \\ H_{ay} = 0, \\ Z_{a} = \frac{μ_{0}}{4 π} \int_{x_{Q}}^{x_{R}} \frac{- 2 M_{n} h}{{(x - ξ)}^{2} + h^{2}} dξ 。 (3) \end{matrix} \end{matrix}$

模型是由一系列平行四边形所组成, 顶、底边与水平坐标轴的夹角由相邻的两个采样点之间的斜率决定, 而侧边与水平坐标轴的夹角恒为90° 。图5是一个截面为平行四边形的水平二度磁性体在起伏地表产生磁异常的示意图, z轴向下为正。该二度体在P点产生的磁异常是6个磁荷面分别在P点产生的磁异常之和。由于二度体走向无限延伸, 则这两个面在P点产生的磁异常视为0。因此, 此二度体的磁场可表示为:

$\begin{matrix} \begin{matrix} Z_{a} = Z_{aAB} + Z_{aCD} + Z_{aAC} + Z_{aBD}, \\ H_{ax} = H_{axAB} + H_{axCD} + H_{axAC} + H_{axBD} 。 (4) \end{matrix} \end{matrix}$

由于4个面与x轴不平行, 为使用式(3)进行面积分, 必须进行坐标转换。由于地形起伏, 测点与平行四边形的位置关系比较复杂, 不同的位置关系对应着不同的积分结果。在求AC与BD面产生的磁异常时, 将测点与平行四边形的位置关系划分为5种情况, 图5中的P₁、P₂、P₃、P₄、P₅分别表示测点位置位于AC左侧、AC之上、AC与BD之间、BD之上、BD右侧; 同理, 在求AB与CD面产生的磁异常时, 将测点与平行四边形的位置关系也划分为5种情况, P₆、P₇、P₈、P₉、P₁₀分别表示测点位置位于AB

上侧、AB之上、AB与CD之间、BD之上、BD下侧。

下面推导截面为平行四边形二度体在不同测点的磁异常正演公式。

	Figure Option View Download New Window
	图4 二度水平磁荷面坐标示意

	Figure Option View Download New Window
	图5 测点与截面平行四边形的位置关系示意

设A、B、C、D的坐标分别为(x_A, z_A)、(x_B, z_B)、(x_C, z_C)、(x_D, z_D), 则AB面与CD面的倾角

$\begin{matrix} δ_{1} = \arctan \frac{z_{B} - z_{A}}{x_{B} - x_{A}} 。 \end{matrix}$

由于AC面、BD面与x轴垂直, 此两面的倾角δ ₂=-90° 。

设有效磁化强度为M_S, 磁化倾角为i_S, 则AB面、CD面与有效磁化强度矢量方向的夹角γ =δ ₁-i_S。AC面、BD面与有效磁化强度矢量方向的夹角γ =δ ₂-i_S。则4个面的磁荷面密度分别为:

$\begin{matrix} \begin{matrix} σ_{AB} = M_{S} sinγ, σ_{CD} = - M_{S} sinγ, \\ σ_{AC} = - M_{S} \cos i_{S}, σ_{BD} = M_{S} \cos i_{S} 。 \end{matrix} \end{matrix}$

设测点坐标为(x, z), 二度体外的任意一点到A、B、C、D的距离为:

$\begin{matrix} \begin{matrix} r_{A} = \sqrt[]{(x - x_{A})^{2} + (z - z_{A})^{2}}, \\ r_{B} = \sqrt[]{(x - x_{B})^{2} + (z - z_{B})^{2}}, \\ r_{C} = \sqrt[]{(x - x_{C})^{2} + (z - z_{C})^{2}}, \\ r_{D} = \sqrt[]{(x - x_{D})^{2} + (z - z_{D})^{2}} 。 \end{matrix} \end{matrix}$

测点到直线AB、CD、AC、BD的距离分别为:

$\begin{matrix} \begin{matrix} h_{AB} = \frac{|(z_{B} - z_{A}) x + (x_{A} - x_{B}) z + x_{B} z_{A} - x_{A} z_{B} |}{\sqrt[]{(z_{B} - z_{A})^{2} + (x_{A} - x_{B})^{2}}}, \\ h_{CD} = \frac{|(z_{D} - z_{C}) x + (x_{C} - x_{D}) z + x_{D} z_{C} - x_{C} z_{D} |}{\sqrt[]{(z_{D} - z_{C})^{2} + (x_{C} - x_{D})^{2}}}, \\ h_{AC} = | x - x_{A} |, \\ h_{BD} = | x - x_{B} |。 \end{matrix} \end{matrix}$

从上式可以看出, 测点到4条直线的距离可能为0, 即式(3)中h=0, 则式(3)可简化为:

$\begin{matrix} \begin{matrix} H_{a x} = \frac{μ_{0}}{4 π} \int_{x_{Q}}^{x_{R}} \frac{2 M_{n}}{(x - ξ)} dξ, \\ H_{ay} = 0, Z_{a} = 0 。 (4) \end{matrix} \end{matrix}$

当x=x_A或x=x_B时, 即测点在AC或BD上时, h_AC=0, h_BD=0, 那么将h=0代入式(3)并积分得:

$\begin{matrix} \begin{matrix} H_{ax} = \frac{- μ_{0} σ_{QR}}{4 π} \ln \frac{(x - x_{R})^{2} + h^{2}}{(x - x_{Q})^{2} + h^{2}}, \\ Z_{a} = \frac{- μ_{0} σ_{QR}}{2 π} (\arctan \frac{x - x_{Q}}{h} - \arctan \frac{x - x_{R}}{h}) 。 (5) \end{matrix} \end{matrix}$

对式(4)积分可得:

$\begin{matrix} \begin{matrix} H_{ax} = \frac{- μ_{0} σ_{QR}}{2 π} \ln |\frac{x - x_{R}}{x - x_{Q}}|, \\ H_{ay} = 0, Z_{a} = 0 。 (6) \end{matrix} \end{matrix}$

由式(5)和式(6)可计算二度体外任意一点的磁场值, 计算AC与BD面产生的磁场时, 将坐标系逆时针旋转δ ₂度; 计算AB面与CD面产生的磁场时, 将坐标系逆时针旋转δ ₁度。计算完磁场分量之后, 再进行坐标反变换, 即可得到测点处的磁场值。

用二度体正演公式计算截面为平行四边形二度体在起伏地表的正演值, 并与眭素文等^[2]研制的二度半体磁异常拟合软件的正演结果进行对比(图6), 可知该正演结果较为可靠, 图6b为起伏地形时模型的横截面, 磁化强度2 A/m, 磁化倾角60° , 磁化强度水平分量矢量与二度体走向夹角为10° 。

3 基于最小支撑函数的聚焦反演原理

一般将地球物理反演问题转化为求解目标函数的极值问题, 若仅考虑数据拟合, 目标函数为:

$\begin{matrix} P (m) = Φ (m) = ‖ A (m) - d ‖^{2} \to \min, (7) \end{matrix}$

其中, m为模型参数, A为正演算子, d为观测数据。当观测值个数少于模型时, 为欠定问题, 无穷多解。即便问题超定, 当用高斯— 牛顿类方法求解方程组时, 由于方程组的病态性、解不稳定, 反演出来的模型m也可能与实际地质情况不符。因此, 常用模型最小二范数作为模型约束, 加强解的稳定性, 即:

$\begin{matrix} \begin{matrix} P (m) & = Φ (m) + αS (m) \\ = ‖ A (m) - d ‖^{2} + α ‖ m ‖^{2} \to \min 。 (8) \end{matrix} \end{matrix}$

	Figure Option View Download New Window
	图6 正演方法准确性验证

任何对模型的约束, 都是基于一定的先验信息或者先验假设。最小模型约束是假设整个模型空间满足高斯分布, 是整个模型空间概率最大的模型。除此之外, 常用的模型约束还有最光滑模型和最平缓模型, 但模型约束在改善方程组的条件数、提高解的稳定性的同时, 也降低了解的分辨率, 让解空间对数据的影响变的不灵敏, 最后得到的解很分散。在中、浅磁性层固体矿产勘查中, 由于磁性矿体的密度和磁化强度与围岩差别很大, 可先假设磁性矿体与围岩有较为清晰的边界, 即可使用聚焦反演方法, 使反演出来的磁性矿体体积有清晰边界, 且密度和磁化强度在合理的范围。

聚焦反演的稳定器种类较多, 常用的有三种:

$\begin{matrix} \begin{matrix} S_{βTV} (m) = \int \sqrt[]{|\nabla m |^{2} + β^{2}} dv, (9) \\ S_{βMS} (m) = \int \frac{m^{2}}{m^{2} + β^{2}} dv, (10) \\ S_{β M GS} (m) = \int \frac{\nabla m \cdot \nabla m}{\nabla m \cdot \nabla m + β^{2}} dv, (11) \end{matrix} \end{matrix}$

其中, S_{β TV}(m)为总变差稳定器, S_{β MS}(m)为最小支撑稳定器, S_{β MGS}(m)为最小梯度支撑稳定器。∇m表示模型梯度, β 为一常数, 式(9)中是改善稳定器在零值处的可微性, 式(10)和式(11)中是防止模型或模型的梯度为0时分母没有意义, 以及聚焦的作用。笔者选择最小支撑稳定器作为模型约束, 是为了使解空间中模型不为0的体积最小。设spt (m)为模型不为0的体积, 即模型参数支撑, 则式(10)为:

$\begin{matrix} \begin{matrix} S_{βMS} (m) & = \int (1 - \frac{β^{2}}{m^{2} + β^{2}}) dv \\ = spt (m) - β^{2} \int \frac{1}{m^{2} + β^{2}} dv, (12) \end{matrix} \end{matrix}$

当β 很小时, S_{β MS}(m)→ spt (m), 即最小支撑稳定器与模型参数支撑成正比, 通过最小支撑稳定器可以让模型非零空间最小, 聚焦反演目标函数可写为:

$\begin{matrix} \begin{matrix} P (m) & = Φ (m) + α S_{βMS} (m) \\ = ‖ A (m) - d ‖^{2} + α \int \frac{m^{2}}{m^{2} + β^{2}} dv \to \min 。 \\ (13) \end{matrix} \end{matrix}$

为适应网格模型空间, 需将模型离散化, 笔者采用共轭梯度法求解目标函数, 将式(10)改写成两个相同函数内积的形式:

$\begin{matrix} S_{βMS} (m) = \int (\frac{m}{\sqrt[]{m^{2} + β^{2}}}, \frac{m}{\sqrt[]{m^{2} + β^{2}}}) dv, (14) \end{matrix}$

记对角阵W_e=diag $\begin{matrix} [\frac{m}{\sqrt[]{m^{2} + β^{2}}}] \end{matrix}$ , 则式(13)可改写为:

$\begin{matrix} \begin{matrix} P (m) & = Φ (m) + α S_{βMS} (m) \\ = ‖ A (m) - d ‖^{2} + α ‖ W_{e} m ‖^{2} \to \min 。 \\ (15) \end{matrix} \end{matrix}$

若直接求解式(15), 虽然有聚焦效果, 但是由于趋肤效应, 反演出的大物性值集中于地表。但是多数情况下, 这样的反演结果与实际地质情况不相符。为了克服趋肤效应, 引入模型加权矩阵W_m:

其中, F_m为正演算子的Fre'chet偏导数矩阵, S为数据对模型参数的综合灵敏度。综合灵敏度为某数据对所有模型参数m的偏导数的平方和, 。

有了模型加权矩阵, 就可以得到加权模型参数m^ω=W_mm, 而加权偏导数矩阵 $\begin{matrix} F_{m}^{ω} \end{matrix}$ =F_m $\begin{matrix} W_{m}^{- 1} \end{matrix}$ , 数据d对 $\begin{matrix} m_{k}^{ω} \end{matrix}$ 的综合灵敏度:

通过引入模型加权矩阵, 目标函数式(15)变为:

$\begin{matrix} \begin{matrix} P (m) & = Φ (m) + α S_{βMS} (m) \\ = ‖ A (m) - d ‖^{2} + α ‖ W_{e} W_{m} m ‖^{2} \to \min 。 \\ (17) \end{matrix} \end{matrix}$

为了加快共轭梯度法收敛速率, 使用再加权模型的方法, 即:

$\begin{matrix} F_{m}^{ω} = F_{m} W_{m}^{- 1} W_{e}^{- 1}, m^{ω} = W_{m} W_{e} m, \end{matrix}$

由于磁场反演物性是线性问题, 正演算子A等价于偏导数矩阵F_m, A^ω=A $\begin{matrix} W_{m}^{- 1} \end{matrix} \begin{matrix} W_{e}^{- 1} \end{matrix}$ , 目标函数又可以改写为:

$\begin{matrix} \begin{matrix} P (m^{ω}) & = Φ (m^{ω}) + 、 α S_{βMS} (m^{ω}) \\ = ‖ A^{ω} (m^{ω}) - d ‖^{2} + α ‖ m^{ω} ‖^{2} \to \min 。 \\ (18) \end{matrix} \end{matrix}$

用共轭梯度法求解此目标函数, 第0步迭代使用梯度法, 第1步迭代时正则化因子取为:

$\begin{matrix} α_{1} = \frac{‖ A^{ω} ({m^{ω}}_{1}) - d ‖^{2}}{‖ {m^{ω}}_{1} ‖^{2}}, \end{matrix}$

随后迭代过程正则化因子按照一定的比例因子衰减, 使前期迭代主要侧重位置聚焦, 后期迭代主要侧重数据拟合。衰减方式为:

$\begin{matrix} α_{k} = α_{1} q^{k - 1}; k = 1, 2, 3, \dots; q > 0 。 \end{matrix}$

迭代结束后, 可反求出未加权的模型参数, 即实际的模型参数m= $\begin{matrix} W_{m}^{- 1} \end{matrix} \begin{matrix} W_{e}^{- 1} \end{matrix}$ m^ω。另外, 在迭代过程中, 由于数值算法原因, 解空间可能出现负值或者超出合理范围的正值, 为此, 需要构造惩罚函数, 对负值和超出合理范围的正值进行截断, 设m(r)为物性参数场, 它描述了地下空间各网格对应的模型参数值, m_ub(r)为允许出现的最大物性值。则截断算法可写为:

$\begin{matrix} \{\begin{matrix} m (r) = m_{ub} (r), 若 m (r) \geq m_{ub} (r), \\ m (r) = 0, 若 m (r) \leq m_{ub} (r) 。 \end{matrix} \end{matrix}$

4 带地形的聚焦反演模型试验

为了说明地形起伏条件下的聚焦反演效果, 将水平地形条件下的聚焦反演和凸地形条件下的光滑反演结果与聚焦反演结果进行对比(图7)。剖面长1 000 m, 深度600 m, 地下空间剖分成20 m× 20 m的网格, 理论模型尺寸100 m× 100 m, 走向延伸无限, 埋深为180 m(凸地形)和100 m(水平地形), 磁化强度1 A/m, 三种情况下反演时的惩罚函数相同。当地形为凸时, 网格是平行四边形(图7b、f); 当地形水平时, 网格是矩形(图7d)。图7a、c、e表示凸地表最大光滑解、水平地表下聚焦解以及凸地表下聚焦解对应的数据拟合情况, 三种情况下数据拟合情况较好。图7b、d、f为对应模型反演磁化强度的分布情况, 凸地形的最大光滑解主要趋近于地表, 而水平地形和凸地形的聚焦解与理论模型吻合较好。

	Figure Option View Download New Window
	图7 凸地形与水平地形单块体反演对比

	Figure Option View Download New Window
	图8 凸地形与水平地形垂向叠加双块体反演对比

图8为双块体水平地形条件下聚焦反演、凸地形条件下光滑反演与凸地形条件下聚焦反演的结果对比, 浅部模型尺寸100 m× 100 m, 埋深80 m, 深部模型尺寸300 m× 150 m, 埋深330 m, 磁化强度1 A/m, 其余与图7的参数相同。图8中可以看出, 最大光滑解很分散, 而且趋近于地表, 完全不能反映实际的块体形状。水平地表下的聚焦解对浅部模型的反映较好, 但对深部模型的反映较差且没有有效地将两个块体分开, 水平地表下聚焦反演垂向分辨力较差。凸地形条件下的聚焦结果能够将两个块体有效分开, 而且块体位置也得到了正确反映, 凸地形条件下聚焦反演方法具有较强的垂向分辨力。

图9为双块体水平地形条件下聚焦反演、凹地形条件下光滑反演与凹地形条件下聚焦反演的结果对比。聚焦反演参数与图8基本相同, 只是地形由凸改为凹、垂向分布的2个模型改为了水平分布, 埋深分别为80 m和90 m。图9b的反演结果也趋于地表, 水平地表和凹地表下的聚焦解能够较好地反演出模型位置、并将两个块体分开, 但物性集中度有所欠缺, 这表明聚焦解具有较好的水平分辨力。

	Figure Option View Download New Window
	图9 凹地形与水平地形水平叠加双块体反演对比

	Figure Option View Download New Window
	图10 凹地形与水平地形垂向叠加双块体反演对比

图10为垂向叠加双块体的水平地形条件下聚焦反演、凹地形条件下光滑反演与凹地形条件下聚焦反演的结果对比。浅部模型尺寸100 m× 100 m, 埋深60 m, 深部模型尺寸300 m× 150 m, 埋深310 m, 磁化强度1 A/m, 其余参数与图9中相同。

从图中可以看出, 最大光滑解趋近于地表, 且十分集中于凹地形的最低点, 完全不能反映实际的块体形状。水平地表下的聚焦解对浅部模型的反映较好, 但对深部模型的反映较差, 且没有有效地将两个块体分开, 这表明水平地表下聚焦反演的垂向分辨力较差。凹地形条件下的聚焦解与水平地表下的聚焦解基本相同, 也未能有效地将两个块体分开, 但反演的强磁化体向下延伸更大且已接近深部模型的上部边界。试验结果表明, 凹地形条件下聚焦反演方法优于水平地表下聚焦方法、但垂向分辨力也有限。

	Figure Option View Download New Window
	图11 凸地形与水平地形向下延伸板状体反演对比

图11为向下延伸长板状体水平地形条件下聚焦反演、凹地形条件下光滑反演结果与凹地形条件下聚焦反演的结果对比。聚焦反演参数与图7中基本相同, 理论模型100 m× 400 m, 埋深100 m(凸地形)和80 m(水平地形), 磁化强度1 A/m。从反演结果可以看出, 凸地形的最大光滑解主要趋近于地表, 而水平地形和凸地形的聚焦解与模型的顶部有一定吻合, 但凸地形的聚焦解与模型的吻合度更高。

5 实际资料处理

为了验证方法的有效性, 对新疆特克斯地区地磁异常进行聚焦反演。

新疆特克斯地区处于中高山地区, 海拔高度为1 900~2 400 m, 属于中深切割、坡陡谷深的地区。出露的地层及其岩性主要为二叠系铁木里克组(P₂t)砾岩夹砂岩、泥质粉砂岩; 下石炭统大哈拉军山组(C₁d)安山玢岩、流纹岩、凝灰岩、凝灰质砂岩、灰岩、砾岩; 侏罗系水西沟群(J_1-2sh)砾岩、砂岩、粘土岩。区内构造较发育, 主要构造为贯穿全区的东西向深大断裂。下石炭统大哈拉军山组(C₁d)岩石的磁化强度为0.1~0.6 A/m, 其他地层的磁性较弱。

图12a为新疆特克斯地区地磁异常平面剖面, 磁异常主要分布在大哈拉军山组和铁木里克组接触带以北的大哈拉军山组出露区, 东西向展布, 选取其中8号测线一段的磁异常和实测地形(图12c、d)进行二维聚焦反演。

	Figure Option View Download New Window
	图12 特克斯磁场平面剖面与其中8号线的地形与磁场剖面

	Figure Option View Download New Window
	图13 测线8的反演结果

聚焦反演的磁化倾角约 63.5° , 磁化偏角4.04° 。30次共轭梯度法迭代之后, 数据拟合差小于2%(图13)。图13b表明, 反演结果显示了三片强磁化模型单元分布区, 分别对应于图13b中的三个红色椭圆圈定的范围; 反演磁化强度的最大值为28 A/m, 与地表岩石物性测量0.1~0.6 A/m有较大差距, 说明该区地下可能存在磁性矿体, 其中南部埋藏较浅为50~100 m, 北部埋藏较深为150~300 m。

6 结论与展望

1) 基于地表起伏的平行四边形二度体网格建模方法能够适应复杂地形的建模需要。

2) 基于二度厚板状体的面积分公式, 推导出的截面为平行四边形的水平二度体的磁异常计算公式对解决起伏地形区磁异常的正、反演问题具有理论与实际意义。

3) 将聚焦反演方法与地形起伏条件下网格建模结合起来, 进行了相应的反演试验, 并与起伏地形条件下最大光滑解以及水平地形条件下的聚焦解进行了对比, 结果表明:聚焦反演结果优于最大光滑反演结果; 凸地形条件下的聚焦反演具有较强的垂向分辨能力。同时也发现, 对于延深较大的地质体或非凸地形条件下的垂向叠加地质体, 聚焦反演结果不能很好地揭示地质体深部或深部地质体的完整形态, 表明此方法还需要进一步研究。

4) 通过新疆特克斯实测磁异常的反演, 认为聚焦反演法能快速获得地下浅部强磁性体的分布情况, 为精细解释提供初始模型。

The authors have declared that no competing interests exist.

参考文献

View Option

[1]	刘沈衡. 磁性起伏地形磁异常改正方法探讨[J]. 地质学刊, 2011, 35(3): 302-306. [本文引用:1]
[2]	眭素文, 于长春, 姚长利. 起伏地形剖面重磁异常半智能处理解释软件及应用[J]. 物探与化探, 2004, 28(1): 65-68. [本文引用:2]
[3]	安玉林. 起伏地形上规则二度体复重磁场正演和直接反演[J]. 物探与化探, 2003, 27(4): 284-291. [本文引用:1]
[4]	Last B J, Kubik K. Compact gravity inversion[J]. Geophysics, 1983, 48(6): 713-721. [本文引用:1]
[5]	Barbosa V C F, Silva J B C. Generalized Compact gravity inversion[J]. Geophysics, 1994, 59(1): 57-68. [本文引用:1]
[6]	Pilkingtong M. 3-D magnetic imaging using conjugate gradient[J]. Geophysics, 1997, 62(4): 1132-1142. [本文引用:1]
[7]	Li Y G, Oldenburg D W. 3-D inversion of gravity data[J]. Geophysics, 1998, 63(1): 109-119. [本文引用:1]
[8]	Li Y G, Oldenburg D W. 3-D inversion of magnetic data[J]. Geophysics, 1996, 61(2): 394-408. [本文引用:1]
[9]	Li Y G. 3-D Inversion of gravity gradiometerdata[J]. Seg Expand ed Abstracts, 1949(1): 1470. [本文引用:1]
[10]	Rudin L I, Osher S, Fatemi E. Nonlinear total variationbased noise removal algorithm[J]. Physica D, 1992, 60: 259-268. [本文引用:1]
[11]	Acar R, Vogel C R. Analysis of total variation penaltymethods[J]. Inverse problems, 1994, 10: 1217-1229. [本文引用:1]
[12]	Portniaguine O, Zhdanov M S. Focusing geophysical inversion images[J]. Geophysics, 1999, 64(3): 874-887. [本文引用:1]
[13]	Oleg Portniaguine, Michael S. Zhdanov. 3-D magnetic inversion with data compressionand image focusing[J]. Geophysics, 2002, 67(5): 1532-1541. [本文引用:1]
[14]	Zhdanov M S, Ellis R, Mukherjee S. Three-dimensional regularized focusing inversion of gravitygradient tensor component data. Geophysics, 2004, 69(4): 925-937. [本文引用:1]
[15]	吴健生, 刘苗. 基于小波的位场数据融合[J]. 同济大学学报: 自然科学版, 2008, 36(8). [本文引用:1]
[16]	杨文采, 施志群, 侯尊泽, 等. 离散小波变换与重力异常多重分解[J]. 地球物理学报, 2001, 44(4): 534-541. [本文引用:1]
[17]	管志宁. 地磁场与磁力勘探[M]. 北京: 地质出版社, 2005. [本文引用:1]
[18]	郑广如, 张玄杰, 范子梁. 新疆西天山特克斯-霍拉山地区1∶5万航磁勘查成果报告[R]. 中国国土资源航空物探遥感中心, 2012. [本文引用:1]

2011

0.0

... 尤其在磁法勘探领域,观测磁异常中包含地形和地表以下各磁性体的磁异常^[1],由于地表可能含有大量火成岩或者磁性矿体,如果直接对数据进行地形改正(曲化平),势必产生很大的误差,这种误差比重力数据的地形改正大得多 ...

2004

0.0

... 眭素文等人采用二度半体磁异常正演公式,研制了最优化反演方法在原起伏地形测量面上直接解释的软件,该方法简单实用,可视化强,但同时具有较强的人为性^[2] ...

... 用二度体正演公式计算截面为平行四边形二度体在起伏地表的正演值,并与眭素文等^[2]研制的二度半体磁异常拟合软件的正演结果进行对比(图6),可知该正演结果较为可靠,图6b为起伏地形时模型的横截面,磁化强度2 A/m,磁化倾角60#cod#x000b0 ...

2003

0.0

... 安玉林等研究了起伏地形上规则二度体复数域的重、磁场正演和直接反演方法,该方法给出了规则二度体复数域正演公式的简单形式,采用线性化方法,建立线性方程组,再结合复高阶方程解法,比较成功地解决了起伏地形二度体重、磁场的直接反演问题^[3] ...

1983

0.0

... Last和Kubik^[4]首先提出最小化场源体积的紧致反演方法,有效地解决了以往最小模型、最光滑模型以及最平缓模型约束所不能解决的地质体边界的问题 ...

1994

0.0

... Silva和Barbosa^[5]也对最小体积反演方法进行了研究,提出了解释复杂场源分布的重、磁交互反演方法 ...

1997

0.0

... Pilkingtong^[6]提出预条件的共轭梯度法反演磁异常数据,提高了反演效率 ...

1998

0.0

... 1998年Li Yaoguo和Oldenburg^[7]采用深度加权函数,构建了三维重力反演框架 ...

1996

0.0

1992

0.0

... Rudin^[10]提出了总差分泛函稳定器,相对于最光滑模型具有相对聚焦作用 ...

1994

0.0

... 1994年Acar和Vogel^[11]改进了Rudin的方法,使模型参数不需要连续而只需分段光滑 ...

1999

0.0

... 1999年Portniaguine和Zhdanov^[12]总结前人结果,提出基于Tikonov正则化方法的聚焦反演成像法,与传统稳定泛函相比,该方法对地质体边界识别效果较好 ...

2002

0.0

... 2002年Portniaguine和Zhdanov^[13]使用聚焦反演方法进行了磁异常三维反演模型试验,并成功应用于加拿大北部实际航磁数据中 ...

2004

0.0

... 2004年Zhdanov和Mukherjee^[14]提出位场数据的上下延拓能有效改善聚焦反演的分辨率 ...

2008

0.0

2001

0.0

2005

0.0

2012

0.0