0 引言

基于岩（矿）石的密度、磁化率、电阻率及波阻抗等物性特点，地球物理方法可分为重、磁、电、震等不同勘探方法。由于观测数据体在空间上的局限性，单凭某一种地球物理方法很难准确推测复杂的地下地质结构^[1]。但不同类型的地球物理数据间具有互补性，可联合多种地球物理数据进行综合解释，降低多解性，提高解释的可靠性^[2-4]。重力场和磁场均是地下所有场源产生异常的总和，具有严重的体积效应，多解性明显。因而，利用重、磁异常进行综合分析与解释，有助于提高地质解释的准确性。

重磁对应分析法是重磁联合定性解释的一种常用方法，该方法基于重磁泊松定理，通过研究重力垂向一阶导数和化极磁异常的相关性开展岩性识别与构造划分^[5]。Garland^[6]首先提出了重磁对应分析方法，应用于计算单一地质体的磁化强度与剩余密度的比值（泊松比），从而确定场源性质；Cordell等^[7]结合重力异常和磁源重力异常振幅谱，估计北极地区的柯尼希斯贝格比（即剩余磁化强度与感应磁化强度之比）及总磁化方向；Chandler等^[8]利用滑动窗口实现了多源重磁场的相关性研究；Chandler等^[9]应用滑动窗口下的重磁对应分析法对美国明尼苏达中东部至佩尼奥克统造山带前寒武系开展地质填图，取得了较好的应用效果。

刘心铸^[10]将重磁对应分析法引入中国，并进行了单一和叠加异常体的模型试验及应用分析；黎益仕等^[11]以重磁对应分析为约束，开展了重磁联合反演，获得了良好的应用效果；范正国等^[12]分析了重磁对应分析法的局限性，提出了自适应重磁对应分析法，在一定程度上提高了方法的可靠性与实用性；宋景明^[13]利用磁源重力异常剔除磁性体引起的重力异常，得到了反映沉积构造特征的剩余重力异常；在此基础上，颜廷杰等^[14]利用多尺度窗口滑动线性回归修正方法将磁场换算为磁源重力异常，得到了纯粹由无磁性地质体产生的剩余重力异常；王彦国等^[15]基于重磁泊松定理推导了总磁化方向的估计方法，并在相山铀矿田获得了良好的应用效果。一些学者还在不同地区开展了重磁对应分析法的应用，取得了很好的应用效果^[16-19]。然而，常规重磁对应分析法需要进行化极处理，且仅适用于重磁异常高度相关的情况。此外，该方法易产生虚假信息，严重限制了该方法的实际应用^{[12, 18]}。鉴于此，本文从重力二阶垂向导数与归一化磁源强度的表达式出发，推导出新形式下的重磁泊松公式。在此基础上，提出了基于重力二阶垂向导数与归一化磁源强度互相关系数的重磁相关性分析方法，并通过噪声扰动处理消除重磁同源体之外的虚假强相关信号。由于本方法使用了重力二阶导数，很大程度上削弱了区域重力场的影响；同时，归一化磁源强度基本不受磁化方向影响。因此，本文方法无需做化极处理，具有较强的适用性和实用性。通过模型试验验证了方法的可行性与优越性。最后，将该方法应用于相山铀矿田重磁资料处理，获得了良好的效果。

1 基本原理

典型地质体的重力异常可表示为^[20-21]

$ g=\frac{G\mathrm{\Delta }\rho {k}_{\mathrm{g}}({z}_{0}-z)}{{\left[(x-{x}_{0}{)}^{2}+(y-{y}_{0}{)}^{2}+(z-{z}_{0}{)}^{2}\right]}^{({N}_{\mathrm{g}}+1)/2}} $

(1)

式中：G是万有引力常数；$ \mathrm{\Delta }\rho $为地质体与围岩的密度差；k_g是与坐标位置无关但与地质体几何形状有关的参数；$ (x, y, z) $和$ ({x}_{0}, {y}_{0}, {z}_{0}) $分别是计算点和场源点坐标；$ {N}_{\mathrm{g}} $是重力源构造指数，与地质体几何形状有关。

对重力异常g求取z方向的二阶导数

$ \begin{array}{l}{g}_{zz}=G\mathrm{\Delta }\rho {k}_{\mathrm{g}}({N}_{\mathrm{g}}+1)({z}_{0}-z)\times \\ \frac{{N}_{\mathrm{g}}({z}_{0}{-z)}^{2}-3\left[(x-{x}_{0}{)}^{2}+(y-{y}_{0}{)}^{2}\right]}{{\left[(x-{x}_{0}{)}^{2}+(y-{y}_{0}{)}^{2}+(z-{z}_{0}{)}^{2}\right]}^{({N}_{\mathrm{g}}+5)/2}}\end{array} $

(2)

典型地质体的磁场可用归一化磁源强度统一表示为^[22-23]

$ \begin{array}{l}{T}_{\mu }=\sqrt{-{\lambda }_{2}^{2}-{\lambda }_{1}{\lambda }_{3}}\\ =\frac{{\mu }_{0}M{k}_{\mathrm{M}}}{4\mathrm{\pi }{\left[(x-{x}_{0}{)}^{2}+(y-{y}_{0}{)}^{2}+(z-{z}_{0}{)}^{2}\right]}^{({N}_{\mathrm{M}}+1)/2}}\end{array} $

(3)

式中：$ {\mu }_{0} $是真空中的磁导率；M是磁源磁化强度；k_M是与坐标无关但与场源形状相关的参数；$ {N}_{\mathrm{M}} $是磁源构造指数，与场源形状有关；$ {\lambda }_{1} $、$ {\lambda }_{2} $和$ {\lambda }_{3} $是磁梯度张量$ \boldsymbol{\varGamma } $的特征值。

磁梯度张量$ \boldsymbol{\varGamma } $^[24]可表示为

$ \boldsymbol{\varGamma }={\mu }_{0}\left[\begin{array}{ccc}\frac{{\partial }^{2}U}{\partial {x}^{2}}& \frac{{\partial }^{2}U}{\partial x\partial y}& \frac{{\partial }^{2}U}{\partial x\partial z}\\ \frac{{\partial }^{2}U}{\partial y\partial x}& \frac{{\partial }^{2}U}{\partial {y}^{2}}& \frac{{\partial }^{2}U}{\partial y\partial z}\\ \frac{{\partial }^{2}U}{\partial z\partial x}& \frac{{\partial }^{2}U}{\partial z\partial y}& \frac{{\partial }^{2}U}{\partial {z}^{2}}\end{array}\right]=\left[\begin{array}{ccc}\frac{\partial {B}_{x}}{\partial x}& \frac{\partial {B}_{x}}{\partial y}& \frac{\partial {B}_{x}}{\partial z}\\ \frac{\partial {B}_{y}}{\partial x}& \frac{\partial {B}_{y}}{\partial y}& \frac{\partial {B}_{y}}{\partial z}\\ \frac{\partial {B}_{z}}{\partial x}& \frac{\partial {B}_{z}}{\partial y}& \frac{\partial {B}_{z}}{\partial z}\end{array}\right] $

(4)

式中：U表示磁位；$ {B}_{x} $、$ {B}_{y} $、$ {B}_{z} $分别表示磁场B的x、y、z三个方向分量，其值可利用波谱关系由磁异常T换算得到。

当观测点位于场源上方时，重力二阶垂向导数和归一化磁源强度公式可分别简化为

$ \left\{\begin{array}{l}{g}_{zz}({x}_{0}, {y}_{0}, z)=G\mathrm{\Delta }\rho {k}_{\mathrm{g}}\frac{{N}_{\mathrm{g}}({N}_{\mathrm{g}}+1)}{({z}_{0}{-z)}^{{N}_{\mathrm{g}}+2}}\\ {T}_{\mu }({x}_{0}, {y}_{0}, z)=\frac{{\mu }_{0}M{k}_{\mathrm{M}}}{4\mathrm{\pi }}\frac{1}{({z}_{0}{-z)}^{{N}_{\mathrm{M}}+1}}\end{array}\right. $

(5)

从常规重磁泊松公式可知，对于重磁同源体，其磁异常与重力异常的一阶导数形式一致，因此$ {N}_{M}={N}_{\mathrm{g}}+1 $，式（5）可转换为

$ {T}_{\mu }({x}_{0}, {y}_{0}, z)=\frac{{\mu }_{0}{k}_{\mathrm{M}}}{4\mathrm{\pi }G{k}_{\mathrm{g}}{N}_{\mathrm{g}}({N}_{\mathrm{g}}+1)}\cdot \frac{M}{\Delta \rho }\cdot {g}_{zz}({x}_{0}, {y}_{0}, z)\\ =k\frac{M}{\Delta \rho }\cdot {g}_{zz}({x}_{0}, {y}_{0}, z) $

(6)

上式即为基于重力二阶垂向导数与归一化磁源强度的重磁泊松公式。对于同一场源，记$ k=\frac{{\mu }_{0}{k}_{\mathrm{M}}}{4\mathrm{\pi }G{k}_{\mathrm{g}}{N}_{\mathrm{g}}({N}_{\mathrm{g}}+1)} $为一常数，此参数仅与场源的几何形状相关。

依据式（6），可引入下面的互相关系数判断重磁的同源性

$ C=\frac{\sum\limits_{i=1}^{W}{g}_{zz}({x}_{i}, {y}_{i}, z){T}_{\mu }({x}_{i}, {y}_{i}, z)}{\sqrt{\sum\limits_{i=1}^{W}{g}_{zz}^{2}({x}_{i}, {y}_{i}, z)\times \sum\limits_{i=1}^{W}{T}_{\mu }^{2}({x}_{i}, {y}_{i}, z)}} $

(7)

式中W表示滑动计算窗口的总点数。上式中包含重力二阶导数和磁场一阶导数，可大幅度降低重磁背景场的影响。

根据柯西不等式可得

$ \begin{array}{l} {\left[\sum\limits_{i=1}^{{W}_{n}}{g}_{zz}({x}_{i}, {y}_{i}, z){T}_{\mu }({x}_{i}, {y}_{i}, z)\right]}^{2}\le \\ \sum\limits_{i=1}^{{W}_{n}}{g}_{zz}^{2}({x}_{i}, {y}_{i}, z)\sum\limits_{i=1}^{{W}_{n}}{T}_{\mu }^{2}({x}_{i}, {y}_{i}, z)\end{array} $

(8)

因此，相关系数C取值范围为[-1, 1]。

同样，依据式（6）可以计算包含场源几何形状影响因素的重磁泊松比

$ k\frac{M}{\mathrm{\Delta }\rho }=\frac{\sum\limits_{i=1}^{W}{T}_{\mu }({x}_{i}, {y}_{i}, z)}{\sum\limits_{i=1}^{W}{g}_{zz}({x}_{i}, {y}_{i}, z)} $

(9)

直接使用重力二阶垂向导数和归一化磁源强度进行互相关计算时，效果并不理想，主要原因是远离场源时，归一化磁源强度太小且窗口内数值变化不明显，此时窗口内各点的$ {T}_{\mu } $趋于窗口内的平均值$ \overline{{T}_{\mu }} $，则式（7）可写作

$ \begin{array}{l}C=\frac{\overline{{T}_{\mu }^{}}\sum\limits_{i=1}^{{W}_{n}}{g}_{zz}({x}_{i}, {y}_{i}, z)}{\sqrt{\sum\limits_{i=1}^{{W}_{n}}{g}_{zz}^{2}({x}_{i}, {y}_{i}, z){\left(\overline{{T}_{\mu }^{}}\right)}^{2}}}\\ =\frac{\sum\limits_{i=1}^{{W}_{n}}{g}_{zz}({x}_{i}, {y}_{i}, z)}{\sqrt{\sum\limits_{i=1}^{{W}_{n}}{g}_{zz}^{2}({x}_{i}, {y}_{i}, z)}}\to \pm 1\end{array} $

(10)

根据上式可知，当远离磁源区且重力二阶垂向导数处于非零值区域时，系数C接近$ \pm 1 $，即完全正相关或负相关，显然这是不合理的。事实上，常规重磁分析方法同样存在这样的问题^[12]，目前对此尚未有理想的解决方案。

为解决这一问题，对重力二阶垂向导数和归一化磁源强度添加两个不相关的噪声干扰，然后再进行互相关计算。虽然添加随机噪声会降低数据信噪比，但互相关系数计算是一个平均化处理过程，可有效削弱噪声对有效信号的干扰。在接近场源区内，场源的有用信号占比更大，互相关系数主要取决于重磁异常的相关性；在远离场源处，噪声占主导地位，且随机分布，此时互相关系数趋于零。由于归一化磁源强度具有非负性，因此互相关系数趋于1时，意味着高密度、强磁性场源的存在；接近-1时，表明场源具有低密度、强磁性的特点。

2 模型试验 2.1 单一异常体模型试验

为了测试方法的可行性与有效性，首先建立一个单一异常体模型，分别模拟重磁完全同源、部分同源及完全不同源三种情况。

模型包含一个正方体异常体，异常体边长为2 km，上顶面埋深为1 km，剩余密度为1.0 g/cm³，磁化强度为1 A/m，磁倾角及磁偏角均为45°。图 1为该模型的重力异常与磁异常。

图 2是该模型的理论重力异常一阶垂向导数和理论垂直磁化磁异常，及含10%随机干扰的异常，以此数据计算常规重磁对应分析法的相关系数。需要指出的是，该模型是重磁完全同源，与重磁不同源的差别是仅异常位置发生偏移，异常形态是一致的，因此文中不对其进行展示。

图 3（上）是基于图 2（上）数据，采用常规重磁分析法分别对重磁异常同源、部分同源、完全不同源三种情况下计算的系数C。可以看出，重磁异常完全同源时，C在整个测区都接近于1，这不利于确定重、磁源的真实位置；当重磁异常部分同源和完全不同源时，C在远离场源区存在大面积的正相关和负相关，在重力源与磁源质心连线上，则存在明显的负相关。这说明无论是重磁同源还是不同源，理论重力一阶垂向导数和垂直磁化磁异常的重磁对应分析法效果并不理想，均不能较准确地反映重磁场的相关性。范正国等^[12]对常规重磁对应分析法开展了详细的模拟实验工作，得到了同样的结论，这也是该方法未能广泛应用的主要原因。

图 3（下）是基于图 2（下）数据，分别对重磁异常同源、部分同源、完全不同源三种情况计算的C。可以看出，由于噪声干扰的引入，C值在远离场源区处于零值上下波动。在重磁完全同源的场源上方接近1，但重磁部分同源及完全不同源时，重磁源质心连线上仍存在明显的负相关。这表明，对重、磁场进行噪声扰动有助于提高重磁对应分析方法在重磁完全同源时的可靠性，但在其他情况（重磁部分同源及完全不同源）下仍不适用，易得到虚假信息，不利于重磁综合解释。

图 4是基于图 1中的重力异常和斜磁化磁异常换算得到的重力二阶垂向导数和归一化磁源强度，以及分别添加10%随机干扰的结果。图 5（上）是对图 4（上）数据分别计算重磁同源、部分同源及完全不同源时的无噪声重力二阶垂向导数与归一化磁源强度的互相关系数C。可以看出，无论重磁是否同源，在重力源内外均存在强相关性，这显然是不合理的。图 5（下）是基于图 4（下）计算的重力二阶垂向导数与归一化磁源强度的互相关系数C。可以看出：添加噪声干扰后，C在重磁同源位置上相关性较强，而在重磁场源外则相关性较弱；当重磁完全不同源时，整个测区内的相关系数基本在零值附近，即重磁完全不相关。显然，通过添加噪声扰动可显著提高方法的有效性与可靠性。相对于常规重磁对应分析法，本文方法能够更好地反映重磁同源体的真实位置。

2.2 叠加模型试验

为了验证方法对复杂情况的应用效果，构建了一个多场源组成的叠加模型，模型参数见表 1。其中正方体与岩脉为重磁同源，棱柱体为部分同源，球体是无磁性的。由于常规重磁对应分析法需要将斜磁化磁异常转换成垂直磁化磁异常，而对于磁化方向不同的多场源，难以进行有效的化极处理，因此这里不再进行常规算法的试验。

图 6a、图 6b是组合模型产生的重力异常和磁异常。可以看出，受规模较大的棱柱体异常的影响，岩脉上方的重力异常表现为等值线同行扭曲；受磁化方向影响，磁异常较复杂，难以直接识别磁源位置。

图 6c、图 6d是无噪声时重力二阶垂向导数及归一化磁源强度。可以看出，重力二阶垂向导数和归一化磁源强的极值均与自身场源位置对应很好。图 6e、图 6f是添加10%随机干扰后的重力二阶垂向导数和归一化磁源强度。

图 7a是利用无噪声时重力二阶垂向导数与归一化磁源强度（图 6c和图 6d）计算的互相关系数。可以看出，在重力场源分布区内存在明显的相关性，在场源外也呈现强相关性，尤其值得注意的是，球体本身是无磁性的，但也反映出明显的负相关，即表现为负密度、强磁性的特征，这显然会给解释带来不利影响。图 7b则是利用噪声扰动后重力二阶垂向导数与归一化磁源强度（图 6e和图 6f）计算的互相关系数。可以看出，互相关系数在正方体上方为明显的负相关（接近-1），在岩脉和棱柱体同源部位表现为强正相关（接近1），在球体上方则显示为完全不相关（接近0），除此之外的区域，互相关系数同样在零值附近波动。这也再次表明在进行重磁相关性计算时，添加噪声扰动有助于准确判断重磁异常是否同源。

3 应用实例

为了检验算法对实际资料的应用效果，选取了江西相山铀矿田地区的地面重、磁数据（比例尺1：50000）进行测试。江西相山火山盆地是中国第一大、世界第三大火山岩型铀矿田^[25]。

相山火山盆地位于华夏板块湘桂赣地块北东缘乐安—抚州断隆带上（图 8）。该区遭受了扬子—加里东期、海西—印支期造山作用，燕山期NEE向赣杭构造火山岩带西南端与近SN向赣中南花岗岩带的交接地带发生了强烈的构造—岩浆—成矿作用^[26-27]。相山铀矿田中心区域出露大面积的下白垩统鹅湖岭组碎斑熔岩，向外依次分布着条带状的下白垩统打鼓顶组流纹英安岩、下白垩统花岗斑岩及青白口系变质岩，西北侧还发育上白垩统的沉积岩。相山地区铀矿化主要受断裂构造、火山构造、变质基底及岩石层间接触面控制，多控矿因素交汇更有利于富矿的形成^[28-31]。因此，推断有利控矿因素的分布位置对相山地区铀矿勘查具有重要的研究意义。

为方便重、磁异常解释及其相关性分析，图 9给出了相山火山盆地主要岩石的密度、磁化率散点分布，其中碎斑熔岩、流纹英安岩、花岗斑岩及变质岩样品分别为822、353、101、75块。从图 9可以看出，碎斑熔岩和花岗斑岩密度均主要为2.60~2.65 g/cm³，均值分别为2.63、2.64 g/cm³，同属于低密度体；流纹英安岩密度主要为2.65~2.75 g/cm³，均值为2.69 g/cm³，为中—高密度体；变质岩密度主要为2.7~2.8 g/cm³，均值为2.76 g/cm³，为高密度体。部分碎斑熔岩岩样的磁化率为（10~100）×10^-6 SI，属于低磁性，但大部分岩样的磁化率分布在（100~1000）×10^-6 SI，相当于中强磁性；流纹英安岩岩样磁化率则多数为（10~100）×10^-6 SI，部分为（100~1000）×10^-6 SI，还有少量样本磁化率大于1000×10^-6 SI，即流纹英安岩磁性变化较大，弱、中、强磁性样品均有；花岗斑岩约一半样品磁化率为（10~100）×10^-6 SI，另一半为（100~1000）×10^-6 SI，说明花岗岩表现为弱、中磁性；变质岩的磁化率分布较集中，多数为（20~30）×10^-6 SI，表现为弱磁性。上述岩石密度、磁化率统计分析结果表明，碎斑熔岩、中—强度磁性花岗斑岩及中—强磁性流纹英安岩这三类岩石都可看作是重、磁同源体。

图 10是相山铀矿田的布格重力异常和总磁力异常图。布格重力异常整体表现为中心区域异常低、外围异常高，这样的分布特征主要与大面积的低密度碎斑熔岩和高密度变质岩有关。磁异常主要以局部异常为主，正磁异常主要分布在研究区的南侧与中东侧，以SN、NE、NEE向为主，这些正磁异常在近正北侧都可以寻找到相对应的负磁异常分布。然而，在邹家山—邹家一带，存在大面积的负磁异常，且这些负磁异常幅值明显大于周围的正异常幅值，表明该区域的岩石有可能具有与当代地磁场方向不一致的剩磁^[15]。

图 11（上）是研究区观测面上的重力二阶垂向导数和归一化磁源强度。可以看出，地面上的重力二阶垂向导数分布凌乱，以小规模条带状异常为主，各个方向异常相互切割严重；地面归一化磁源强度主要在中西部、中南部及西北部存在明显的异常，多数铀矿点位于异常高—低分界位置。

图 11（下）是向上延拓500 m的重力二阶垂向导数和归一化磁源强度。可见此时重力导数及归一化磁源强度的局部异常规模较大，走向明显，均以NS、NE、NEE向为主。铀矿点主要位于重力导数梯度带（或零值线）附近，及归一化磁源强度异常高—低分界线位置这些位置往往与断裂带或岩性接触带关系密切。

图 12是基于添加了10%噪声扰动下相山地区重力二阶垂向导数与归一化磁源强度计算的互相关系数。可以看出，地面上互相关系数C（图 12a）仅在牛头岭北侧存在一个呈近东西向展布的条带状高值区，该区对应大面积出露的流纹英安岩，其他位置基本都在-0.5~0.5，即大部分地区重磁相关性并不明显，这是由于对重力场进行了二阶导数计算，因而异常较为凌乱。向上延拓500 m后的互相关系数C（图 12b）的分布反映出多个重磁强相关区域：①牛头岭北侧、淮头南侧这两个正强相关区域反映了高密度、强磁性的流纹英安岩具有一定规模；②云际—尧岗西沿线、芙蓉山西侧、居隆庵西侧存在强正相关区，推测为碎斑熔岩之下存在一定规模的流纹英安岩，其中居隆庵地区的大量钻孔资料已证实了这一推测^[32-33]，这为相山地区在层间接触面寻找铀矿提供了参考，尤其相山铀矿田西部的云际—尧岗西沿线一带；③牛头岭、居隆庵、相山、尧岗北等地区存在大面积负强相关区域，应是低密度、强磁性碎斑熔岩的反映，何家—戴家村一带存在近东西向的条带状负相关区域，应与低密度的花岗斑岩有关，不过何家附近的负相关性更加显著，区域更大，表明这里的花岗斑岩应具有较大规模，而位于相山铀矿田北部的沙洲矿床产于花岗岩内外接触带^[34-35]，因此重磁相关性分析结果为沙洲矿床开展深部铀矿勘查也提供了依据；④云际东、凤岗南及芙蓉山等区域的强负相关性区域，可能反映了低密度、强磁性的花岗斑岩分布；⑤几乎所有的铀矿点都位于强正相关与强负相关的过渡带，反映了相山地区铀矿点与岩性接触面关系密切；⑥牛头岭、芙蓉山、相山、尧岗北分布的古火山口处于大面积的强负相关区域，表明这些火山口喷发的主要物质是碎斑熔岩^[36]，而牛头岭北、淮头南两地的火山口则以流纹英安岩喷出为主，另外居隆庵东侧大面积的负强相关则表明该区可能还存在一个碎斑熔岩的古火山口。显然，这些丰富的重磁相关性信息为相山铀矿田地质研究及铀矿勘查提供了一定参考依据。

4 结论

基于重力二阶垂向导数与归一化磁源强度表达式关系，本文提出了一种新的重磁相关分析方法，削弱了背景场与磁化方向的影响。为提高相关分析结果的可靠性，本文对重力二阶垂向导数和归一化磁源强度进行较强的随机噪声扰动。模型试验表明，此方法可以有效地识别重磁同源体的位置，且不存在虚假信息；同时，该方法适用于重磁场源位置不同、磁源磁化方向不一致的复杂环境。将本文方法应用于相山铀矿田，重磁相关性分析结果较好地揭示了鹅湖岭组碎斑熔岩、打鼓顶组流纹英安岩及花岗斑岩等重磁同源岩石的分布位置及规模，也反映了不同位置古火山喷出的主要岩石类型，同时揭示了铀矿与岩性接触面的密切关系。这些结论为研究相山地区深部铀矿成矿环境提供了数据基础。