0 引言

在位场数据解释中，场源边界识别是不可或缺的内容^[1]。事实上，构造单元划分、断裂构造带位置的确定、岩性识别与地层分布解释、物性填图等，均需确定边界位置^[2]。

边界识别方法的发展可分为三个阶段。第一阶段的研究大多基于位场异常表面不同阶次、不同方向上的导数及其组合，包括垂直导数^[3]、水平总导数和二维解析信号^[4]等。由于位场表面的导数随场源深度的增加衰减较快，且处理实际数据时求导计算会放大噪声，所以该阶段方法的主要缺点是对异常幅值较小的场源识别能力较差且抗噪能力差。第二阶段的方法主要针对边界均衡识别，包括倾斜角法(Tilt angle)^[5]、倾斜角的水平总导数法^[6]、自动增益控制法^[7]等。这类方法能同时识别浅部与深部的地质体，但是提取的边界比较发散，定位不够准确。在第三阶段，学者们致力于在均衡识别的基础上增强所提取的边界，同时提高方法的抗噪能力，如Theta map法^[8]、归一化标准差(NSTD)法^[9]、增强型均衡滤波器法^[10]、水平总导数的倾斜角(TAHG)法^[11]等。由此可见，边界识别方法向更均衡、更收敛和更稳定的方向发展。

各阶谱矩及其统计不变量能够详细地描述表面几何特征，早先被广泛应用于工程领域的表面形貌识别，目前已被引入地球物理学领域并取得了良好的应用效果。孙艳云等^[12]和杨文采等^[13-14]基于二阶谱矩定义了脊形化系数和边界脊形化系数，从重力场成功地提取了地壳变形带信息，其中边界脊形化系数在构造边界定位中效果显著。付丽华等^[15]基于二阶谱矩提出了一种描述表面数据粗糙度的特征因子，可有效反映数据波动与变异特征。付丽华等^[16]基于四阶谱矩，应用算术平均顶点曲率提取磁异常数据特征，并将所提取的信息用于场源深度反演，取得了很好的效果。因此，利用谱矩能够从位场数据中提取出与场源边界、埋深等有关的几何特征信息。

相较于其他方法，基于谱矩的位场几何特征分析方法能够更均衡、收敛地识别地质体的几何特征。二阶谱矩可用于提取线形刻痕^[12]，但不能很好地识别其他刻痕。应用四阶谱矩的几何特征分析方法则可以进一步完善对弧形等刻痕信息的识别。为此，提出四阶弧形刻痕系数，理论和真实数据试验证明了利用此方法提取弧形和交会线形刻痕的有效性。最后结合二阶谱矩的脊形化系数对塔里木盆地航磁数据进行特征提取，更为全面地解释了该区域的地壳构造信息。

1 方法原理 1.1 表面谱矩

在笛卡尔直角坐标系中，可用二维随机过程函数z(x, y)表示位场表面，其中z为异常强度，(x, y)是平面直角坐标位置。连续位场表面的p+q阶谱矩定义为^[17-19]

$ {m_{pq}} = \int_{ - \infty }^{ + \infty } {\int_{ - \infty }^{ + \infty } G } ({f_x},{f_y})f_x^pf_y^q{\rm{d}}{f_x}{\rm{d}}{f_y} $

(1)

式中：G(f_x, f_y)是位场表面z(x, y)的二维功率谱密度；f_x、f_y分别表示x、y方向上的空间频率。

实际勘探中得到的位场数据是离散的、有限的。假定研究区域共有M×N个采样点，在x、y方向上分别以Δx、Δy等间距采样，此时位场表面函数可表示为z(x_i, y_j)(i=1, 2, …, M; j=1, 2, …, N)，对这样的离散位场表面，需要用滑动窗口提取局部特征信息。孙艳云等^[12]推导出离散位场表面的p+q阶谱矩的表达式为

$ \begin{array}{l} {m_{pq}} = \sum\limits_{v = 1}^{{N_1}} {\sum\limits_{u = 1}^{{M_1}} G } ({f_u},{f_v})f_u^pf_v^q\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \frac{1}{{{N_1}{M_1}\Delta x\Delta y}}\sum\limits_{v = 1}^{{N_1}} {\sum\limits_{u = 1}^{{M_1}} {[F({f_u},{f_v})f_u^{p/2}f_v^{q/2}] \times } } \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} [{F^*}({f_u},{f_v})f_u^{p/2}f_v^{q/2}] \end{array} $

(2)

式中：M₁、N₁分别为滑动窗口在x、y方向包含的点数；F(f_u, f_v)为z(x, y)的二维离散傅里叶变换，F^*(f_u, f_v)是其共轭。F(f_u, f_v)的傅里叶逆变换有如下性质

$ \begin{array}{*{20}{l}} {{\rm{IDFT}} [F({f_u},{f_v}){f_u}] = \frac{{\partial z({x_u},{y_v})}}{{\partial x}}}\\ {u = 1,2, \cdots ,{M_1};v = 1,2, \cdots ,{N_1}} \end{array} $

(3)

$ \begin{array}{*{20}{l}} {{\rm{IDFT}} [F({f_u},{f_v}){f_v}] = \frac{{\partial z({x_u},{y_v})}}{{\partial y}}}\\ {u = 1,2, \cdots ,{M_1};v = 1,2, \cdots ,{N_1}} \end{array} $

(4)

式中IDFT表示离散傅里叶逆变换。因此窗口面元的二阶谱矩元可用褶积形式表示为

$ \begin{array}{*{20}{l}} {{m_{20}} = {z_x}({x_u},{y_v}) * {z_x}( - {x_u}, - {y_v})}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \sum\limits_{k = 1}^{{N_1}} {\sum\limits_{i = 1}^{{M_1}} {z_x^2} } ({x_i},{y_k})} \end{array} $

(5)

$ {m_{11}} = \sum\limits_{k = 1}^{{N_1}} {\sum\limits_{i = 1}^{{M_1}} {{z_x}} } ({x_i},{y_k}){z_y}({x_i},{y_k}) $

(6)

$ {m_{02}} = \sum\limits_{k = 1}^{{N_1}} {\sum\limits_{i = 1}^{{M_1}} {z_y^2} } ({x_i},{y_k}) $

(7)

且有

$ {{z_x}(x,y) = \frac{{\partial z(x,y)}}{{\partial x}}} $

$ {{z_y}(x,y) = \frac{{\partial z(x,y)}}{{\partial y}}} $

式中：m₂₀为表面x方向斜率的方差；m₀₂为y方向斜率的方差；m₁₁是x方向斜率与y方向斜率的协方差。类似地，可推导四阶谱矩元m₄₀、m₃₁、m₂₂、m₁₃和m₀₄的离散空间域表达式

$ {{m_{40}} = \sum\limits_{k = 1}^{{N_1}} {\sum\limits_{i = 1}^{{M_1}} {z_{xx}^2} } ({x_i},{y_k})} $

(8)

$ {{m_{31}} = \sum\limits_{k = 1}^{{N_1}} {\sum\limits_{i = 1}^{{M_1}} {{z_{xx}}} } ({x_i},{y_k}){z_{xy}}({x_i},{y_k})} $

(9)

$ {{m_{22}} = \sum\limits_{k = 1}^{{N_1}} {\sum\limits_{i = 1}^{{M_1}} {z_{xy}^2} } ({x_i},{y_k})} $

(10)

$ {{m_{13}} = \sum\limits_{k = 1}^{{N_1}} {\sum\limits_{i = 1}^{{M_1}} {{z_{yy}}} } ({x_i},{y_k}){z_{xy}}({x_i},{y_k})} $

(11)

$ {{m_{04}} = \sum\limits_{k = 1}^{{N_1}} {\sum\limits_{i = 1}^{{M_1}} {z_{yy}^2} } ({x_i},{y_k})} $

(12)

式中

$ {{z_{xx}}(x,y) = \frac{{{\partial ^2}z(x,y)}}{{\partial {x^2}}}} $

$ {{z_{xy}}(x,y) = \frac{{{\partial ^2}z(x,y)}}{{\partial x\partial y}}} $

$ {{z_{yy}}(x,y) = \frac{{{\partial ^2}z(x,y)}}{{\partial {y^2}}}} $

四阶谱矩元中，m₄₀为x方向二阶导数$\frac{{{\partial ^2}z}}{{\partial {x^2}}} $的方差，m₀₄为y方向二阶导数$ \frac{{{\partial ^2}z}}{{\partial {y^2}}}$的方差，m₃₁为$\frac{{{\partial ^2}z}}{{\partial {x^2}}} $与$ \frac{{{\partial ^2}z}}{{\partial x\partial y}}$的协方差，m₁₃为$ \frac{{{\partial ^2}z}}{{\partial {y^2}}}$与$ \frac{{{\partial ^2}z}}{{\partial x\partial y}}$的协方差，m₂₂为$ \frac{{{\partial ^2}z}}{{\partial x\partial y}}$的方差或者$\frac{{{\partial ^2}z}}{{\partial {x^2}}} $与$ \frac{{{\partial ^2}z}}{{\partial {y^2}}}$的协方差。将计算结果定位于窗口中心点，逐点移动滑动窗口，即可得到每一个数据点的二阶和四阶谱矩元，从而提取到离散位场表面的特征信息。

1.2 基于谱矩的边界识别方法

对于二阶谱矩，位场表面可用统计不变量M₂和Δ₂表征^{[17, 20-21]}

$ {{M_2} = {m_{20}} + {m_{02}}} $

(13)

$ {{\varDelta _2} = \left| {\begin{array}{*{20}{c}} {{m_{20}}}&{{m_{11}}}\\ {{m_{11}}}&{{m_{02}}} \end{array}} \right| = {m_{20}}{m_{02}} - m_{11}^2} $

(14)

式中M₂表示位场表面斜率的方差，能够反映位场异常的刻痕强弱，因此定义为刻痕的强度系数；Δ₂主要与场在局部区域的各向异性有关，反映表面刻痕脊形化的程度，且通常Δ₂≥0^{[12, 22-23]}。

为从位场表面准确识别地壳变形带，孙艳云等^[12]提出了刻痕的脊形化系数

$ {\varLambda _2} = \frac{{2\sqrt {{\varDelta _2}} }}{{{M_2}}} = \frac{{2\sqrt {{m_{20}}{m_{02}} - m_{11}^2} }}{{{m_{20}} + {m_{02}}}} $

(15)

Λ₂变化区间为[0, 1]，与位场表面各向异性程度呈正相关，与刻痕强度系数M₂成反比，可增强对异常幅值较小的场源的识别能力^{[12, 22]}。

类似地，四阶谱矩的统计不变量M₄和Δ₄的计算公式^[20]为

$ {{M_4} = {m_{40}} + 2{m_{22}} + {m_{04}}} $

(16)

$ {{\varDelta _4} = \left| {\begin{array}{*{20}{l}} {{m_{40}}}&{{m_{31}}}&{{m_{22}}}\\ {{m_{31}}}&{{m_{22}}}&{{m_{13}}}\\ {{m_{22}}}&{{m_{13}}}&{{m_{04}}} \end{array}} \right|} $

(17)

式中：M₄表示表面曲率的方差，可反映表面曲率变化的剧烈程度，数值大说明表面明显地尖锐，峰顶呈尖峰状，数值小说明表面有较小的粗糙度变化，峰顶较圆滑；Δ₄度量表面曲率的方向性效应，对应位场表面呈现圆弧形或者交会线形刻痕的几何特征，值小表明位场表面这种特征不明显。

与二阶谱矩的脊形化系数不同，四阶谱矩反映的是表面曲率变化，可提取圆弧形或者交会线形刻痕的几何特征，这种特征不受位场幅度大小的影响。因此，可以用四阶谱矩提取位场中的弧形刻痕信息，尤其是弧形的弱信号。为此，弧形刻痕系数必须与统计不变量Δ₄正相关，与M₄负相关。根据归一化的要求，用四阶谱矩的统计不变量M₄和Δ₄定义弧形刻痕系数

$ {\varLambda _4} = \frac{{3\sqrt[3]{{{\Delta _4}}}}}{{{M_4}}} $

(18)

Λ₄的变化区间亦为[0, 1]。Λ₄与 $ \sqrt[3]{\mathit{\Delta}_{4}}$ 成正比，表示与位场表面圆弧形刻痕呈正相关；与M₄成反比，可增强对弱刻痕的识别能力。

2 理论模型试验 2.1 球体模型

为检验谱矩方法提取弧形刻痕信息的有效性，设计图 1a所示模型。在大小为120km×120km、网格间距为1km×1km的区域内有两个半径均为15km、中心埋深均为20km的球状磁源体。球体A和B的球心坐标分别为(-20km，-20km)和(20km，20km)，磁化强度分别为0.8A/m和0.3A/m。该模型的理论磁异常见图 1b。

球体在地面的磁异常具有完整的弧形刻痕特征，是应用四阶谱矩提取位场中弧形刻痕信息的理想模型。采用3×3大小的滑动窗口对理论磁异常数据(图 1b)进行特征提取，结果如图 1c~图 1f所示。图 1c为刻痕的强度系数M₂，由式(13)可知水平导数越大其值越大，主要反映较强磁异常的边界刻痕，对弱异常的识别能力较差，且提取的刻痕信息比较发散，不利于边界的准确定位。图 1d为脊形化系数Λ₂，由图可知，Λ₂能均衡识别强弱磁异常，可以确定两球体的中心位置，但是无法准确定位其弧形边界。图 1e为表面曲率的方差M₄，主要反映较强磁异常的形状与位置信息，对弱异常的识别能力较差。由图 1f可知，弧形刻痕系数Λ₄能均衡识别强弱磁异常，其高值清晰收敛地刻画出磁源体的边界。图中球体B因磁化强度较小而产生的异常较弱，但在Λ₄的提取结果中，它的弧形刻痕信息和球体A的一样清晰，证明了应用四阶谱矩能够提取位场中的弧形刻痕的弱信号。

实际情况中地质体的埋深是未知的，因此边界识别方法对深度越不敏感越好，说明地质体埋深对提取效果的影响不大。为进一步检验弧形刻痕系数Λ₄识别弧形边界的有效性，改变图 1a所示模型中球体的埋深，测试Λ₄对深度的敏感程度。图 2a~图 2d是球体中心埋深分别为20、22、24、26km时的理论磁异常和Λ₄提取结果。由于计算结果的范围变化较大，图中根据场值的大小采用不同的色标。

由图 2可知，随着深度增加，原本独立的两磁异常范围渐渐扩大并产生交叠，场源强度逐渐变小，增加了边界识别的难度。分析Λ₄的提取结果可知，弧形分布信息受场源强度影响较小，能稳定有效地反映不同深度地质体的弧形边界信息，对深度的敏感度较低。

为对比弧形刻痕系数Λ₄与常用边界识别方法的效果，本文选取水平总导数(TDX)^[9]、水平总导数的倾斜角TAHG^[11]、Theta map^[8]和归一化标准差(NSTD)方法^[9]对模型(图 1a)进行边界信息提取，结果如图 3a~图 3d所示。

由图 3可知，四种传统方法都能提取到两个球状磁源体的边界，但是均存在边界发散的问题，除TDX外的另外三种方法的提取结果中还存在干扰刻痕，难以准确分辨地质体的边界位置。对比图 1f可知，Λ₄相较于传统方法能更准确、收敛地刻画球状磁源体的边界，且没有干扰刻痕。试验结果表明，Λ₄的边界识别效果明显优于传统边界识别方法。

2.2 长方体模型

为检验谱矩方法提取线形刻痕信息的有效性，设计图 4a所示包含两个长方体状磁源体(C和D)的组合模型进行理论试验。研究范围为300km× 300km、网格间距为1km×1km。磁源体C的中心坐标为(0，5.32km，9.00km)，顶面埋深为4km，长、宽、高分别为150.00、20.22、10.00km，磁化强度为0.5A/m；磁源体D的中心坐标为(0，-79.79km，13.00km)，顶面埋深为8km，长、宽、高分别为179.78、20.21、10.00km，磁化强度为0.2A/m。图 4b为该模型理论磁异常。

据图 4可见，应用Λ₂能够提取位场中的弱板形刻痕信号，应用Λ₄也能够提取位场中的弱板形刻痕信号，不过边界不如Λ₂清晰，这是因为Λ₂主要反映的是磁异常表面斜率的变化，而Λ₄主要反映的是曲率的变化，在长方体模型的平直边界处，磁异常斜率变化剧烈，曲率变化相对较小。仅仅利用Λ₄，则可清晰地反映板形体的端角位置，而不是边界。板形体的端角呈现交会线形刻痕的几何特征，能够作为位场中的弧形刻痕信号被提取出来。

3 实际应用

杨文采等^[22]曾经介绍过塔里木盆地重力场的高阶谱矩分析，现在利用四阶谱矩的弧形刻痕系数分析该地区的航磁场。航空物探数据已经过标定和同化处理^[24]，经度和纬度方向点距分别为0.00886°和0.03380°，采样点数为678×356。磁力化极异常如图 5a所示。采用3×3的滑动窗口对其进行特征提取，Λ₂和Λ₄的提取结果分别见图 5b和图 5c。对比Λ₂与Λ₄的提取结果可见，二者有明显的区别，Λ₂主要反映盆地基底中的线形构造，而Λ₄主要反映出盆地内部主体为等轴状构造，如二叠纪古火山群等。

由于异常叠加及反演固有的多解性，准确解释区域磁异常非常困难，计算磁异常源的埋藏深度可以对地质成果的解释提供更准确的依据^[25-27]。杨文采等^[24]应用三维欧拉反褶积计算了区域磁异常源的埋藏深度。由于克拉通盆地地层通常具有上新下老的规律，将磁异常源分解为三个深度层次，分别为浅层2~5km、中层5~10km和深层10~20km。浅层磁异常源分布在盆地边缘，由中生代构造运动和海西期玄武岩侵位共同作用而成，在图 5a和图 5b中盆地边缘的线性构造中有明显反映。中层磁异常源分布在盆地内部，主要是海西期玄武岩侵位形成的，盆地内部5~10km的海相地层中常有发现。深层磁异常源主要是古老基底变质作用形成的，反映为图 5a和图 5b中盆地南部的线性构造。在磁异常图(图 5a)中，盆地内部深度为5~10km的海相地层中玄武岩侵位围绕古火山群分布，表现为等轴状构造的弱异常，淹没在古老变质基底的深层磁异常中，难以定位。在图 5c中Λ₄强异常表现为红色圆圈状刻痕和圈点(如图中箭头所示)，可能反映古火山群和其他等轴形状的地质体，说明应用四阶谱矩分析和信息提取技术成功地揭示了新的地壳组构。

4 结论

(1) 基于谱矩的位场几何特征分析方法是近些年提出的一种新的自动边界识别方法，相较于传统方法能够更均衡、收敛地识别地质体的几何特征。二阶谱矩可用于提取线形刻痕，但不能很好地识别其他刻痕。应用四阶谱矩的几何特征分析方法可以进一步完善弧形等刻痕信息的识别。

(2) 与刻痕的脊形化系数不同，弧形刻痕系数反映的是圆弧形或者交会线形刻痕的几何特征，这种特征独立于位场的幅度。因此，可以应用四阶谱矩提取位场中的弧形刻痕信息，尤其是弧形的弱信号。

(3) 球体在地面的磁异常具有完整的弧形刻痕特征，是应用四阶谱矩来提取位场中弧形刻痕信息的理想模型。理论模型试验证明了应用四阶谱矩能够提取位场中的弧形刻痕的弱信号。

(4) 在塔里木盆地磁异常图中，盆地内部深度为5~10km的海相地层中玄武岩侵位围绕古火山群分布，表现为等轴状构造的弱异常，淹没在古老变质基底的深层磁异常中，难以定位。弧形刻痕系数强异常可反映古火山群和其他等轴形状的地质体，证明应用四阶谱矩分析和信息提取技术可以揭示特殊的地壳组构。