计算场源深度是重磁勘探数据解释、定量计算及反演的目标之一，获取较准确的深度信息对重磁勘探具有重要意义.目前计算场源深度的方法可归结为两种，一种以求取场源的深度为主，不求取场源的具体形态，如特征点、切线法、功率谱法(Spector and Grant，1970；侯重初和李保国，1988)、Naudy法(Naudy，1971)、解析信号法(Nabighian，1972，1974)、温纳反褶积法(Ku and Sharp，1983)、欧拉反褶积法(Thompson，1982；Reid et al.，1990)、SPI法或iSPI(improved source parameter imaging)法(Thurston and Smith，1997；Smith et al.，1998)、连续小波变换法(Sailhac et al.，2000；Boukerbout et al.，2003)、tilt-depth法(Salem et al.，2007)等，这类方法已经在重磁勘探解释中得到了广泛应用；另一类方法同时考虑场源的深度及具体形态分布特征，从场源的分布获取深度信息，如物性反演(Li and Oldenburg，1996，1998；Pilkington，1997，2009；管志宁等，1998；Boulanger and Chouteau，2001； Portniaguine and Zhdanov，1999，2002；王妙月等，2004；孟小红等，2012)以及形态反演(Talwani and Ewing，1960；Barnett，1976).物性反演将场源区域划分成具有不同物性的规则单元组合，通过反演求得单元组合的物性参数分布；形态反演以几何形体的形状、角点的坐标等作为求解参数，以几何形体表示目标场源分布.

本文以等效源技术(Dampney，1969)以及物性反演方法为基础，提出一种新的自动求取场源深度方法，通过判断最小反演拟合差求取场源深度.利用位场的多解性特点，在场源区域设置一层规则单元组合体(等效源)，等效源层按照一定间隔由浅至深逐次放置，每一深度反演一次，通过反演求得等效源的物性参数.理论上，等效源的深度越接近实际场源，等效源在观测点上产生的场值与实际场值越吻合.当反演拟合差最小时，等效源的位置达到或接近实际场源位置，由此获取实际场源的深度信息.由于仅需要反演一层等效源，比传统的物性反演计算时间大大减少，并且不需要进行深度加权约束.本文以剖面数据的场源深度计算为例，所述方法可扩展到平面数据.

方法的总体思路是在观测剖面下方某一深度均匀布置一层等效场源，采用物性反演得到等效源物性参数.等效源层以一定的间隔从浅部向深部移动，每移动一次进行固定次数的迭代反演，求得每一深度的反演拟合差.在等效源层由浅变深的过程中，拟合差通常由大变小，直到等效源层到达实际场源中心位置时拟合差最小；当拟合差最小或显著增大时，等效源层停止移动，反演结束，此时等效源层的深度即视为实际场源的深度.可以根据异常特征及实际需求选取不同的移动间隔，设置的间隔越小，等效源向下移动的次数越多，所需的计算时间越长，计算流程见图 1.

图 1

Fig. 1

图 1 深度计算流程Fig. 1 Flow chart of depth estimation

以剖面数据的深度计算为例，将观测面下正演区域剖分成矩形单元，设N_X、N_Z分别为X轴和Z轴上单元个数，正演计算可表示为

下面以理论模型数据进行计算，模型包括不同长宽比例、深度的厚板、薄板，并在部分理论数据中加入了随机噪声.为了便于讨论，本文的理论模型数据算例中皆使用垂直磁化的磁性体.磁异常的正演计算皆采用长方体无解析“奇点”的理论表达式(郭志宏等，2004).

设计的厚板模型形态及位置如图 2所示，形态参数为150 m×150 m×150 m，磁化强度为10 A·m^-1； 模型的中心深度为95 m；观测面离地高度为200 m(以航空磁测数据为例)，点距为25 m.

图 2

Fig. 2

图 2 厚板模型剖面Fig. 2 Profile of thick plate model

根据等效源的设置原则，在每个测点下放置一个宽度为点距的等效源，其横截面的长度与宽度相当，即每个等效源的形态参数都为25 m×20 m×150 m，初始磁化强度为0.01 A·m^-1，反演迭代固定次数为25次.

将反演异常与理论异常求差(见图 3)，结果显示，当等效源底深为94 m及95 m时，差值不大；但当等效源底深移动到96 m(超过模型的中心深度)时，差值显著增大，说明当等效源底深超过实际场源中心深度时，反演异常与理论异常的差值显著增大.理论上，当等效源的深度比实际场源的深度大时，等效源反演无法较好地拟合异常，使反演拟合精度显著变低.本算例中，当等效源的底深为95 m、96 m及105 m时，拟合差分别为0.0004 nT、0.0088 nT及0.3436 nT，这说明等效源的底深比实际场源的中心深度越深，拟合差越大.

图 3

Fig. 3

图 3 反演异常与理论异常的差值Fig. 3 Differences between theoretical anomalies and inversion anomalies at different bottom depth

计算结果显示(见图 4；设置的等效源初始底深为24 m，移动间隔为10 m)，等效源深度由浅至深接近实际场源时，拟合差由大变小，直到等效源的底深超过实际场源的中心深度时，拟合差显著增大.此外，在观测数据中加10%的随机噪声，计算表明该方法对含噪声数据同样适用.

图 4

Fig. 4

图 4 反演拟合差与深度图Fig. 4 Inversion RMS and depth curves

为了说明不同等效源移动间隔对反演结果的影响，采用图 2所示模型的异常曲线进行深度计算.分别选用1 m，5 m，10 m，15 m及20 m的向下移动间隔，计算结果见图 5、表 1(不同移动间隔的突变拐点拟合差对比)，得出以下结论：

图 5

Fig. 5

图 5 不同移动间隔的反演拟合差与深度图Fig. 5 Inversion RMS and depth curves of different intervals

表 1(Table 1)

表 1 不同移动间隔的突变拐点拟合差 Table 1 Fitting errors of break points with different intervals 深度(m) 移动间隔对应的突变拐点拟合差(×10 -4nT) 1 m 5 m 10 m 15 m 20 m 84 8.46 8.46 94 2.77 2.77 95 4.14 96 87.97 99 299.44 299.44 104 2251.43 2251.43 计算深度范围 95~96 94~99 94~104 84~99 84~104 精度 1 5 10 15 20

表 1 不同移动间隔的突变拐点拟合差 Table 1 Fitting errors of break points with different intervals

(1)1，5，10，15 m及20 m拟合差的突变拐点分别为95，94，94，84 m及84 m，计算获得的深度结果范围分别为95~96 m，94~99 m，94~104 m，84~99 m及84~104 m，精度分别为1，5，10，15 m及20 m.说明移动间隔越小，计算的深度精度越高.

(2)5 m及10 m的拟合差突变拐点一致，15 m及20 m的拟合差突变拐点一致，说明等效源移动到同样的深度位置上时的反演拟合差一致.

为了说明该方法对于深部场源同样适用，将该模型放置到中心深度为800 m处.对于深部场源，为了节省计算时间，等效源向下移动的间隔可以适当加大，本算例中采用50 m的向下移动间隔.计算结果显示(图 6)，当等效源的底深超过800 m后，拟合差显著增大，该方法对于深部场源同样适用.

图 6

Fig. 6

图 6 反演拟合差与深度图(模型中心深度800 m)Fig. 6 Inversion RMS and depth curves (model′s central depth is 800 m)

另外计算了宽厚比分别为0.5、0.8、1.2及1.5的厚板状体(见图 7)，等效源的厚度设置为点距，计算结果显示(图 8)，对于厚板状体，计算的深度误差小于1个点距.

图 7

Fig. 7

图 7 不同宽厚比的场源示意图Fig. 7 Models with different ratio of width and thickness

图 8

Fig. 8

图 8 不同厚度场源的反演拟合差与深度曲线Fig. 8 Inversion RMS and depth curves of models with different thickness

水平薄板模型如图 9所示，模型长宽比为7.5，磁化强度为10 A·m^-1；模型的中心深度为95 m； 观测面离地高度为200 m，点距为25 m.分别以10 m、 20 m及30 m厚度的等效源进行计算，递增深度间隔为10 m.结果显示(图 10)，对于该水平薄板，采用10 m厚度等效源的计算结果准确性最高，但采用20 m、30 m厚度等效源的计算结果与实际深度误差都在1个点距左右.

图 9

Fig. 9

图 9 水平薄板模型剖面Fig. 9 Profile of horizontal thin model

图 10

Fig. 10

图 10 不同等效源厚度的反演拟合差与深度曲线(水平薄板模型)Fig. 10 Inversion RMS and depth curves of equivalent sources with different thickness (horizontal thin model)

垂直薄板模型如图 11所示，模型长宽比为7.5，磁化强度为10 A·m^-1；模型的中心深度为95 m； 观测面离地高度为200 m，点距为25 m.分别以20 m、 40 m及60 m厚度的等效源进行计算，递增深度为 10 m.结果显示(图 12)，对于该垂直薄板，采用60 m 厚度等效源的计算结果准确性最高，但采用20 m、40 m厚度等效源的计算结果与实际深度误差都在1个点距左右.

图 11

Fig. 11

图 11 垂直薄板模型剖面Fig. 11 Profile of vertical thin model

图 12

Fig. 12

图 12 不同等效源厚度的反演拟合差与深度曲线(垂直薄板模型)Fig. 12 Inversion RMS and depth curves of equivalent sources with different thickness (vertical thin model)

从水平及垂直薄板状体的算例结果看，等效源的厚度对计算深度的准确度有影响，即等效源的厚度越接近实际场源的厚度，计算结果越精确.但通常情况下，我们并不知道地下异常体的厚度，计算时无法选取与地下场源相符的等效源厚度.当选用与点距相当厚度的等效源时，计算的深度误差为1个点距左右.实际计算时，可选用与点距相当的等效源厚度进行计算.不同的等效源厚度对计算结果的影响主要对于浅部的场源而言(算例中的场源深度都在100 m左右)，对于深部场源，这种影响较小，可忽略不计.

在实际深度反演计算中，在无法确认地下场源形状的前提下，通常采用正方体的等效源，也可采用不同厚度的等效源进行试验，根据计算结果的范围选取合适的深度.

下面讨论组合模型异常的深度计算，为不失一般性，以大小相同(形态参数均为150 m×150 m×150 m)、埋深不同的两个模型引起的异常为例.存在以下两种情况：

(1)浅部模型的磁化强度大于或等于深部模型的磁化强度，则浅部异常大于深部异常，如图 13a所示.该组合模型中深部、浅部模型的磁化强度均为10 A·m^-1，模型的中心深度分别为200 m、150 m；

图 13

Fig. 13

图 13 组合模型异常及深度反演结果 (a)、(b)为不同磁化强度及深度的模型产生的异常；(c)、(d)分别为(a)图、(b)图异常的深度计算结果；(e)、(f)分别为(a)图、(b)图异常深度计算中不同底深的等效源的异常差值. Fig. 13 Combined anomalies and inversion results (a) and (b) Anomalies of two models with different depth and magnetization; (c) Depth estimation result of (a); (d) Depth estimation result of (b); (e) and (f) Corresponding differences between theoretical anomalies and inversion anomalies at different bottom depth.

(2)浅部模型的磁化强度小于深部模型的磁化强度，则深部异常有可能大于浅部异常，如图 13b所示.该组合模型中深部、浅部模型的磁化强度分别为20 A·m^-1、10 A·m^-1，模型的中心深度分别为200 m、150 m.

深度计算时设置等效源的参数与厚板模型算例中的等效源参数一致，向下移动的间隔为15 m.

反演深度计算结果见图 13c、13d，图中反演拟合差的突变拐点深度都为150 m，精确计算得到了浅部模型的中心深度，由于拟合差达到了反演结束的条件，反演停止.将反演异常与理论异常求差，结果显示(图 13e、13f)，当等效源深度超过浅部模型的中心深度时，浅部模型的理论异常与反演异常的差值显著增大，深部模型的差值变化不大.因此在组合模型中，不管浅部模型的异常幅值大于或小于深部异常，本方法都能够较准确地计算得到浅部模型的中心深度，即该方法还可用于识别深浅异常.对于组合模型中深部的场源深度计算，可利用本方法识别浅部异常及计算获得其场源深度后，对异常进行分离后再计算深部场源深度.

对于峰值位置十分接近且叠加的组合模型异常(图 14a)，深度计算拟合差(见图 14b)的突变拐点仍为150 m，同样计算获得了浅部模型的中心深度.

图 14

Fig. 14

图 14 组合模型异常(a)及深度反演结果(b)Fig. 14 Combined anomalies (a) and inversion result (b)

上述理论模型的算例中皆使用了垂直磁化的模型数据，对于斜磁化的磁异常，可以采用两种解决办法，一是在已知磁性体总磁化方向的情况下，可以直接在反演中设置等效源的总磁化方向进行深度计算；二是在未知磁性体总磁化方向的情况下，可将磁异常转化成磁梯度数据，并构建二维解析信号(analytic signal)进行深度计算，二度体磁异常的解析信号模不受斜磁化的影响(Nabighian，1972，1974)，因此可对解析信号模采用本文方法进行深度计算.近年来，随着国内航磁全轴梯度测量的发展(熊盛青，2009；骆遥等，2011，2012；谢汝宽等，2013)，我国已经获取了越来越多的高精度实测航磁全轴梯度数据，可直接用于构建解析信号进行深度计算.

本次实际数据采用2011年在河北某地的航磁全轴梯度数据.航磁全轴梯度测量通过计算磁力仪间的地磁场差值获得沿机身方向的纵向水平梯度、沿机翼方向的横向水平梯度和垂直向下的垂向梯度，反映地磁场在空间三个方向的变化率.本次采用垂向和纵向水平梯度进行深度计算，剖面见图 15.该异常整体上由两个异常组合而成，异常沿测线方向展布约3 km，构建的解析信号剖面图见图 16.利用该解析信号数据进行深度计算，计算时等效源深度向下递增间隔为50 m.深度计算结果见图 17，结果显示随着深度的增加，反演拟合曲线越来越接近实测曲线(图 17a)，反演拟合差越来越小(图 17b、17c).

图 15

Fig. 15

图 15 航磁梯度剖面图Fig. 15 Profile of aeromagnetic gradients

图 16

Fig. 16

图 16 解析信号剖面图Fig. 16 Profile of analytic signal

图 17

Fig. 17

图 17 实测异常与反演异常解释信号剖面图(a),每点的误差曲线(b)及反演拟合差(c)Fig. 17 (a) Measured and inversion analytic signal;(b) inversion fitting error of every measured point;(c) inversion fitting erro

当等效源底深为200 m时，反演拟合差达到最小，至250 m时反演拟合差增大，此时反演拟合曲线已经偏离实测曲线，因此认为计算的磁源中心深度在200~250 m之间.此处已经实施过一口钻井，在钻井岩心的200~300 m处见闪长岩，磁性较大，可能含磁铁矿，磁化率为(12000~25000)×10^-5SI，解释认为该异常由闪长岩磁性体引起，因此计算的深度与实际较吻合.

本文基于等效源及物性反演方法，提出了一种新的自动求取场源深度的方法，并以理论及实际数据进行了算例验证.计算结果表明，等效源的底深超过实际场源的中心深度后，反演拟合差显著增大，据 此可以较准确地获得场源的中心深度.由于仅需要

反演一层等效源，比传统的物性反演计算时间大大减少，并且不需要进行深度加权约束.本方法的计算结果也可以为物性反演、形态反演提供深度信息.计算时可以选用不同的等效源移动间隔，移动间隔越小，计算的深度精度越高.

本方法不仅适用于孤立异常，也适用于叠加了局部异常的组合异常，其计算结果为该组合异常中的浅部异常体的深度，可以进行异常分离后计算深部场源的深度.该方法可以扩展到平面数据的深度计算中，设置的等效源可为一平面，具体计算时可采用滑动窗口等方法选取局部异常进行计算.