用于地壳弧形构造信息提取的四阶谱矩分析

2020-08-18付丽华杨文采李燕梅

石油地球物理勘探 2020年4期

付丽华曾诚杨文采李燕梅

（①中国地质大学（武汉）数学与物理学院，湖北武汉430074；②西安电子科技大学计算机科学与技术学院，陕西西安710071；③浙江大学地球科学学院，浙江杭州310027；④哈尔滨工业大学数学学院，黑龙江哈尔滨150001）

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］

式中：G（fx，fy）是位场表面z（x，y）的二维功率谱密度；fx、fy分别表示x、y方向上的空间频率。

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

式中：M1、N1分别为滑动窗口在x、y方向包含的点数；F（fu，fv）为z（x，y）的二维离散傅里叶变换，F*（fu，fv）是其共轭。F（fu，fv）的傅里叶逆变换有如下性质

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

式中：m20为表面x方向斜率的方差；m02为y方向斜率的方差；m11是x方向斜率与y方向斜率的协方差。类似地，可推导四阶谱矩元m40、m31、m22、m13和m04的离散空间域表达式

四阶谱矩元中，m40为x方向二阶导数的方差，m04为y方向二阶导数的方差，m31为与的协方差，m13为的协方差，m22为的方差或者的协方差。将计算结果定位于窗口中心点，逐点移动滑动窗口，即可得到每一个数据点的二阶和四阶谱矩元，从而提取到离散位场表面的特征信息。

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

对于二阶谱矩，位场表面可用统计不变量M2和Δ2表征［17，20-21］

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

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

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

类似地，四阶谱矩的统计不变量M4和Δ4的计算公式［20］为

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

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

Λ4的变化区间亦为［0，1］。Λ4与成正比，表示与位场表面圆弧形刻痕呈正相关；与M4成反比，可增强对弱刻痕的识别能力。

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为刻痕的强度系数M2，由式（13）可知水平导数越大其值越大，主要反映较强磁异常的边界刻痕，对弱异常的识别能力较差，且提取的刻痕信息比较发散，不利于边界的准确定位。图1d为脊形化系数Λ2，由图可知，Λ2能均衡识别强弱磁异常，可以确定两球体的中心位置，但是无法准确定位其弧形边界。图1e为表面曲率的方差M4，主要反映较强磁异常的形状与位置信息，对弱异常的识别能力较差。由图1f可知，弧形刻痕系数Λ4能均衡识别强弱磁异常，其高值清晰收敛地刻画出磁源体的边界。图中球体B因磁化强度较小而产生的异常较弱，但在Λ4的提取结果中，它的弧形刻痕信息和球体A的一样清晰，证明了应用四阶谱矩能够提取位场中的弧形刻痕的弱信号。

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

图1 球体模型及谱矩方法不同异常参数提取结果

图2 不同深度球体模型的理论磁异常（左）和弧形刻痕系数Λ4的提取结果（右）

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

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

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

2.2 长方体模型

为检验谱矩方法提取线形刻痕信息的有效性，设计图4a所示包含两个长方体状磁源体（C和D）的组合模型进行理论试验。研究范围为300 km×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可见，应用Λ2能够提取位场中的弱板形刻痕信号，应用Λ4也能够提取位场中的弱板形刻痕信号，不过边界不如Λ2清晰，这是因为Λ2主要反映的是磁异常表面斜率的变化，而Λ4主要反映的是曲率的变化，在长方体模型的平直边界处，磁异常斜率变化剧烈，曲率变化相对较小。仅仅利用Λ4，则可清晰地反映板形体的端角位置，而不是边界。板形体的端角呈现交会线形刻痕的几何特征，能够作为位场中的弧形刻痕信号被提取出来。

图3 常用边界识别方法对球状地壳模型的提取结果（图中黑色虚线是球体边界）

图4 长方体模型及谱矩方法不同异常参数提取结果

3 实际应用

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

图5 塔里木盆地磁异常（a）、Λ2（b）和Λ4（c）分布图

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