CN104316972A - 一种磁源重力视强度反演成像方法 - Google Patents

Abstract

本发明涉及一种磁源重力视强度反演成像方法，其包括以下步骤：将磁异常转换成空间域的磁源重力异常；选取不同的切割半径，采用插值切割法对空间域的磁源重力异常进行不同深度空间层的磁源重力异常分离，获取各空间层的剩余磁源重力异常；采用迭代法将各空间层的剩余磁源重力异常向下延拓到各空间层的磁性体的顶界面上，获取磁性体顶界面的磁源重力异常；对磁性体顶界面的磁源重力异常的视强度参数进行反演，获取各空间层的视强度参数。本发明具有反演精度高且能够降低反演多解性的优点，本发明可以广泛应用于深部矿产勘探和深部地质构造研究中。

Description

一种磁源重力视强度反演成像方法

技术领域

本发明涉及一种反演成像方法，特别是关于应用地球物理学领域中的一种磁源重力视强度反演成像方法。

背景技术

目前，国内外用于磁场物性(磁化强度)反演的方法分为两类：线性反演和非线性反演。具体的算法大致分为广义最小二乘法(Gauss法)、最速下降法、阻尼最小二乘法、奇异值分解法、人工神经网络法BP算法、遗传算法、模拟退火算法和基于Radon变换、小波多尺度边缘的算法等。上述反演方法在理论创新及小面积实际资料反演应用中都取得了一定的成绩，但在用于大面积的三维物性反演的实际生产当中没有发挥到应有的作用。其主要存在问题：1)物性参数与几何参数全局反演时，需要解大型线性代数方程组，计算时间长；2)求解方程组数学病态十分严重，解的稳定性很差；3)磁源偶极正负异常相互叠加效应导致异常反演解精度差、分辨率低。

发明内容

针对上述问题，本发明的目的是提供一种反演精度高且能够降低反演多解性的磁源重力视强度反演成像方法。

为实现上述目的，本发明采取以下技术方案：一种磁源重力视强度反演成像方法，其包括以下步骤：1)将磁异常△T(x,y)转换成空间域的磁源重力异常△g(x,y)；2)选取不同的切割半径，采用插值切割法对空间域的磁源重力异常△g(x,y)进行不同深度空间层的磁源重力异常分离，获取各空间层的剩余磁源重力异常△g_j(x,y,0)；3)采用迭代法将各空间层的剩余磁源重力异常△g_j(x,y,0)向下延拓到各空间层的磁性体的顶界面上，获取磁性体顶界面的磁源重力异常△g(x,y,h_j)；4)对磁性体顶界面的磁源重力异常△g(x,y,h_j)的视强度参数进行反演，获取各空间层的视强度参数。

所述步骤1)中，将磁异常△T(x,y)转换成空间域的磁源重力异常△g(x,y)，其包括以下步骤：(I)将地面观测到的磁异常△T(x,y)转换成频率域的磁异常即

$\mathrm{\Δ}\tilde{T}(u,v)=F\left[\mathrm{\ΔT}(x,y)\right],$

式中，u和v分别表示x方向和y方向的波数；(II)将频率域的磁异常转换成化极磁异常即

${\tilde{Z}}_{\⊥}(u,v)=\frac{u^2{+v}^2}{[i(P_{0}u+Q_{0}v)+R_0\sqrt{u^2+v^2}]}\·\frac{1}{[i(P_{1}u+Q_{1}v)+R_1\sqrt{u^2+v^2}]}\·\tilde{T}(u,v),$

式中，i表示虚分量，(P₀,Q₀,R₀)表示磁性体所在地层的方向余弦，(P₁,Q₁,R₁)表示磁性体在地表的方向余弦，(P₂,Q₂,R₂)表示磁性体在磁倾角为90°的磁化方向下的方向余弦；(III)将化极磁异常的偶极磁源作类重力的单极源处理，即对沿磁性体南北极方向积分，得到频率域的磁源重力异常即

$\mathrm{\Δ}\tilde{g}(u,v)=\frac{\mathrm{G\ρ}}{2\mathrm{\πM}}\frac{\sqrt{u^2+v^2}}{i(P_{0}u+Q_{0}v)+R_0\sqrt{u^2+v^3}}\frac{1}{i(P_{1}u+Q_{1}v)+R_1\sqrt{u^2+v^2}}\·\mathrm{\Δ}\tilde{T}(u,v),$

式中，M表示磁化强度矢量，G表示万有引力常量，ρ表示磁性体的剩余密度；(IV)对频率域的磁源重力异常进行反傅里叶变换，得到空间域的磁源重力异常△g(x,y)，即

$\mathrm{\Δg}(x,y)=F^{-1}\left[\mathrm{\Δ}\tilde{g}(u,v)\right].$

所述步骤2)中，获取各空间层的剩余磁源重力异常△g_j(x,y,0)，其具体包括以下步骤：(I)选取不同的切割半径r₁,r₂,r₃…r_n，其中，n表示切割半径的个数，采用插值切割方法对磁源重力异常△g(x,y)进行不同深度空间层的磁源重力异常分离，获取不同切割半径下的浅源异常：

L₁(x,y),L₂(x,y),L₃(x,y)…L_n(x,y)，

(II)根据不同切割半径下的浅源异常，得到各空间层的剩余磁源重力异常：

$\left\{\begin{array}{c}\mathrm{\Δ}{g}_0(x,y,0)=L_1(x,y)\\ \mathrm{\Δ}{g}_1(x,y,0)=L_2(x,y)-L_1(x,y)\\ \mathrm{\Δ}{g}_2(x,y,0)=L_3(x,y)-L_2(x,y)\\ .\\ .\\ .\\ \mathrm{\Δ}{g}_j(x,y,0)=L_{j+1}(x,y)-L_j(x,y)\\ .\\ .\\ .\\ \mathrm{\Δ}{g}_{n-1}(x,y,0)=L_n(x,y)-L_{n-1}(x,y)\\ \mathrm{\Δ}{g}_n(x,y,0)=R_n(x,y)\end{array}\right.,$

式中，△g_j(x,y,0)表示第j层磁性体在地面反映的磁源重力异常，即剩余磁源重力异常。

所述步骤3)中，采用迭代法向下延拓，获取磁性体顶界面的磁源重力异常△g(x,y,h_j)，其具体包括以下步骤：(I)将地面ΓA上的异常值△g(x,y,0)中第s点的异常值G_As放到平面ΓB的垂直投影点s上，作为平面ΓB上的异常值即并通过傅里叶变换得到ΓB上异常值的频谱(II)当地面ΓA与平面ΓB之间没有场源时，重力异常满足拉普拉斯方程，采用向上延拓公式：

$\mathrm{\Δg}(x,y,h)=F^{-1}\left[e^{-2\mathrm{\π}\left|h\right|\sqrt{u^2+v^2}}\mathrm{\Δ}\tilde{g}(u,v,0)\right],$

将平面ΓB上的异常值的频谱代入上式求出地面ΓA上s点的异常值即：

$G_{\mathrm{As}}^{\left(1\right)}(x,y,h)=F^{-1}\left[e^{2\mathrm{\π}\left|h\right|\sqrt{u^2+v^2}}{G}_{\mathrm{Bs}}^{\left(1\right)}(u,v,0)\right],$

式中，表示地面ΓA(z＝0)上的磁源重力异常值△g(x,y,0)的傅里叶变换；F^-1表示反傅里叶变换，h表示延拓高度；(III)利用G_As与的差对平面ΓB上的异常值进行修正，得到新的异常值即：

$G_{\mathrm{Bs}}^{\left(2\right)}=G_{\mathrm{Bs}}^{\left(1\right)}+l*(G_{\mathrm{As}}-G_{\mathrm{As}}^{\left(1\right)}),$

式中，l表示步长；(IV)重复步骤(II)和步骤(III)，得到如下迭代公式：

$G_{\mathrm{Bs}}^{(n+1)}=G_{\mathrm{Bs}}^{\left(n\right)}+l*(G_{\mathrm{As}}-G_{\mathrm{As}}^{\left(n\right)}),$

当 $|G_{\mathrm{As}}-G_{\mathrm{As}}^{\left(n\right)}|<\mathrm{\&epsiv;}$ 时，由式(16)可知 $|G_{\mathrm{Bs}}^{(n+1)}-G_{\mathrm{Bs}}^{\left(n\right)}|<\mathrm{k\&epsiv;},$ ε表示一个给定的趋近于零的数，则迭代次数n取20～50次；(V)采用与步骤(I)～(IV)相同的方法，计算得到平面ΓB上全部点的磁源重力异常值并利用向上延拓公式：

$\mathrm{\&Delta;g}(x,y,h)=F^{-1}\left[e^{-2\mathrm{\&pi;}\left|h\right|\sqrt{u^2+v^2}}\mathrm{\&Delta;}\tilde{g}(u,v,0)\right],$

作反傅里叶变换得到△g(x,y,0)向下延拓h_j深度的磁源重力异常△g(x,y,h_j)。

所述步骤4)中，获取各空间层的视强度参数，其具体包括以下步骤：(I)将引起磁源重力异常△g(x,y,h_j)的非规则的磁性体分割为M*N个大小规则、质心位置不同的细长垂直棱柱体；设棱柱体的质心点坐标为(x_0k,y_0k)，其剩余密度为ρ_k，棱柱体的底面边长分别为a和b，棱柱体高度为dh，则在频率域中各棱柱体在顶界面(x,y)点产生的重力异常△g_k(x,y)的频谱为：

$\mathrm{\&Delta;}{\tilde{g}}_k(u,v)=2\mathrm{\&pi;G}{\mathrm{\&rho;}}_k\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})\&CenterDot;e^{-i({\mathrm{ux}}_{0k}+{\mathrm{vy}}_{0k})},$

式中，e表示自然常数，k表示棱柱体序号，k＝1,2,3,…,M*N；(II)将各棱柱体在顶界面(x,y)点产生的重力异常的频谱进行叠加，得到整个深度层在顶界面(x,y)点产生的重力异常△g(x,y)的频谱为：

$\mathrm{\&Delta;}\tilde{g}(u,v)=2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}}){\mathrm{\&Sigma;}}_{k=1}^{\mathrm{MN}}{\mathrm{\&rho;}}_k\&CenterDot;e^{-i({\mathrm{ux}}_{0k}+{\mathrm{vy}}_k)},$

令

$\widetilde{\mathrm{\&rho;}}(u,v)={\mathrm{\&Sigma;}}_{k=1}^{\mathrm{MN}}{\mathrm{\&rho;}}_j\&CenterDot;e^{-i({\mathrm{ux}}_{0k}+{\mathrm{vy}}_k)},$

则磁源重力异常在顶界面的频谱响应为：

$\mathrm{\&Delta;}\tilde{g}(u,v)=2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})\&CenterDot;\widetilde{\mathrm{\&rho;}}(u,v),$

则频域中空间层的剩余磁源重力异常顶界面的视密度为：

$\widetilde{\mathrm{\&rho;}}(u,v)=\frac{\mathrm{\&Delta;}\tilde{g}(u,v)}{2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})};$

(III)对频域中空间层的剩余磁源重力异常顶界面的视密度进行反傅里叶变换，得到该深度层顶面各点(x,y)的剩余密度值为：

$\mathrm{\&rho;}(x,y)=F^{-1}\left[\frac{\mathrm{\&Delta;}\tilde{g}(u,v)}{2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})}\right];$

(IV)将剩余密度值ρ(x,y)代入式

$\mathrm{\&kappa;}(x,y)=\mathrm{\&rho;}(x,y)\&CenterDot;\frac{{\mathrm{\&kappa;}}^{\&prime;}}{{\mathrm{\&rho;}}^{\&prime;}},$

计算得到磁异常的视磁化率κ(x,y)；将计算得到的磁异常的视磁化率κ(x,y)代入式

M＝κ(x,y)*T₀，

计算得到磁异常的视磁化强度M；式中，κ(x,y)为磁异常的视磁化率，M为磁异常的视磁化强度，ρ(x,y)为磁源重力的视密度，T₀为正常磁场强度，κ'和ρ'分别是泊松公式转换时需要的已知磁化率值和密度值。

本发明由于采取以上技术方案，其具有以下优点：1、现有反演方法是基于空间域的优化算法，通常是物理参数及几何参数全局反演，反演过程复杂，计算精度低；与现有反演方法相比，本发明的反演方法只对物性强度参数进行反演，大大提高了反演精度，降低了反演的多解性。2、现有反演方法需求解大型线性方程组，计算机内存消耗巨大，反演时间长，因此不能进行大面积快速的三维成像反演；本发明的反演方法是在频率域进行的，运算速度快，无需求解大型线性方程组，计算机内存消耗小，因此能进行大范围、快速的三维视磁化强度反演成像。3、现有技术中偶极磁源采用传统的类重力反演方法(即用重力异常反演方法直接反演磁异常)时，会反演出两个不同的正负磁化强度参数；而本发明的磁源重力视强度反演成像方法直观，且能够反映磁性体的中心位置。4、现有反演方法大多是针对单一模型的，而对于多源模型反演，则因磁异常相互叠加，异常的等效性使得反演更加困难，反演精度差；而本发明的磁源重力视强度反演成像方法采用多源偶极磁源，能进行类似重力异常的单极叠加异常分离，并能进行基于磁源外(含磁源顶界面)任意深度空间磁异常信息获取的高分辨率视磁化强度反演求解。5、本发明的反演方法能够建立三维空间的视磁化强度属性分布。基于以上优点，本发明可以广泛应用于深部矿产勘探和深部地质构造研究中。

附图说明

图1是本发明磁源重力视强度反演成像方法的流程图；

图2是采用迭代法向下延拓的示意图；

图3是理论模型空间分布及地面投影图；其中，图(a)是理论模型空间分布图，图(b)是理论模型的地面投影图；

图4是理论模型地面磁异常示意图；

图5是地面磁异常化向地磁极异常示意图；

图6是模型地面磁源重力异常示意图；

图7是0～200米空间层剩余磁源重力异常图；

图8是200～500米空间层剩余磁源重力异常图；

图9是200米深顶界面磁源重力异常图；

图10是500米深顶界面磁源重力异常图；

图11是200米深顶界面磁源重力视密度反演图；

图12是500米深顶界面磁源重力视密度反演图；

图13是200米深顶界面磁源重力视磁化强度反演图；

图14是500米深顶界面磁源重力视磁化强度反演图。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

如图1所示，本发明磁源重力视强度反演成像方法包括以下步骤：

1)将磁异常△T(x,y)转换成空间域的磁源重力异常△g(x,y)；

地下地质体间存在密度差异或磁性差异时，地面上可观测到其相应的重力异常△g(x,y)或磁异常△T(x,y)，由于重磁场满足拉普拉斯位场理论，因此数学上可做相互转换，即重力异常可以用梯度的形式换算成重力源磁异常、磁异常也可以用积分的形式转换成磁源重力异常。下面具体说明磁异常转换成磁源重力异常的实现过程：

(I)将地面观测到的磁力异常△T(x,y)转换成频率域的磁异常

对地面观测到的磁异常△T(x,y)进行傅里叶变换，得到地面磁异常△T(x,y)的频谱即

$\mathrm{\&Delta;}\tilde{T}(u,v)=F\left[\mathrm{\&Delta;T}(x,y)\right]---\left(1\right)$

式中，u和v分别表示x方向和y方向的波数。

(II)将频率域的磁异常换算成具有重力方向的磁异常，即转换成化极磁异常

设定地表实测磁力异常为△T(x,y)的磁性体其初始磁倾角为I₁，磁偏角为D₁，该磁性体在地表的方向余弦为(P₁,Q₁,R₁)；该磁性体所在地层的地磁倾角为I₀，地磁偏角为D₀，该磁性体所在地层的方向余弦为(P₀,Q₀,R₀)；该磁性体的磁化方向从磁倾角为I₁转换到磁倾角为I₂时，磁偏角为D₂，该磁性体在该磁化方向下的方向余弦为(P₂,Q₂,R₂)，则将磁化方向从磁倾角为I₁转换成到磁倾角为I₂时，磁化方向转化公式为：

${\tilde{Z}}_a(u,v)=\frac{\sqrt{u^2{+v}^2}}{[i(P_{0}u+Q_{0}v)+R_0\sqrt{u^2+v^2}]}\&CenterDot;\frac{[i(P_{2}u+Q_{2}v)+R_2\sqrt{u^2+v^2}]}{[i(P_{1}u+Q_{1}v)+R_1\sqrt{u^2+v^2}]}\&CenterDot;\tilde{T}(u,v)---\left(2\right)$

式中，i表示虚分量，即i²＝-1。该磁性体在地表的方向余弦(P₁,Q₁,R₁)为：

$\left\{\begin{array}{c}{P}_1=\cos I_1\cos D_1\\ Q_1=\cos I_1\sin D_1\\ R_1=\sin I_1\end{array}\right.,$

该磁性体所在地层的方向余弦(P₀,Q₀,R₀)为：

$\left\{\begin{array}{c}{P}_0=\cos I_0\cos D_0\\ Q_0=\cos I_0\sin D_0\\ R_0=\sin I_0\end{array}\right.,$

该磁性体在磁倾角为I₂的磁化方向下的方向余弦(P₂,Q₂,R₂)为：

$\left\{\begin{array}{c}{P}_2=\cos I_2\cos D_2\\ Q_2=\cos I_2\sin D_2\\ R_2=\sin I_2\end{array}\right.---\left(3\right)$

当该磁性体的磁化方向从磁倾角为I₁转换到垂直磁化方向(化向地磁极)时，即I₂＝90°，将I₂＝90°代入式(3)中，则该磁性体在磁倾角为I₂的磁化方向下的方向余弦(P₂,Q₂,R₂)为：

$\left\{\begin{array}{c}{P}_2=0\\ Q_2=0\\ R_2=1\end{array}\right.---\left(4\right)$

将式(4)带入式(2)中，得到化向地磁极的频谱为：

${\tilde{Z}}_{\&perp;}(u,v)=\frac{u^2{+v}^2}{[i(P_{0}u+Q_{0}v)+R_0\sqrt{u^2+v^2}]}\&CenterDot;\frac{1}{[i(P_{1}u+Q_{1}v)+R_1\sqrt{u^2+v^2}]}\&CenterDot;\tilde{T}(u,v)---\left(5\right)$

(III)将化极磁异常的偶极磁源作类重力的单极源处理，即对沿磁性体南北极方向积分，得到频率域的磁源重力异常

由于重磁位梯度值分别对应梯度方向的磁场强度和重力场强度值，即

$\left\{\begin{array}{c}\frac{\&PartialD;U}{\&PartialD;z}={\tilde{Z}}_{\&perp;}(u,v)\\ \frac{\&PartialD;V}{\&PartialD;z}=\mathrm{\&Delta;}\tilde{g}(u,v)\end{array}\right.---\left(6\right)$

且磁位U与引力位V之间的关系满足泊松公式，即

$U=-\frac{1}{4\mathrm{\&pi;G\&rho;}}M\&CenterDot;\mathrm{grad}V---\left(7\right)$

式中，M表示磁化强度矢量，G表示万有引力常量，ρ表示磁性体的剩余密度。

由式(5)～式(7)计算得到频率域的磁源重力异常为：

$\mathrm{\&Delta;}\tilde{g}(u,v)=\frac{\mathrm{G\&rho;}}{2\mathrm{\&pi;M}}\frac{\sqrt{u^2+v^2}}{i(P_{0}u+Q_{0}v)+R_0\sqrt{u^2+v^3}}\frac{1}{i(P_{1}u+Q_{1}v)+R_1\sqrt{u^2+v^2}}\&CenterDot;\mathrm{\&Delta;}\tilde{T}(u,v)---\left(8\right)$

(IV)对频率域的磁源重力异常进行反傅里叶变换，得到空间域的磁源重力异常△g(x,y)，即

$\mathrm{\&Delta;g}(x,y)=F^{-1}\left[\mathrm{\&Delta;}\tilde{g}(u,v)\right]---\left(9\right)$

2)选取不同的切割半径，采用插值切割法对空间域的磁源重力异常△g(x,y)进行不同深度空间层的磁源重力异常分离，获取各空间层的剩余磁源重力异常△g_j(x,y,0)，其具体包括以下步骤：

(I)选取不同的切割半径r₁,r₂,r₃…r_n，其中，n表示切割半径的个数，采用插值切割方法对磁源重力异常△g(x,y)进行不同深度空间层的磁源重力异常分离，获取不同切割半径下的浅源异常：

L₁(x,y),L₂(x,y),L₃(x,y)…L_n(x,y)          (10)

(II)根据不同切割半径下的浅源异常，得到各空间层的剩余磁源重力异常：

$\left\{\begin{array}{c}\mathrm{\&Delta;}{g}_0(x,y,0)=L_1(x,y)\\ \mathrm{\&Delta;}{g}_1(x,y,0)=L_2(x,y)-L_1(x,y)\\ \mathrm{\&Delta;}{g}_2(x,y,0)=L_3(x,y)-L_2(x,y)\\ .\\ .\\ .\\ \mathrm{\&Delta;}{g}_j(x,y,0)=L_{j+1}(x,y)-L_j(x,y)\\ .\\ .\\ .\\ \mathrm{\&Delta;}{g}_{n-1}(x,y,0)=L_n(x,y)-L_{n-1}(x,y)\\ \mathrm{\&Delta;}{g}_n(x,y,0)=R_n(x,y)\end{array}\right.---\left(11\right)$

式中，△g_j(x,y,0)表示第j(j＝0,1,2,3…n)层磁性体在地面反映的磁源重力异常，即剩余磁源重力异常。

3)采用迭代法将各空间层的剩余磁源重力异常△g_j(x,y,0)向下延拓到各空间层的磁性体的顶界面上，获取磁性体顶界面的磁源重力异常△g(x,y,h_j)，其中，△g(x,y,h_j)表示将j层剩余磁源重力异常△g_j(x,y,0)向下延拓至地面以下h_j深度，获取磁性体顶界面的磁源重力异常，深度h_j为：

$\left\{\begin{array}{c}{h}_0=0,\\ h_1=r_1,\\ .\\ .\\ .,\\ h_j=r_j,\\ .\\ .\\ .,\\ h_n=r_n\end{array}\right.---\left(12\right)$

采用迭代法向下延拓，获取磁性体顶界面的磁源重力异常△g(x,y,h_j)的具体步骤包括：

(I)如图2所示，将地面ΓA上的异常值△g(x,y,0)中第s点的异常值G_As放到平面ΓB的垂直投影点s上，作为平面ΓB上的异常值即并通过傅里叶变换得到ΓB上异常值的频谱

(II)当地面ΓA与平面ΓB之间没有场源时，重力异常满足拉普拉斯方程，采用向上延拓公式：

$\mathrm{\&Delta;g}(x,y,h)=F^{-1}\left[e^{-2\mathrm{\&pi;}\left|h\right|\sqrt{u^2+v^2}}\mathrm{\&Delta;}\tilde{g}(u,v,0)\right]---\left(13\right)$

将平面ΓB上的异常值的频谱代入上式求出地面ΓA上s点的异常值即：

$G_{\mathrm{As}}^{\left(1\right)}(x,y,h)=F^{-1}\left[e^{2\mathrm{\&pi;}\left|h\right|\sqrt{u^2+v^2}}{G}_{\mathrm{Bs}}^{\left(1\right)}(u,v,0)\right]---\left(14\right)$

式(13)中，表示地面ΓA(z＝0)上的磁源重力异常值△g(x,y,0)的傅里叶变换；F^-1表示反傅里叶变换，h表示延拓高度。

(III)利用G_As与的差对平面ΓB上的异常值进行修正，得到新的异常值即：

$G_{\mathrm{Bs}}^{\left(2\right)}=G_{\mathrm{Bs}}^{\left(1\right)}+l*(G_{\mathrm{As}}-G_{\mathrm{As}}^{\left(1\right)})---\left(15\right)$

式中，l表示步长，一般取值为1。

(IV)重复步骤(II)和步骤(III)，得到如下迭代公式：

$G_{\mathrm{Bs}}^{(n+1)}=G_{\mathrm{Bs}}^{\left(n\right)}+l*(G_{\mathrm{As}}-G_{\mathrm{As}}^{\left(n\right)})---\left(16\right)$

当 $|G_{\mathrm{As}}-G_{\mathrm{As}}^{\left(n\right)}|<\mathrm{\&epsiv;}$ 时，由式(16)可知 $|G_{\mathrm{Bs}}^{(n+1)}-G_{\mathrm{Bs}}^{\left(n\right)}|<\mathrm{k\&epsiv;},$ ε表示一个给定的趋近于零的数，则 $G_{\mathrm{Bs}}^{(n+1)}\&ap;G_{\mathrm{Bs}}^{\left(n\right)}.$

一般迭代次数n取20～50次。

(V)采用与步骤(I)～(IV)相同的方法，计算得到平面ΓB上全部点的磁源重力异常值并利用向上延拓公式(13)得到△g(x,y,0)向下延拓hj深度的磁源重力异常△g(x,y,h_j)。

4)对顶界面磁源重力异常△g(x,y,h_j)的视强度参数进行反演，获取各空间层的视强度参数，其具体包括以下步骤：

(I)将引起磁源重力异常△g(x,y,h_j)的非规则的磁性体分割为M*N个大小规则、质心位置不同的细长垂直棱柱体。

设棱柱体的质心点坐标为(x_0k,y_0k)，其剩余密度为ρ_k(k＝1,2,3,…,M*N)，棱柱体的底面边长分别为a和b，棱柱体高度为dh，则在频率域中各棱柱体在顶界面(x,y)点产生的重力异常△gk(x,y)的频谱为：

$\mathrm{\&Delta;}{\tilde{g}}_k(u,v)=2\mathrm{\&pi;G}{\mathrm{\&rho;}}_k\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})\&CenterDot;e^{-i({\mathrm{ux}}_{0k}+{\mathrm{vy}}_{0k})}---\left(17\right)$

式中，e表示自然常数，k表示棱柱体序号。

(II)将各棱柱体在顶界面(x,y)点产生的重力异常的频谱进行叠加，得到整个深度层在顶界面(x,y)点产生的重力异常△g(x,y)的频谱为：

$\mathrm{\&Delta;}\tilde{g}(u,v)=2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}}){\mathrm{\&Sigma;}}_{k=1}^{\mathrm{MN}}{\mathrm{\&rho;}}_k\&CenterDot;e^{-i({\mathrm{ux}}_{0k}+{\mathrm{vy}}_k)}---\left(18\right)$

令

$\widetilde{\mathrm{\&rho;}}(u,v)={\mathrm{\&Sigma;}}_{k=1}^{\mathrm{MN}}{\mathrm{\&rho;}}_j\&CenterDot;e^{-i({\mathrm{ux}}_{0k}+{\mathrm{vy}}_k)}---\left(19\right)$

则磁源重力异常在顶界面的频谱响应为：

$\mathrm{\&Delta;}\tilde{g}(u,v)=2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})\&CenterDot;\widetilde{\mathrm{\&rho;}}(u,v)---\left(20\right)$

由式(20)得到频域中空间层剩余磁源重力异常顶界面的视密度为：

$\widetilde{\mathrm{\&rho;}}(u,v)=\frac{\mathrm{\&Delta;}\tilde{g}(u,v)}{2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})}---\left(21\right)$

(III)对频域中空间层剩余磁源重力异常顶界面的视密度进行反傅里叶变换，得到该深度层顶面各点(x,y)的剩余密度值为：

$\mathrm{\&rho;}(x,y)=F^{-1}\left[\frac{\mathrm{\&Delta;}\tilde{g}(u,v)}{2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})}\right]---\left(22\right)$

(IV)根据泊松公式，磁化率与密度存在如下线性关系：

$\mathrm{\&kappa;}(x,y)=\mathrm{\&rho;}(x,y)\&CenterDot;\frac{{\mathrm{\&kappa;}}^{\&prime;}}{{\mathrm{\&rho;}}^{\&prime;}}---\left(23\right)$

M＝κ(x,y)*T₀          (24)

式中，κ(x,y)为磁异常的视磁化率，M为磁异常的视磁化强度，ρ(x,y)为磁源重力的视密度，T₀为正常磁场强度，κ'和ρ'分别是泊松公式转换时需要的已知磁化率值和密度值。

将计算得到的剩余密度值ρ(x,y)，即磁源重力的视密度，代入式(23)中，计算得到磁异常的视磁化率κ(x,y)。将计算得到的磁异常的视磁化率κ(x,y)代入式(24)中，计算得到磁异常的视磁化强度M。

实施例1：为了详细描述本发明反演成像方法的过程及技术上的优势，创建一个如图3所示的柱体模型组，各模型柱的正演几何参数和物理参数如表1所示。

表1 模型的具体参数

采用本发明反演成像方法对该模型表示的地质体进行反演，其具体实施过程为：

(1)获取如图4所示的地面磁场值，并视其为地面观测数据；

(2)对地面观测磁异常数据做重力方向转换，获取如图5所示的化极磁异常；

(3)对化极偶极磁源做单极源转换，获取如图6所示的磁源重力异常；

(4)分别利用200米和500米的切割半径对地面磁源重力进行异常分离，获取如图7所示的0～200米空间层及如图8所示的200～500米空间层的地面剩余磁源重力异常；

(5)采用迭代法对0～200米空间层剩余磁源重力进行向下延拓，获取如图9所示的200米深磁性体顶界面磁源重力异常，采用迭代法对200～500米空间层剩余磁源重力进行向下延拓，获取如图10所示的500米深磁性体顶界面磁源重力异常；

(6)分别对200米深磁性体顶界面磁源重力异常和500米深磁性体顶界面磁源重力做视密度反演，获取如图11所示的200米深顶界面磁源重力视密度反演图和如图12所示的500米深顶界面磁源重力视密度反演图；

(7)分别对200米深视密度和500米深磁性体的视密度做视磁化强度转换，获取如图13所示的200米深顶界面磁源重力视强度反演图和如图14所示的500米深顶界面磁源重力视强度反演图。

从图4～图14可以看出，采用本发明磁源重力视强度成像方法，反演得到的如图13和图14所示的参数与如图3所示的模型参数相比，有很好的精确度。其中，深源反演磁化强度为75，与理论数据80比较接近，浅源磁化强度反演值为100，与理论值100吻合，这说明方法理论可行。在实际应用中，由于选取切割半径时存在一定的随意性，因此分离出来的层间剩余异常不能真实反映实际磁性体的空间位置，采用迭代法向下延拓后的“磁性体顶界面”可能与真实磁性体的顶界面并不重合，使得反演出来的强度参数与真实参数值存在一定的差异。然而，反演参数并不影响横向多个磁性体的磁性定性比较，因此，本发明将所反演的参数叫做视强度参数反演。

上述各实施例仅用于说明本发明，其中各部件的结构、连接方式和方法步骤等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (5)

1.一种磁源重力视强度反演成像方法，其包括以下步骤：

1)将磁异常△T(x,y)转换成空间域的磁源重力异常△g(x,y)；

2)选取不同的切割半径，采用插值切割法对空间域的磁源重力异常△g(x,y)进行不同深度空间层的磁源重力异常分离，获取各空间层的剩余磁源重力异常△g_j(x,y,0)；

3)采用迭代法将各空间层的剩余磁源重力异常△g_j(x,y,0)向下延拓到各空间层的磁性体的顶界面上，获取磁性体顶界面的磁源重力异常△g(x,y,h_j)；

4)对磁性体顶界面的磁源重力异常△g(x,y,h_j)的视强度参数进行反演，获取各空间层的视强度参数。

2.如权利要求1所述的一种磁源重力视强度反演成像方法，其特征在于：所述步骤1)中，将磁异常△T(x,y)转换成空间域的磁源重力异常△g(x,y)，其包括以下步骤：

(I)将地面观测到的磁异常△T(x,y)转换成频率域的磁异常即

$\mathrm{\&Delta;}\tilde{T}(u,v)=F\left[\mathrm{\&Delta;T}(x,y)\right],$

式中，u和v分别表示x方向和y方向的波数；

(II)将频率域的磁异常转换成化极磁异常即

${\tilde{Z}}_{\&perp;}(u,v)=\frac{u^2+v^2}{[i(P_{0}u+Q_{0}v)+R_0\sqrt{u^2+v^2}]}\&CenterDot;\frac{1}{[i(P_{1}u+Q_{1}v)+R_1\sqrt{u^2+v^2}]}\&CenterDot;\tilde{T}(u,v),$

式中，i表示虚分量，(P₀,Q₀,R₀)表示磁性体所在地层的方向余弦，(P₁,Q₁,R₁)表示磁性体在地表的方向余弦，(P₂,Q₂,R₂)表示磁性体在磁倾角为90°的磁化方向下的方向余弦；

(III)将化极磁异常的偶极磁源作类重力的单极源处理，即对沿磁性体南北极方向积分，得到频率域的磁源重力异常即

$\mathrm{\&Delta;}\tilde{g}(u,v)=\frac{\mathrm{G\&rho;}}{2\mathrm{\&pi;M}}\frac{\sqrt{u^2+v^2}}{i(P_{0}u+Q_{0}v)+R_0\sqrt{u^2+v^2}}\frac{1}{i(P_{1}u+Q_{1}v)+R_1\sqrt{u^2+v^2}}\&CenterDot;\mathrm{\&Delta;}\tilde{T}(u,v),$

式中，M表示磁化强度矢量，G表示万有引力常量，ρ表示磁性体的剩余密度；

(IV)对频率域的磁源重力异常进行反傅里叶变换，得到空间域的磁源重力异常△g(x,y)，即

$\mathrm{\&Delta;g}(x,y)=F^{-1}\left[\mathrm{\&Delta;}\tilde{g}(u,v)\right].$

3.如权利要求1或2所述的一种磁源重力视强度反演成像方法，其特征在于：所述步骤2)中，获取各空间层的剩余磁源重力异常△g_j(x,y,0)，其具体包括以下步骤：

(I)选取不同的切割半径r₁,r₂,r₃…r_n，其中，n表示切割半径的个数，采用插值切割方法对磁源重力异常△g(x,y)进行不同深度空间层的磁源重力异常分离，获取不同切割半径下的浅源异常：

L₁(x,y),L₂(x,y),L₃(x,y)…L_n(x,y)，

(II)根据不同切割半径下的浅源异常，得到各空间层的剩余磁源重力异常：

$\left\{\begin{array}{c}{\mathrm{\&Delta;g}}_0(x,y,0)=L_1(x,y)\\ {\mathrm{\&Delta;g}}_1(x,y,0)=L_2(x,y)-L_1(x,y)\\ {\mathrm{\&Delta;g}}_2(x,y,0)=L_3(x,y)-L_2(x,y)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{\&Delta;g}}_j(x,y,0)=L_{j+1}(x,y)-L_j(x,y)\\ \&CenterDot;\\ \&CenterDot;\\ \&CenterDot;\\ {\mathrm{\&Delta;g}}_{n-1}(x,y,0)=L_n(x,y)-L_{n-1}(x,y)\\ {\mathrm{\&Delta;g}}_n(x,y,0)=R_n(x,y)\end{array}\right.,$

式中，△g_j(x,y,0)表示第j层磁性体在地面反映的磁源重力异常，即剩余磁源重力异常。

4.如权利要求1或2所述的一种磁源重力视强度反演成像方法，其特征在于：所述步骤3)中，采用迭代法向下延拓，获取磁性体顶界面的磁源重力异常△g(x,y,h_j)，其具体包括以下步骤：

(I)将地面ΓA上的异常值△g(x,y,0)中第s点的异常值G_As放到平面ΓB的垂直投影点s上，作为平面ΓB上的异常值即并通过傅里叶变换得到ΓB上异常值的频谱

(II)当地面ΓA与平面ΓB之间没有场源时，重力异常满足拉普拉斯方程，采用向上延拓公式：

$\mathrm{\&Delta;g}(x,y,h)=F^{-1}\left[e^{-2\mathrm{\&pi;}\left|h\right|\sqrt{u^2+v^2}}\mathrm{\&Delta;}\tilde{g}(u,v,0)\right],$

将平面ΓB上的异常值的频谱代入上式求出地面ΓA上s点的异常值即：

$G_{\mathrm{As}}^{\left(1\right)}(x,y,h)=F^{-1}\left[e^{-2\mathrm{\&pi;}\left|h\right|\sqrt{u^2+v^2}}{G}_{\mathrm{Bs}}^{\left(1\right)}(u,v,0)\right],$

式中，表示地面ΓA(z＝0)上的磁源重力异常值△g(x,y,0)的傅里叶变换；F^-1表示反傅里叶变换，h表示延拓高度；

(III)利用G_As与的差对平面ΓB上的异常值进行修正，得到新的异常值即：

$G_{\mathrm{Bs}}^{\left(2\right)}=G_{\mathrm{Bs}}^{\left(1\right)}+l*(G_{\mathrm{As}}-G_{\mathrm{As}}^{\left(1\right)}),$

式中，l表示步长；

(IV)重复步骤(II)和步骤(III)，得到如下迭代公式：

$G_{\mathrm{Bs}}^{(n+1)}=G_{\mathrm{Bs}}^{\left(n\right)}+l*(G_{\mathrm{As}}-G_{\mathrm{As}}^{\left(n\right)}),$

当时，由式(16)可知ε表示一个给定的趋近于零的数，则迭代次数n取20～50次；

(V)采用与步骤(I)～(IV)相同的方法，计算得到平面ΓB上全部点的磁源重力异常值并利用向上延拓公式：

$\mathrm{\&Delta;g}(x,y,h)=F^{-1}\left[e^{-2\mathrm{\&pi;}\left|h\right|\sqrt{u^2+v^2}}\mathrm{\&Delta;}\tilde{g}(u,v,0)\right],$

作反傅里叶变换得到△g(x,y,0)向下延拓h_j深度的磁源重力异常△g(x,y,h_j)。

5.如权利要求1或2所述的一种磁源重力视强度反演成像方法，其特征在于：所述步骤4)中，获取各空间层的视强度参数，其具体包括以下步骤：

(I)将引起磁源重力异常△g(x,y,h_j)的非规则的磁性体分割为M*N个大小规则、质心位置不同的细长垂直棱柱体；

设棱柱体的质心点坐标为(x_0k,y_0k)，其剩余密度为ρ_k，棱柱体的底面边长分别为a和b，棱柱体高度为dh，则在频率域中各棱柱体在顶界面(x,y)点产生的重力异常△g_k(x,y)的频谱为：

$\mathrm{\&Delta;}{\tilde{g}}_k(u,v)=2\mathrm{\&pi;G}{\mathrm{\&rho;}}_k\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})\&CenterDot;e^{-i({\mathrm{ux}}_{0k}+{\mathrm{vy}}_{0k})},$

式中，e表示自然常数，k表示棱柱体序号，k＝1,2,3,…,M*N；

(II)将各棱柱体在顶界面(x,y)点产生的重力异常的频谱进行叠加，得到整个深度层在顶界面(x,y)点产生的重力异常△g(x,y)的频谱为：

$\mathrm{\&Delta;}\tilde{g}(u,v)=2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}}){\mathrm{\&Sigma;}}_{k=1}^{\mathrm{MN}}{\mathrm{\&rho;}}_k\&CenterDot;e^{-i({\mathrm{ux}}_{0k}+{\mathrm{vy}}_k)},$

令

$\widetilde{\mathrm{\&rho;}}(u,v)={\mathrm{\&Sigma;}}_{k=1}^{\mathrm{MN}}{\mathrm{\&rho;}}_j\&CenterDot;e^{-i({\mathrm{ux}}_{0k}+{\mathrm{vy}}_k)},$

则磁源重力异常在顶界面的频谱响应为：

$\mathrm{\&Delta;}\tilde{g}(u,v)=2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{by}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})\&CenterDot;\widetilde{\mathrm{\&rho;}}(u,v),$

则频域中空间层的剩余磁源重力异常顶界面的视密度为：

$\widetilde{\mathrm{\&rho;}}(u,v)=\frac{\mathrm{\&Delta;}\tilde{g}(u,v)}{2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{by}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})};$

(III)对频域中空间层的剩余磁源重力异常顶界面的视密度进行反傅里叶变换，得到该深度层顶面各点(x,y)的剩余密度值为：

$\mathrm{\&rho;}(x,y)=F^{-1}\left[\frac{\mathrm{\&Delta;}\tilde{g}(u,v)}{2\mathrm{\&pi;G}\&CenterDot;\frac{4}{\mathrm{uv}}\sin \left(\frac{\mathrm{au}}{2}\right)\sin \left(\frac{\mathrm{bv}}{2}\right)\&CenterDot;\frac{1}{r}\&CenterDot;(1-e^{-\mathrm{dhr}})}\right];$

(IV)将剩余密度值ρ(x,y)代入式

$\mathrm{\&kappa;}(x,y)=\mathrm{\&rho;}(x,y)\&CenterDot;\frac{{\mathrm{\&kappa;}}^{\&prime;}}{{\mathrm{\&rho;}}^{\&prime;}},$

计算得到磁异常的视磁化率κ(x,y)；将计算得到的磁异常的视磁化率κ(x,y)代入式

M＝κ(x,y)*T₀，

计算得到磁异常的视磁化强度M；

式中，κ(x,y)为磁异常的视磁化率，M为磁异常的视磁化强度，ρ(x,y)为磁源重力的视密度，T₀为正常磁场强度，κ'和ρ'分别是泊松公式转换时需要的已知磁化率值和密度值。