一种基于磁偶极子等效源法的地磁异常数据网格化方法

技术领域

本发明属于地磁异常数据处理领域，具体涉及一种基于磁偶极子等效源法的地磁异常数 据网格化方法。

背景技术

地磁异常数据由地壳浅层及表面的岩石、矿产及其他磁性特征产生，特征丰富、稳定性 强，在矿产探测、地质研究、运载器导航等众多方面有广泛的应用，由于直接磁场测量结果 往往分布不均匀，为了便于进行数据分析和工程应用，首先需要对测量数据进行网格化处理。

当前已有很多方法用于空间散乱数据的网格化，可以分为三类，一是采用插值方法对网 格点位置进行估值，针对某个特定的网格点，将其周围附近的已知点当作参考点，对该位置 进行估计，代表方法有克里金法、径向基函数法、反距离加权法、样条函数法、谢别德法等； 二是模型法，即利用所有已知的数据对研究区域建立统一的数学模型，并利用所建立的模型 对网格点位置的磁场值进行解算，代表方法有泰勒多项式法、矩谐模型方法、球谐或冠谐模 型方法，但这类方法的缺点是采用一个统一的表达式很难对整片区域的磁场值进行准确的描 述，尤其是对于随位置变化比较剧烈的地磁异常场，当研究区域的磁异常变化丰富并且数据 量较大时，该方法的有效性将大大减弱；三是采用智能预测算法，即将已知数据作为训练集， 对智能算法模型进行训练，然后利用训练好的模型对网格点位置的磁场值进行预测，这类方 法的优点是对边界区域的数据处理效果比较稳定，但是整体的预测精度比较有限，并且严重 依赖于训练数据的数量和质量。

此外，以上所有方法有一个共同特点，它们都是完全基于测量数据本身的数学运算过程， 而忽略了从磁场的产生原理上进行分析。地磁异常场是由地壳浅层及表面的磁性物质产生， 且任何磁性物体都可以看做是由磁偶极子组合而成，因此本文提出一种基于磁偶极子等效源 法的地磁异常场网格化方法。首先在研究区域均匀布置一定数量的磁偶极子作为磁源，根据 实测地磁异常数据，利用本专利发明的磁偶极子等效源法反演这些磁偶极子的磁矩，然后再 进行正演，计算目标网格点位置的磁场。本方法直接从地磁异常的产生机理入手，遵循磁场 分布的基本规律，网格化处理的精度和稳定性相比传统方法具有一定的优势；通过设定磁偶 极子在目标平面上的有效区域，避免了传统等效源方法中对大量数据解算时的复杂度过高的 问题，也避免了分块计算方法带来的反演精度降低现象；是当测量数据比较稀疏，或者某些 特定区域数据量较小时，传统方法将遇到巨大的瓶颈，由于无法获取足够的参考信息，在这 些位置将产生很大的误差甚至发散，但本方法依然能够通过此场正演准确计算这些位置的磁 场值；同时本文方法可以有效地抑制处理边界区域数据的边界效应，从而有效地对地磁异常 观测数据进行表达，并进行网格化。

发明内容

本发明的目的在于提供一种提高计算精度以及运算效率的基于磁偶极子等效源法的地磁 异常数据网格化方法。

本发明的目的是这样实现的：

（1）确定磁偶极子产生磁场的有效区域；

磁场观测平面距水平面高度为h₁，磁源平面在水平面下方深度为h₀，磁源平面与磁场观 测平面的距离为h=h₁+h₀，磁源平面上磁偶极子的磁矩大小为p(ξ_i,η_i,ζ_i)，磁矩方向和当地 的磁化方向一致，磁偶极子在上方任意一点(x,y,z)处产生的ΔT磁异常为，

ΔT(x,y,z)=p(ξ_i,η_i,ζ_i)·H(x,y,z,ξ_i,η_i,ζ_i,α,β,γ)

其中，

H(x,y,z,ξ_i,η_i,ζ_i,α,β,γ)=-U_xi(x,y,z)cosα-U_yi(x,y,z)cosβ-U_zi(x,y,z)cosγ

$U_{\mathrm{xi}}(x,y,z)=-\frac{1}{4\pi R_i^2}{A}_i(x,y,z),U_{\mathrm{yi}}(x,u,z)=-\frac{1}{4\pi R_i^2}{B}_i(x,y,z),$

$U_{\mathrm{zi}}(x,y,z)=-\frac{1}{4\pi R_i^2}{C}_i(x,y,z),$

A_i(x,y,z)=??2(x-ξ_i)²-(y-η_i)²-(z-ζ_i)²]cosα+3(x-ξ_i)(y-η_i)cosβ+3(x-ξ_i)(z-ζ_i)cosγ ，

B_i(x,y,z)=3(x-ξ_i)(y-η_i)cosα+[2(y-η_i)²-(x-ξ_i)²-(z-ζ_i)²]cosβ+3(y-η_i)(z-ζ_i)cosγ ，

C_i(x,y,z)=3(x-ξ_i)(z-ζ_i)cosα+3(y-η_i)(Z-ζ_i)cosβ+[2(z-ζ_i)²-(x-ξ_i)²-(y-η_i)²]cosγ

$,R_i=\sqrt{{(x-{\xi }_i)}^2+{(y-{\eta }_i)}^2+{(z-{\zeta }_i)}^2}$

其中，α，β，γ分别为磁化方向与当地地理坐标系X,Y,Z轴的夹角；

磁源平面上一磁偶极子为p，计算磁偶极子p在磁场观测平面上产生的磁场，确定该磁 偶极子的有效区域，即磁偶极子p在有效区域之外任意一点产生的磁场是磁偶极子在观测面 上产生磁场最大值的1/20；

（2）在磁源平面上布置磁偶极子；

在磁场观测区域水平面下方的平面上，以网格形式均匀布置磁偶极子组，磁偶极子的间 隔为观测数据平均间隔的5倍；磁偶极子组的分布范围大于观测数据区域，并使得最外面的 磁偶极子的有效区域内至少包含一个磁场观测值最终布置的磁偶极子的数量为N；

（3）利用布置的磁偶极子描述磁场观测值；

磁场观测值序列为T₁,T₂,……,T_M，对于其中任意一个观测点T_i(1≤i≤M)，确定有效区域 包含该观测点的所有磁偶极子序列为p_i1,p_i2,……，p_is，下标表示有效区域包含第i个观测点的 所有磁偶极子的序列，对观测点建立一个包含s个未知量的线性方程，

$T_i=\underset{j=1}{\overset{s}{\Sigma }}{p}_{\mathrm{ij}}·H(x_i,y_i,z_i,{\xi }_{\mathrm{ij}},{\eta }_{\mathrm{ij}},{\zeta }_{\mathrm{ij}},\alpha ,\beta ,\gamma )$

其中，x_i，y_i，z_i表示观测点T的坐标，ξ_ij，η_ij，ζ_ij表示磁偶极子p_ij(1≤j≤s)的坐标值， 根据所有磁场观测值，建立一个共包含M个方程，N个未知量的线性方程组，

$\left(\begin{array}{c}{T}_1=\underset{j=1}{\overset{s1}{\Sigma }}{p}_{1j}·H(x_1,y_1,z_1,{\xi }_{1j},{\eta }_{1j},{\zeta }_{1j},\alpha ,\beta ,\gamma )\\ T_2=\underset{j=1}{\overset{s2}{\Sigma }}{p}_{2j}·H(x_2,y_2,z_2,{\xi }_{2j},{\eta }_{2j},{\zeta }_{2j},\alpha ,\beta ,\gamma )\\ ·\\ ·\\ ·\\ T_M=\underset{j=1}{\overset{\mathrm{sM}}{\Sigma }}{p}_{\mathrm{Mj}}·H(x_M,y_M,z_M,{\xi }_{\mathrm{Mj}},{\eta }_{\mathrm{Mj}},{\zeta }_{\mathrm{Mj}},\alpha ,\beta ,\gamma )\end{array}\right)$

其中，s1,s2,…,sM分别表示相应磁场观测值对应的磁偶极子的个数，M>N；

（4）求解磁偶极子的磁矩；

$\left(\begin{array}{cccc}{H}_{11}& H_{12}& ···& H_{1N}\\ H_{21}& H_{22}& ···& H_{2N}\\ ·& ·& ·& ·\\ ·& ·& ·& ·\\ ·& ·& ·& ·\\ H_{M1}& H_{M2}& ···& H_{\mathrm{MN}}\end{array}\right)\left(\begin{array}{c}{p}_1\\ p_2\\ ·\\ ·\\ ·\\ p_N\end{array}\right)=\left(\begin{array}{c}{T}_1\\ T_2\\ ·\\ ·\\ ·\\ T_M\end{array}\right)$

其中，p₁,p₂,……,p_N表示N个磁偶极子的磁矩，H_i1,H_i2,……,H_iN表示计算得到的第i个 磁场观测点与每个磁偶极子之间的距离，对于p_j(1≤j≤N)，第i个观测点不在其有效区域内， H_ij=0(1≤i≤M,1≤j≤N)；

采用LSQR方法进行求解，得到所有磁偶极子的磁矩；

（5）利用磁偶极子求解观测面上目标网格点的磁场值；

已求得所有磁偶极子的磁矩，根据磁偶极子产生磁场的原理，计算每个磁偶极子在其有 效区域内包含的每个待求网格点上的磁场，将每个网格点上的磁场相加，完成对磁场观测数 据的网格化。

本发明的有益效果在于：

一、通过利用观测磁场反演磁源，然后正演目标位置磁场的方法，遵循了磁场产生和分 布的基本规律，避免了单纯数据处理方法在处理实际物理场数据的不足；

二、通过对磁偶极子在磁场观测平面上设定有效区域，有效的降低了对大量数据处理过 程中的复杂度过高问题，也避免了传统等效源方法中对数据进行分块带来的运算精度降低的 现象。

三、当观测面上的部分区域出现数据量过少或缺失现象时，本方法仍然能够进行该区域 稳定的正演计算，避免了传统方法由于参考数据缺失带来的处理结果误差过大或发散想象。

附图说明

图1是本文方法流程图；

图2是磁偶极子布置示意图；

图3是对于某一磁偶极子在磁场观测平面上在y=0直线上的磁场变化规律图；

图4是对于某一磁偶极子在磁场观测平面上在x=0直线上的磁场变化规律图；

图5是对于某一磁偶极子在整个磁场观测平面上的磁场变化规律图；

图6是磁偶极子在观测平面上的有效区域图。

具体实施方式

下面结合附图1-6对本发明做进一步描述。

根据地磁异常场产生的基本原理，通过在研究区域均匀布置磁偶极子，利用发明的磁偶 极子等效源法对磁矩进行反演，然后通过正演计算目标位置的磁场，从而提出一种基于磁偶 极子等效源法的地磁异常数据网格化方法。

一种基于磁偶极子等效源法的地磁异常数据网格化方法，具体有以下步骤，

步骤一：确定磁偶极子产生磁场的有效区域。

根据磁偶极子产生磁场的基本原理，它在某个观测平面上产生的磁场变化规律明显。对 于每个磁偶极子，在观测平面上设定有效区域，使得磁偶极子在该区域外任意一点的磁场强 度的绝对值小于磁场最大值绝对值的1/20。

步骤二：在磁源平面上布置磁偶极子。

在磁场观测区域水平面下方的某一平面上，以网格形式均匀布置磁偶极子组，并保证， 一：磁偶极子的间隔为观测数据平均间隔的5倍；二：磁偶极子的分布范围应大于观测数据 区域，并使得每个磁偶极子的有效区域内至少包含一个磁场观测值，设磁偶极子的数量为N。

步骤三：利用布置的磁偶极子描述磁场观测值。

设磁场观测序列为T₁,T₂,……,T_M，对于其中任意一个观测点T_i(1≤i≤M)，确定有效区域 包含该观测点的所有磁偶极子，设磁矩分别为p_i1,p_i2,……,p_it,根据磁偶极子产生磁场原理， 可建立一个包含t个未知量的线性方程，进一步，针对所有M个观测点，可建立一个共包含M 个方程，N个未知量的线性方程组。

步骤四：求解磁偶极子的磁矩。

由于观测点数大于布置的磁偶极子数量，因此步骤三中所建方程组是一个超定线性方程 组，同时由于每个方程中只包含有限个未知数，首先在每个方程中将所有未知量补齐，并将 相应的系数值设为零。将方程组写作矩阵形式，并采用LSQR方法求解线性方程组，从而求 得所有磁偶极子的磁矩。

步骤五：利用磁偶极子求解观测面上目标网格点的磁场值。

对于每个目标网格点，确定有效区域包含该网格点的所有磁偶极子，根据磁偶极子产生 磁场原理和已求得的所有磁偶极子的磁矩，求解所有目标网格点的磁场值。

步骤一：确定磁偶极子产生磁场的有效区域。

设磁场观测平面距水平面高度为h₁，确定磁源平面在水平面下方深度为h₀，则磁源与磁 场观测平面的距离为h=h₁+h₂。设磁源平面上一磁偶极子的磁矩大小为p(ξ_i,η_i,ζ_i)，磁矩方 向和当地的磁化方向一致，则该磁偶极子在其上方任意一点(x,y,z)处产生的ΔT磁异常可以 写作，

ΔT(x,y,z)=p(ξ_i,η_i,ζ_i)·H(x,y,z,ξ_i,η_i,ζ_i,α,β,γ)

其中，

H(x,y,z,ξ_i,η_i,ζ_i,α,β,γ)=-U_xi(x,y,z)cosα-U_yi(x,y,z)cosβ-U_zi(x,y,z)cosγ   (2)

$U_{\mathrm{xi}}(x,y,z)=-\frac{1}{4\pi R_i^2}{A}_i(x,y,z),U_{\mathrm{yi}}(x,u,z)=-\frac{1}{4\pi R_i^2}{B}_i(x,y,z),$

$U_{\mathrm{zi}}(x,y,z)=-\frac{1}{4\pi R_i^2}{C}_i(x,y,z)$

A_i(x,y,z)=[2(x-ξ_i)²-(y-η_i)²-(z-ζ_i)²]cosα+3(x-ξ_i)(y-η_i)cosβ+3(x-ξ_i)(z-ζ_i)cosγ

B_i(x,y,z)=3(x-ξ_i)(y-η_i)cosα+[2(y-η_i)²-(x-ξ_i)²-(z-ζ_i)²]cosβ+3(y-η_i)(z-ζ_i)cosγ

C_i(x,y,z)=3(x-ξ_i)(z-ζ_i)cosα+3(y-η_i)(Z-ζ_i)cosβ+[2(z-ζ_i)²-(x-ξ_i)²-(y-η_i)²]cosγ

$R_i=\sqrt{{(x-{\xi }_i)}^2+{(y-{\eta }_i)}^2+{(z-{\zeta }_i)}^2}$

其中，α，β，γ分别为磁化方向与当地地理坐标系X,Y,Z轴的夹角。

设磁源平面上一磁偶极子为p，根据式(1)计算其在磁场观测平面上产生的磁场，若磁偶 极子在某个范围之外任意一点产生的磁场其在观测面上产生磁场最大值的1/20，则这个范围 即为该磁偶极子的有效区域。

设定某一磁偶极子的坐标为（0，0，0），研究区域的磁化方向的余弦值分别为 cosα=0.3714，cosβ=0.5571，cosγ=0.7428，磁场观测平面距离磁源平面高度设为200m， 则磁偶极子在观测平面的对应点坐标为（0，0，200），在观测平面上，在y=0直线上，该磁 偶极子产生磁场的变化规律如附图3所示，在x=0直线上，磁场的变化规律如附图4所示， 综合两个方向上的磁场变化，该磁偶极子在整个磁场观测平面上的磁场变化如附图5所示， 从而可计算有效区域如附图6所示。

步骤二：在磁源平面上布置磁偶极子。

在磁场观测区域水平面下方的某一平面上，以网格形式均匀布置磁偶极子组，并且满足 以下两点，

（1）磁偶极子的间隔为观测数据平均间隔的5倍；

（2）磁偶极子组的分布范围应大于观测数据区域，并使得最外面的磁偶极子的有效区 域内至少包含一个磁场观测值。

设最终布置的磁偶极子的数量为N。

步骤三：利用布置的磁偶极子描述磁场观测值。

设磁场观测值序列为T₁,T₂,……,T_M，对于其中任意一个观测点T_i(1≤i≤M)，确定有效区 域包含该观测点的所有磁偶极子序列为p_i1,p_i2,……,p_is(下标表示有效区域包含第i个观测点 的所有磁偶极子的序列)，根据式(1)，对该观测点可建立一个包含s个未知量的线性方程，

$T_i=\underset{j=1}{\overset{s}{\Sigma }}{p}_{\mathrm{ij}}·H(x_i,y_i,z_i,{\xi }_{\mathrm{ij}},{\eta }_{\mathrm{ij}},{\zeta }_{\mathrm{ij}},\alpha ,\beta ,\gamma )---\left(3\right)$

其中，x_i，y_i，z_i表示观测点T的坐标，ξ_ij，η_ij，ζ_ij表示磁偶极子pi_j(1≤j≤s)的坐标值。 进而，根据所有磁场观测值，可建立一个共包含M个方程，N个未知量的线性方程组，

$\left(\begin{array}{c}{T}_1=\underset{j=1}{\overset{s1}{\Sigma }}{p}_{1j}·H(x_1,y_1,z_1,{\xi }_{1j},{\eta }_{1j},{\zeta }_{1j},\alpha ,\beta ,\gamma )\\ T_2=\underset{j=1}{\overset{s2}{\Sigma }}{p}_{2j}·H(x_2,y_2,z_2,{\xi }_{2j},{\eta }_{2j},{\zeta }_{2j},\alpha ,\beta ,\gamma )\\ ·\\ ·\\ ·\\ T_M=\underset{j=1}{\overset{\mathrm{sM}}{\Sigma }}{p}_{\mathrm{Mj}}·H(x_M,y_M,z_M,{\xi }_{\mathrm{Mj}},{\eta }_{\mathrm{Mj}},{\zeta }_{\mathrm{Mj}},\alpha ,\beta ,\gamma )\end{array}\right)---\left(4\right)$

其中，s1,s2,…,sM分别表示相应磁场观测值对应的磁偶极子的个数，同时由于M>N， 该方程组是一个超定方程组。

此外，每个磁场观测值对应的磁偶极子序列之间是有重复的，为了便于书写和区分，在 方程组(4)的不同方程内，即使是相同的磁偶极子，编号也不相同。

步骤四：求解磁偶极子的磁矩。

首先将方程组(4)写作矩阵方程的形式，如式(5),

$\left(\begin{array}{cccc}{H}_{11}& H_{12}& ···& H_{1N}\\ H_{21}& H_{22}& ···& H_{2N}\\ ·& ·& ·& ·\\ ·& ·& ·& ·\\ ·& ·& ·& ·\\ H_{M1}& H_{M2}& ···& H_{\mathrm{MN}}\end{array}\right)\left(\begin{array}{c}{p}_1\\ p_2\\ ·\\ ·\\ ·\\ p_N\end{array}\right)=\left(\begin{array}{c}{T}_1\\ T_2\\ ·\\ ·\\ ·\\ T_M\end{array}\right)---\left(5\right)$

其中，p₁,p₂,……,p_N表示N个磁偶极子的磁矩，H_i1,H_i2,……,H_iN表示根据公式(2)计算 得到的第i个磁场观测点与每个磁偶极子之间的距离，并且对于p_j(1≤j≤N)，如果第i个观 测点不在其有效区域内，则此时H_ij=0(1≤i≤M,1≤j≤N)。

采用LSQR方法对(5)进行求解，从而得到所有磁偶极子的磁矩。

步骤五：利用磁偶极子求解观测面上目标网格点的磁场值。

已求得所有磁偶极子的磁矩，根据步骤一中磁偶极子产生磁场的原理，计算每个磁偶极 子在其有效区域内包含的每个待求网格点上的磁场，最后，将每个网格点上的磁场相加，则 完成了对磁场观测数据的网格化。