0 引言

近年来，强磁性体的退磁效应逐渐成为磁场正、反演领域的研究热点。当磁性体磁化率低于0.1 SI时，可忽略退磁效应^[1]；当磁性体磁化率高于0.1 SI时，退磁效应会显著影响磁异常的幅值、形态等，进而影响磁测数据的处理和解释结果，所以在强磁性体磁场正、反演过程中必须考虑退磁效应的影响。考虑到三维模型的计算量太大，且实际勘探过程中，当地质体沿走向的尺度远大于垂直走向的尺度时(如断层、接触带等)，这类地质体可用走向方向无限延伸的二度体近似，既能满足勘探要求，也能大大提高计算效率^[2]。

目前，对于二度体磁异常的研究多是基于弱磁情况，可忽略退磁效应。Bhattacharyya等^[3-4]、Peder-sen^[5]、吴宣志^[6]给出了频率域不同二度体模型的弱磁异常场表达式；Talwani等^[7]、Shuey等^[8]、Plouff^[9]、Singh等^[10]、Won等^[11]、贾真^[12]提出了空间域二度体数值模拟算法；对于频率域弱磁正演计算方法，柴玉璞^[13]基于偏移抽样提高了反傅里叶变换数值精度；Tontinif等^[14]研究了快速傅里叶变换扩边法与误差的关系；Wu等^[15]在偏移采样的基础上提出了高斯傅里叶变换方法，有效提高了弱磁场数值模拟的精度和效率。

考虑退磁效应的情况下，强磁场数值模拟方法主要有以下三类：积分方程法、有限体积法和有限单元法。Eskola等^[16]基于积分方程提出了考虑退磁效应的静磁方程，其解的准确性与代数方程组的规模有关；Ouyang等^[17]基于强磁场积分方程，采用棱柱体剖分方法实现了强磁性体的正演计算，采用Gauss-FFT方法减小了截断效应；Kostrov^[18]提出了一种采用三角单元的体积分方法求解考虑退磁效应的二维磁场正演问题；付文祥^[19]基于体积分方法提出一种采用三角棱柱单元拟合二度体的方法，并给出了考虑退磁影响的无限长水平圆柱体的总磁场强度异常的理论解析表达式；欧洋等^[20]利用有限体积法实现了同时考虑退磁和剩磁的正演模拟；王书惠^[21]在不引入均匀磁化设定、考虑退磁效应影响的情况下，利用有限元计算了强磁性复杂形体的磁异常；刘双等^[22]利用FlexPDE软件实现了二度圆柱体、板状体的有限元强磁场正演，总结了退磁作用对磁异常幅值影响的规律；另外，刘鑫磊^[23]计算了具有各向异性、非均匀磁化率分布的任意形态地质体的磁场，计算过程考虑了退磁效应，但计算结果可能会不连续；Purssm等^[24]采用迭代方法推导了复杂地形下高磁化率地质体的磁场表达式，但计算结果误差较大。

现有的强磁场数值模拟方法大多基于有限单元法、有限体积法等，这些方法最终都归结为大型线性方程组的求解，对于复杂条件下的大规模二维强磁场数值的模拟存在计算效率低的问题。此外，由重磁位场的研究结论可知，与总磁场相比，梯度张量对异常体边界的识别能力更强、分辨率更高，而现有关于强磁场梯度张量的数值模拟研究较少。因此，本文提出一种适用于大规模复杂条件(复杂磁化率分布和起伏地形)的二度体强磁场及其张量梯度的正演方法——空间—波数混合域二维强磁场及梯度张量数值模拟方法。该方法利用一维傅里叶变换将空间域二维强磁性体磁位满足的二维偏微分方程转换为空间—波数混合域磁位满足的一维常微分方程，利用有限单元法对每个波数下的空间—波数混合域磁位进行计算，并利用追赶法求解定带宽线性方程组，得到空间—波数混合域磁位，利用位场与张量梯度之间的关系得到相应的空间—波数域磁场分量，再经过反傅里叶变换求得空间域场。该算法通过迭代对强磁场进行数值逼近，确保了算法稳定和收敛。

1 方法原理 1.1 控制方程

强磁性体二维磁位方程^[25]为

$ \nabla^2 U(x, z)=\nabla \cdot \boldsymbol{M}(x, z) $

(1)

式中：U为空间域磁位；M =[M_x M_z]^T为空间域磁化强度，其中M_x、M_z分别表示x和z方向的磁化强度；散度算子$\nabla=\boldsymbol{e}_x \frac{\partial}{\partial x}+\boldsymbol{e}_z \frac{\partial}{\partial z}$，其中e_x和e_z分别表示x和z方向的单位矢量。

关于磁介质受到外部磁场的作用，可用磁化强度M衡量磁化程度，它与磁场强度H之间的关系为

$ \boldsymbol{M}=\chi \boldsymbol{H}=\chi\left(\boldsymbol{H}_0+\boldsymbol{H}_{\mathrm{a}}\right) $

(2)

式中：χ表示磁化率；H₀表示背景磁场；H_a表示与区域磁性体相关的磁异常场。对于χ＜0.1 SI的弱磁性体，H_a$ \ll $H₀，因此H_a可忽略不计；对于χ ≥ 0. 1 SI的强磁性体，H_a不可忽略。

对式(1)沿x方向进行一维傅里叶变换，得到空间—波数混合域磁位方程

$ \frac{\partial^2 \widetilde{U}(k, z)}{\partial z^2}-k^2 \widetilde{U}(k, z)=\mathrm{i} k \widetilde{M}_x+\frac{\partial \widetilde{M}_z}{\partial z} $

(3)

式中：k为x方向的波数；$\widetilde{U}$为空间—波数混合域磁位；$\widetilde{M}_x、\widetilde{M}_z$分别为x和z方向上的空间—波数混合域磁化强度。

1.2 边界条件

式(3)的通解为

$ \widetilde{U}=a \mathrm{e}^{k z}+b \mathrm{e}^{-k z} $

(4)

式中a和b为常数。在笛卡尔坐标系中，令z轴向下为正，模型的上边界z=z_min，下边界z=z_max，根据上行波和下行波相关理论，上、下边界条件分别为

$ \left.\frac{\partial \widetilde{U}}{\partial z}\right|_{z=z_{\min }}=k \widetilde{U} $

(5)

$ \left.\frac{\partial \widetilde{U}}{\partial z}\right|_{z=z_{\max }}=-k \widetilde{U} $

(6)

1.3 磁场和磁场梯度张量

对磁位求导可得到相应的磁异常场及其梯度张量，空间域磁异常场B_a、磁异常场梯度张量T和总强度磁异常ΔT^[26]的计算公式分别为

$ \boldsymbol{B}_{\mathrm{a}}=\left(\begin{array}{l} B_{\mathrm{a} x} \\ B_{\mathrm{a} z} \end{array}\right)=\mu\left(\begin{array}{l} H_{\mathrm{a} x} \\ H_{\mathrm{a} z} \end{array}\right)=-\mu\left(\begin{array}{l} \frac{\partial U}{\partial x} \\ \frac{\partial U}{\partial z} \end{array}\right) $

(7)

$ \boldsymbol{T}=\left(\begin{array}{ll} B_{\mathrm{a} x x} & B_{\mathrm{a} x z} \\ B_{\mathrm{a} z x} & B_{\mathrm{a} z z} \end{array}\right)\left(\begin{array}{ll} \frac{\partial B_{\mathrm{a} x}}{\partial x} & \frac{\partial B_{\mathrm{a} x}}{\partial z} \\ \frac{\partial B_{\mathrm{a} z}}{\partial x} & \frac{\partial B_{\mathrm{a} z}}{\partial z} \end{array}\right) $

(8)

$ \Delta T=\left\|\boldsymbol{B}_0+\boldsymbol{B}_{\mathrm{a}}\right\|-\left\|\boldsymbol{B}_0\right\| $

(9)

式中：B_ax、B_az分别表示磁异常场B_a的x和z分量；H_ax、H_az分别表示磁异常场强度的H_a的x和z分量；B_axx、B_axz、B_azx、B_azz表示磁异常场梯度张量的各个分量；μ=μ₀(1+ χ)表示介质磁导率，其中μ₀表示真空磁导率；B₀为背景磁场，本文令B₀= |B₀|= 50000 nT。

由式(7)可得空间—波数混合域磁位$\widetilde{U}$与$\widetilde{\boldsymbol{B}}_{\mathrm{a}}$磁场的关系为

$ \widetilde{\boldsymbol{B}}_{\mathrm{a}}=\left(\begin{array}{l} \widetilde{B}_{\mathrm{ax}} \\ \widetilde{B}_{\mathrm{a} z} \end{array}\right)=-\mu\left(\begin{array}{l} \mathrm{i} k \widetilde{U} \\ k \widetilde{U} \end{array}\right) $

(10)

式中$\widetilde{B}_{\mathrm{a} x}$和$\widetilde{B}_{\mathrm{a} z}$分别表示空间—波数域磁异常场$\widetilde{\boldsymbol{B}}_{\mathrm{a}}$的x和z分量。

由式(8)可得空间—波数混合域磁位$\widetilde{U}$与磁异常场梯度张量$\widetilde{\boldsymbol{T}}$的关系式

$ \widetilde{\boldsymbol{T}}=\left(\begin{array}{ll} \widetilde{B}_{\mathrm{a} x x} & \widetilde{B}_{\mathrm{a} x z} \\ \widetilde{B}_{\mathrm{a} z x} & \widetilde{B}_{\mathrm{a}zz} \end{array}\right)=\mu\left(\begin{array}{cc} k^2 \widetilde{U} & k(-\mathrm{i} k \widetilde{U}) \\ -\mathrm{i} k(k \widetilde{U}) & -k^2 \widetilde{U} \end{array}\right) $

(11)

式中：$\widetilde{B}_{\mathrm{a} x x}=-\widetilde{B}_{\mathrm{a} z z}; \widetilde{B}_{\mathrm{a} x z}=\widetilde{B}_{\mathrm{a} z x}$。

1.4 空间—波数域常微分方程的求解

联立式(3)、式(5)、式(6)，可求解空间—波数混合域磁位的边值问题。

对于二维强磁性体，磁位在空间—波数混合域满足边值问题，可采用基于二次插值的一维有限单元法对其进行求解。该方法保留垂向方向为空间域，可根据实际需求灵活改变垂向网格密度，兼顾了计算精度和计算效率；利用有限元法得到的五对角线性方程组，可利用追赶法^[27]实现方程组的高效、高精度求解。

基于变分原理^[25]，得到与边值问题等价的变分问题

$ \begin{aligned} F(\widetilde{U})= & \int_{z_{\min }}^{z_{\max }}\left(\frac{\partial \widetilde{U}}{\partial z}\right)^2 \mathrm{~d} z+ \\ & 2 \int_{z_{\min }}^{z_{\max }} \widetilde{U}\left(\mathrm{i} k \widetilde{M}_x+\frac{\partial \widetilde{M}_x}{\partial z}\right) \mathrm{d} z+ \\ & \int_{z_{\min }}^{z_{\max }}(k \widetilde{U})^2 \mathrm{~d} z+k\left(\widetilde{U}_{z_{\max }}^2+\widetilde{U}_{z_{\min }}^2\right) \end{aligned} $

(12)

令$\delta F(\widetilde{U})=0$，在笛卡尔坐标系中，沿z方向进行单元剖分，采用二次插值函数，即$\widetilde{U}$在每个单元内是二次变化的。对变分问题(式(12))逐项进行单元分析、总体合成，得到扩展矩阵或列阵

$ F(\widetilde{U})=\boldsymbol{u}^{\mathrm{T}} \boldsymbol{K} \boldsymbol{u}-2 \boldsymbol{u}^{\mathrm{T}} \boldsymbol{P} $

(13)

式中：u为某个波数在z方向各个节点的磁位$\widetilde{U}$组成的列向量；$\boldsymbol{P}=\sum \boldsymbol{P}_{1 q} \boldsymbol{s} \boldsymbol{z}_q$是源项，q表示积分单元；五对角矩阵$\boldsymbol{K}=\sum \boldsymbol{K}_{1 q}+\boldsymbol{b}_1+\boldsymbol{b}_2$。单元矩阵K_1q、b₁、b₂、P_1q和列阵sz_q的求解见附录A。

对F($\widetilde{U}$)求变分，得到

$ \delta F(\widetilde{U})=\delta \boldsymbol{u}^{\mathrm{T}}(2 \boldsymbol{K} \boldsymbol{u}-2 \boldsymbol{P})=0 $

(14)

由于δu ≠ 0，故

$ K \boldsymbol{u}=\boldsymbol{P} $

(15)

基于磁异常场、磁异常场梯度张量与磁位之间的关系(式(10)、式(11))，可得到空间—波数混合域磁异常场$\widetilde{\boldsymbol{B}}_{\mathrm{a}}$和磁异常场梯度张量$\widetilde{\boldsymbol{T}}$，最后利用一维傅里叶反变换，可得空间域的磁场B_a和磁异常场梯度张量T。

1.5 紧算子

求解空间—波数混合域磁位的边值问题时，磁化强度M的表达式(式2)中存在未知项H_a，若直接求解，不能得到五对角线性方程组。本文采用迭代法对真解进行逐次逼近。在电磁法数值模拟中，Torres-Verdín等^[28]提出了一种电磁场扩展的Bron近似法，这是一种基于内部场的新型近似法，通过将背景电场投影到散射张量上对电磁场进行近似表示，该方法适用于异常体体积较大、电阻率差异明显及宽频的情形。Zhdanov等^[29]构建了一个收敛的Bron级数，Gao^[30]进一步修改得到新的收敛Bron近似。Ouyang等^[17]基于前人研究，根据电磁场紧算子推导过程构造了适用于强磁场稳定收敛的迭代格式，其迭代公式为

$ \boldsymbol{H}^{(i+1)}=\frac{2\left[\boldsymbol{H}_0+\boldsymbol{H}_{\mathrm{a}}^{(i+1)}\right]+\chi \boldsymbol{H}^{(i)}}{2+\chi} $

(16)

式中：i表示迭代次数; H⁽ⁱ⁾和H⁽ⁱ⁺¹⁾分别表示第i次和第i+1次迭代得到的总场。

1.6 算法流程

本文提出的适用于大规模复杂条件(复杂磁化率分布和起伏地形)的二度体强磁场及其张量梯度的正演方法，即空间—波数混合域二维强磁场及梯度张量数值模拟方法流程见图 1，具体步骤如下：

(1) 输入背景场磁场强度H₀，即第一次迭代时总磁场强度的初始值H⁽⁰⁾；

(2) 利用式(2)计算空间域磁化强度M，并通过一维傅里叶变换得到空间—波数域磁化强度$\widetilde{\boldsymbol{M}}$；

(3) 利用有限单元法求解空间—波数域磁位的边值问题，得到空间—波数域异常场磁位$\widetilde{U}_{\mathrm{a}}$；

(4) 根据空间—波数域磁异常场强度与异常场磁位的关系求解空间—波数域异常磁场强度$\widetilde{\boldsymbol{H}}_{\mathrm{a}}$，并通过一维反傅里叶变换，得到空间域磁异常强度H_a；

(5) 采用空间域紧算子(式(16))，得到校正后的总磁场H⁽ⁱ⁺¹⁾, 定义两次迭代的强磁总场的相对误差error $=\frac{\left|\sum\right| \boldsymbol{H}^{(i+1)}\left|-\sum\right| \boldsymbol{H}^{(i)}||}{\sum\left|\boldsymbol{H}^{(i+1)}\right|}$。若error < ε(ε为期望值，本文设置为1×10^-4)，则输出H_a，并基于空间域磁异常场B_a计算空间域磁异常场梯度张量T；若error≥ε，则返回第(2)步继续循环计算，直至满足精度要求。

2 算例分析

首先设计一个二维圆柱体模型和一个棱柱体模型分别验证本文算法的正确性和高效性，然后设计一个带起伏地形的、包含一个顺磁性异常体和一个强磁性异常体的模型，验证正演算法对复杂条件的适应性。算例采用串行计算，算法中的傅里叶变换通过Gauss-FFT实现。算法基于Fortran95语言编程，电脑配置为：Intel Core i3-4150 CPU，主频为3.50 GHz，内存为12 GB。

2.1 正确性验证

设计图 2所示二维圆柱体模型，圆柱体沿y方向无限延伸，顶部埋深为300 m。模型的研究区域为：x=[-2000 m，2000 m]，z=[0，1000 m]；网格数为800(x)×400(y)，水平采样间隔为5 m，垂向采样间隔为2.5 m；观测面为z=0。圆柱体磁异常场的x分量B_ax、z分量B_az、梯度张量水平分量∂B_ax/∂x和垂向分量∂B_ax/∂z的解析表达式详见附录B。该模型的总强度磁异常ΔT、磁异常场B_a及其梯度张量T的解析解和数值解及观测面上各点的相对误差分别见图 3、图 4和图 5。

由图 3~图 5可以看出，ΔT的相对误差小于0.5%，B_ax和B_az的相对误差小于0.5%，∂B_ax/∂x和∂B_ax/∂z的相对误差小于0.1%，即ΔT、B_ax、B_az、∂B_ax/∂x和∂B_ax/∂z的解析解和数值解基本吻合，仅在极少数零值点附近相对误差较大。该圆柱体模型的正演结果验证了本算法的正确性。

图 6为该模型的磁异常水平分量B_ax和垂直分量B_az的相对误差随迭代次数的变化曲线。可以看出，经9次迭代后，计算结果即可满足迭代终止条件。

改变圆柱体的磁化率χ (0~100 SI)，统计采用本文算法进行正演的迭代次数，结果见图 7。可以看出，随着χ的增大，所需迭代次数呈近似线性增大趋势。

图 6和图 7表明本文迭代算法适应于不同磁化率的模型，收敛稳定。

2.2 算法效率分析

引用文献[31]中的二度体模型(图 8)，对比本文算法与COMSOL软件^[32]的计算效率。模型包含一个板状异常体，顶面埋深100 m，大小为100 m(x)× 200 m(z); 模拟研究区域为：x=[-500 m，500 m]，z=[0，500 m]，剖分网格数为200×100，水平采样间隔和垂向采样间隔均为5 m。

参照文献[31]中利用有限体积法计算得到的数值解，衡量COMSOL软件和本文算法的精度。引入相对均方根误差rrms^[33]统计整条观测线上总强度磁异常ΔT的相对误差

$ \operatorname{rrms}=\frac{\sqrt{\sum\limits_{i=1}^N\left(X_i-Y_i\right)^2}}{\sqrt{\sum\limits_{i=1}^N Y_i^2}} \times 100 \% $

(17)

式中：N表示观测总点数；X_i表示COMSOL或者本文算法计算得到的第i个观测点的总强度磁异常；Y_i表示有限体积法计算得到的第i个观测点的总强度磁异常。

当网格数为200×100时，采用本文算法的rrms为2.84%，计算时间0.75 s，占用内存12.0 MB。表 1为不同网格剖分方案下基于COMSOL软件的计算精度、耗时和内存占用。基于COMSOL计算时需对研究范围进行扩边处理，当网格扩边4倍时，rrms为2.92%，计算时间618 s，占用内存2387.1 MB。因此，与COMSOL软件相比，在获得相似精度的情况下，本文算法耗时更短，占用内存更少，计算效率更高。

有限体积法、COMSOL软件(扩边4倍)和本文算法的总强度磁异常见图 9a，以有限体积法的数值解为参照，采用COMSOL软件和本文方法计算结果的相对误差见图 9b。可以看出，采用COMSOL软件和空间—波数域方法计算得到的总强度磁异常与有限体积法的计算结果图形均大致吻合。

对于图 8模型，改变模型网格剖分个数，分析本文算法计算时间和占用内存随网格节点数的变化，结果见图 10和表 2。对比表 1与表 2，可见本文算法计算时间更短、占用内存更少，占用内存和计算时间与网格节点数均呈现近似线性正相关变化趋势。

2.3 适应性验证

为了验证本文正演算法对复杂模型的适应性，设计一个起伏地形模型(图 11)，模型同时包含顺磁性异常体和强磁性异常体。模型的背景磁场T₀=50000 nT，磁倾角α=45°，磁偏角β=5°。研究区域：x=[-2000 m，2000 m]，网格数为400；z=[-320 m，820 m]，其中z方向区间[320 m，820 m]为起伏地形最低点以下的计算区域，网格数为100。起伏地形高程满足函数

$ \begin{gathered} f(x)=\frac{1}{5}(2000-|x|) \sin \frac{x}{100 \pi} \\ x \in[-2000, 2000] \end{gathered} $

(18)

基于该模型研究起伏地形条件下插值平面个数对数值解计算精度的影响。起伏地形在z方向的区域为[-320 m，320 m]，该方向采用均匀网格剖分，将插值平面个数分别设为5、9、13、17、21、25、29和33。采用式(17)计算整条起伏观测线上的总强度磁异常ΔT、磁异常场B_a及梯度张量T的解析解。

图 12为该起伏地形模型的总强度磁异常、磁异常场及梯度张量的解析解和数值解的相对均方根误差随插值平面个数的变化曲线。可以看出，总强度磁异常、磁异常场与其梯度张量的rrms均随插值平面个数的增加而降低。采用5个平面进行插值时，起伏地形上总强度磁异常、磁异常场与其梯度张量的rrms均大于1%；当采用25个及以上不同高度平面进行插值时，起伏地形上总强度磁异常、磁异常场与其梯度张量的rrms均低于1%，约为0.9%，说明采用25个插值平面即可达到满意的计算精度。

另外，采用25个插值平面时，耗时约2.330 s，占用内存52.8 MB，表明本文算法计算速度快、占用内存小，适用于复杂地形下二度体磁异常场及其梯度张量的计算，计算效率较高。

采用25个高度平面进行插值，基于本文方法得到的总强度磁异常、磁异常场及其梯度张量分别见图 13、图 14和图 15。可以看出，解析解和数值解的吻合度很高，表明了本文算法对复杂地形模型的计算精度较高。

3 结论

针对强磁性体的退磁效应不能忽略的问题，本文提出了一种空间—波数混合域二度体磁异常场及其梯度张量正演方法。该算法充分利用了傅里叶变换的快速性、追赶法求解的高效性和迭代算法的稳定性的优势，实现了二维强磁性体磁异常场及其梯度张量的高效、高精度数值模拟。设计了水平地形圆柱体模型进行数值试验，通过对比解析解与数值解验证了本算法的正确性和较高的精度；以有限体积法数值解为参照，对比了本文算法和COMSOL软件的计算效率，证明了本算法的高效性；利用起伏地形下同时包含顺磁性异常体和强磁性异常体的复杂模型验证了本算法对复杂地质条件的适应性。另外，在本文算法的计算环节中，正、反傅里叶变换和求解常微分方程这两个环节均有很好的并行性，采用并行方式可进一步提高计算效率。本文为地下二度体磁异常场及其张量梯度的正演提供了一种高效、高精度算法，对于二度体磁异常的反演和人机交互正、反演解释具有重要意义。

附录A 三维强磁场变分问题求解

求解变分问题(式(12))时，可将其写为

$ \begin{aligned} F(\widetilde{U})= & \sum \int_q\left(\frac{\partial \widetilde{U}}{\partial z}\right)^2 \mathrm{~d} z+\sum \int_q(k \widetilde{U})^2 \mathrm{~d} z+ \\ & 2 \sum \int_q \widetilde{J} \widetilde{U} \mathrm{~d} z+k\left(\widetilde{U}_{z_{\max }}^2+\widetilde{U}_{z_{\min }}^2\right) \end{aligned} $

(A-1)

式中

$ \widetilde{J}=\mathrm{i} k \widetilde{M}_x+\frac{\partial \widetilde{M}_x}{\partial z} $

(A-2)

令

$ \left\{\begin{array}{l} \boldsymbol{N}=\left(N_j, N_p, N_m\right)^{\mathrm{T}} \\ \boldsymbol{u}_q=\left(u_j, u_p, u_m\right)^{\mathrm{T}} \end{array}\right. $

(A-3)

式中：j、m分别为单元两端节点坐标；p为单元中点坐标；积分区域q为单元范围。可得

$ \boldsymbol{u}=\boldsymbol{N}^{\mathrm{T}} \boldsymbol{u}_q=\boldsymbol{u}_q^{\mathrm{T}} \boldsymbol{N} $

(A-4)

式(A-1)中第一项单元积分为

$ \begin{aligned} \int_q\left(\frac{\partial \boldsymbol{u}}{\partial z}\right)^2 \mathrm{~d} z & =\boldsymbol{u}_q^{\mathrm{T}}\left(\int_q \frac{\partial \boldsymbol{N}}{\partial z} \frac{\partial \boldsymbol{N}^{\mathrm{T}}}{\partial z} \mathrm{~d} z\right) \boldsymbol{u}_q \\ & =\boldsymbol{u}_q^{\mathrm{T}} \boldsymbol{K}_{1 q} \boldsymbol{u}_q \end{aligned} $

(A-5)

其中

$ \boldsymbol{K}_{1 q}=\int_q \frac{\partial \boldsymbol{N}}{\partial z} \frac{\partial \boldsymbol{N}^{\mathrm{T}}}{\partial z} \mathrm{~d} z=\frac{1}{3 l}\left(\begin{array}{ccc} 7 & -8 & 1 \\ -8 & 16 & -8 \\ 1 & -8 & 7 \end{array}\right) $

(A-6)

式中l为单元长度。

式(A-1)中第二项单元积分为

$ \int_q k^2 \boldsymbol{u}^2 \mathrm{~d} z=\boldsymbol{u}_q^{\mathrm{T}}\left(\int_q k^2 \boldsymbol{N} \boldsymbol{N}^{\mathrm{T}} \mathrm{d} z\right) \boldsymbol{u}_q=\boldsymbol{u}_q^{\mathrm{T}} \boldsymbol{K}_{2 q} \boldsymbol{u}_q $

(A-7)

其中

$ \boldsymbol{K}_{2 q}=\int_q k^2 \boldsymbol{N N}^{\mathrm{T}} \mathrm{d} z=\frac{l k^2}{30}\left(\begin{array}{ccc} 4 & 2 & -1 \\ 2 & 16 & 2 \\ -1 & 2 & 4 \end{array}\right) $

(A-8)

式(A-1)中第三项单元积分为

$ \int_q \widetilde{J} \boldsymbol{u} \mathrm{d} z=\int_q J_{\mathrm{a} z} \boldsymbol{u} \mathrm{d} z $

(A-9)

利用形函数可将上式中的J_az表示为

$ \begin{aligned} J_{\mathrm{a}z} & =J_{\mathrm{a} zj} N_j+J_{\mathrm{a}zp} N_p+J_{\mathrm{a} z m} N_m \\ & =\boldsymbol{N}^{\mathrm{T}} \boldsymbol{s z}_q=\boldsymbol{s z}_q^{\mathrm{T}} \boldsymbol{N} \end{aligned} $

(A-10)

其中

$ \boldsymbol{s} \boldsymbol{z}_q=\left(J_{\mathrm{a} z j}, J_{\mathrm{a} z p}, J_{\mathrm{a}zm}\right)^{\mathrm{T}} $

(A-11)

则式(A-9)中右端项可表示为

$ \int_q J_{\mathrm{a} z} \boldsymbol{u} \mathrm{d} z=\boldsymbol{u}_q^{\mathrm{T}}\left(\int_q \boldsymbol{N} \boldsymbol{N}^{\mathrm{T}} \mathrm{d} z\right) \boldsymbol{s} \boldsymbol{z}_q=\boldsymbol{u}_q^{\mathrm{T}} \boldsymbol{P}_{1 q} \boldsymbol{s} \boldsymbol{z}_q $

(A-12)

其中

$ \boldsymbol{P}_{1 q}=\int_q \boldsymbol{N N}^{\mathrm{T}} \mathrm{d} z=\frac{l}{30}\left(\begin{array}{ccc} 4 & 2 & -1 \\ 2 & 16 & 2 \\ -1 & 2 & 4 \end{array}\right) $

(A-13)

式(A-1)中，z=z_min、z=zmax分别为第一个和最后一个节点，可将其分别扩展为

$ \left\{\begin{array}{l} k u_{z_{\min }}^2=\boldsymbol{u}^{\mathrm{T}} \boldsymbol{b}_1 \boldsymbol{u} \\ k u_{z_{\max }}^2=\boldsymbol{u}^{\mathrm{T}} \boldsymbol{b}_2 \boldsymbol{u} \end{array}\right. $

(A-14)

其中

$ \boldsymbol{b}_1=\left(\begin{array}{llll} k & & & \\ & 0 & & \\ & & \cdots & \\ & & & 0 \end{array}\right) $

(A-15)

$ \boldsymbol{b}_2=\left(\begin{array}{llll} 0 & & & \\ & 0 & & \\ & & \cdots & \\ & & & k \end{array}\right) $

(A-16)

综上所述，式(A-1)可表示为

$ F(\boldsymbol{u})=\boldsymbol{u}^{\mathrm{T}} \boldsymbol{K} \boldsymbol{u}-2 \boldsymbol{u}^{\mathrm{T}} \boldsymbol{P} $

(A-17)

附录B 二维圆柱体强磁场及其梯度张量解析解

假设沿走向无限长圆柱体的磁化率为χ，半径为a，则该圆柱体外的H_a可以表示为

$ \begin{aligned} & \boldsymbol{H}_{\mathrm{a}}= \\ & \frac{C_{\mathrm{m}}}{r^4}\left[\begin{array}{cc} \left(x-x_0\right)^2-\left(z-z_0\right)^2 & 2\left(x-x_0\right)\left(z-z_0\right) \\ 2\left(x-x_0\right)\left(z-z_0\right) & \left(z-z_0\right)^2-\left(x-x_0\right)^2 \end{array}\right] \boldsymbol{H}_0 \end{aligned} $

(B-1)

$ C_{\mathrm{m}}=\frac{\chi}{2+\chi} a^2 $

(B-2)

式中$r=\sqrt{\left(x-x_0\right)^2+\left(z-z_0\right)^2}$表示观测点(x，z)到磁异常体中心(x₀，z₀)的距离。

沿走向无限长圆柱体的磁场梯度张量

$ \boldsymbol{T}=\left(\begin{array}{ll} \frac{\partial B_{\mathrm{a} x}}{\partial x} & \frac{\partial B_{\mathrm{a} x}}{\partial z} \\ \frac{\partial B_{\mathrm{a} z}}{\partial x} & \frac{\partial B_{\mathrm{a} z}}{\partial z} \end{array}\right)=\mu\left(\begin{array}{ll} \frac{\partial H_{\mathrm{a} }x}{\partial x} & \frac{\partial H_{\mathrm{a}x}}{\partial z} \\ \frac{\partial H_{\mathrm{a} z}}{\partial x} & \frac{\partial H_{\mathrm{a} z}}{\partial z} \end{array}\right) $

(B-3)

其中

$ \begin{aligned} \frac{\partial H_{\mathrm{a} x}}{\partial x} & =-\frac{\partial H_{\mathrm{a} z}}{\partial z} \\ & =\frac{2 C_{\mathrm{m}}}{r^6}\left[H_{0 x}\left(x-x_0\right)\left[3\left(z-z_0\right)^2-\left(x-x_0\right)^2\right]+\right. \\ & \left.H_{0 z}\left(z-z_0\right)\left[\left(z-z_0\right)^2-3\left(x-x_0\right)^2\right]\right] \end{aligned} $

(B-4)

$ \begin{aligned} \frac{\partial H_{\mathrm{a}x}}{\partial z} & =\frac{\partial H_{\mathrm{a} z}}{\partial x} \\ & =\frac{2 C_{\mathrm{m}}}{r^6}\left[H_{0 x}\left(z-z_0\right)\left[\left(z-z_0\right)^2-3\left(x-x_0\right)^2\right]+\right. \\ & \left.H_{0 z}\left(x-x_0\right)\left[\left(x-x_0\right)^2-3\left(z-z_0\right)^2\right]\right] \end{aligned} $

(B-5)

式中H_0x、H_0z分别为背景磁场H₀的x分量和z分量。