CN110414060B - 一种基于四阶谱矩的位场边界识别方法 - Google Patents

Abstract

本发明公开了一种基于四阶谱矩的位场边界识别方法，首先通过滑动窗口计算出位场每一个数据点的四阶谱矩元，然后根据本发明定义的四阶曲率弧刻痕系数公式计算其特征提取结果，最后由该系数的高值指示出场源的边界位置。本发明提出的基于四阶谱矩的位场边界识别方法，是一种适用于提取弧形边界的方法，克服了常规传统方法和表面二阶谱矩不能有效识别弧形边界的缺陷。重磁异常的仿真试验结果表明，本发明提出的基于四阶谱矩的位场边界识别方法能够均衡收敛地识别弧形边界，且对深度的敏感性较低，可用于提取位场中的弱曲率弧刻痕信号。

Description

一种基于四阶谱矩的位场边界识别方法

技术领域

本发明属于地球物理学位场数据处理技术领域，具体涉及一种基于四阶谱矩的位场边界识别方法。

背景技术

揭示地下地质体的边缘位置是位场数据解释中的重要内容。各地质体边缘是具有一定密度或磁性差异的，在地质体的边缘附近，重磁异常的变化率较大，因此可以利用这一特点对其进行边界识别。目前，存在三大类位场源边界识别方法，它们分别是数值类计算方法、数理统计类计算方法及其它计算方法。其中，数值类计算方法是现阶段研究最多的一类方法，且其应用也最为广泛，数理统计类计算方法次之，其他计算方法应用较少。

数值类计算方法大多基于位场异常不同方向、不同阶次的导数及其组合，如总水平导数(TotalHorizontal Derivative,THDR)、垂向导数(Vertical Derivative,VDR)、解析信号振幅(Analytical Signal Amplitude,ASM)、倾斜角(Tilt Angle,TA)和Theta map。数理统计类方法主要包括两种，小子域滤波和归一化标准差(Normalized StandardDeviation,NSTD)，该类方法通过选取一定尺寸大小的滑动窗口来进行计算。上述方法均属于传统的边界识别方法，其中，目前比较常用的方法有以下几种：

总水平导数法(THDR)是一种应用广泛的位场边界识别方法，计算公式为

其中，u表示位场异常。

解析信号振幅(ASM)是由Nabighian提出的，该方法也叫总梯度模量法，由三个方向上的一阶导数组合得到，表达式为

式中，THDR代表总水平导数，而VDR则表示垂向导数，u表示位场异常；ASM计算结果的极大值表示地质体的位置，零值刻画的是地质体边界信息。

倾斜角法(TA)，即斜导数法，是一种具有均衡意义的场源边界识别方法，它在1994年由Miller和Singh提出。倾斜角法实质上是各一阶导数组合的比值，具体计算公式如下，

总水平导数的倾斜角(TAHG)法是Ferreira于2013年提出的，在总水平导数的结果上计算倾斜角，以达到边界增强的效果。其计算式子如下

Theta map方法是由Wijns等于2005年首次提出，它实质上是THDR与ASM的比值，表达式如下，

归一化标准差(NSTD)方法属于数理统计类边界识别方法，由Cooper和Cowan于2008年提出，它是垂向坐标方向的一阶导数的标准差与三个坐标方向一阶导数标准差的比值，计算式子如下

其中，u代表位场异常，σ表示标准差。这种方法是利用一定大小的滑动窗口来计算窗口中心点的归一化标准差值，其计算结果的极大值位置代表地质体的边界位置。

上述常规的边界识别方法可以识别位场的线形边界，但对于弧形边界提取效果较差。谱矩方法以随机过程理论为基础，早先被广泛应用于工程学的表面形貌识别，因其能够细致地描述表面的几何特征而被引入地球物理学领域，用于位场表面数据的解释。目前，孙艳云等基于表面二阶谱矩，定义了位场的脊形化系数Λ₂，计算公式如下

式中，Δ₂和M₂均为二阶谱矩的统计不变量，其表达式如下

M₂＝m₂₀+m₀₂ (9)

m₂₀、m₁₁、m₀₂是二阶谱矩元，分别为x和y方向上一阶导数的方差或协方差。由于二阶谱矩表征位场的斜率信息，故该系数主要提取的是地质体的线形边界，而不能有效识别弧形边界。

发明内容

针对传统方法以及表面二阶谱矩识别场源边界的不足，考虑表面四阶谱矩在重磁异常数据处理中的应用。由于四阶谱矩主要描述位场表面与曲率相关的几何特征，因此可利用其提取弧形的边界信息。本发明提出一种基于四阶谱矩的边界识别方法，理论试验结果表明，该方法可以进一步完善弧形边界信息的识别，同时对深度的敏感性较低，可用于提取位场中的弱曲率弧刻痕信号。

本发明所采用的技术方案是：一种基于四阶谱矩的位场边界识别方法，其特征在于，包括以下步骤：

步骤1：在x、y方向上分别以间隔Δx、Δy等间距采样，采样点数分别为M和N，得到研究区域的位场异常数据矩阵A；

步骤2：采用滑动窗口圈定窗口中心点及其邻近数据点的异常值，计算窗口面元的四阶谱矩，将结果定位于面元中心点；

步骤3：逐点移动滑动窗口，得到矩阵A上每个数据点对应的四阶谱矩元m₄₀、m₃₁、m₂₂、m₁₃、m₀₄的值；

步骤4：根据矩阵A中每一个数据点处四阶谱矩元的值计算该点处统计不变量M₄、Δ₄的值：

M₄＝m₄₀+2m₂₂+m₀₄

步骤5：计算矩阵A中每一个数据点处的四阶曲率弧刻痕系数值：

得到矩阵B，即为四阶曲率弧刻痕系数对位场数据的提取结果，其高值指示场源的边界位置。

作为优选，步骤1中，对研究区域的位场异常进行测量，在x、y方向上分别以间隔Δx、Δy等间距采样，采样点数分别为M和N，得到N×M的位场异常数据矩阵A；选取尺寸为N₁×M₁的滑动窗口，沿x方向包含M₁个数据点，沿y方向包含N₁个数据点，M₁、N₁均为奇数。

作为优选，步骤2中，将N×M的矩阵A的列向外延拓(M₁-1)/2列，行向外延拓(N₁-1)/2行，得到大小为[N+(N₁-1)]×[M+(M₁-1)]的新矩阵

对矩阵

中每个数据点求二阶导数，得到

在矩阵

中，选定延拓前第一行第一列的数据点作为滑动窗口的中心点，利用窗口圈定窗口中心点及其邻近数据点的异常值，形成大小为N₁×M₁的局部面元，通过如下公式计算其四阶谱矩元的值：

其中，z_xx(x_i,y_k)、z_xy(x_i,y_k)、z_yy(x_i,y_k)分别为数据点(x_i,y_k)处位场异常的二阶导数；将上述各谱矩元的值赋给窗口中心的数据点。

作为优选，步骤3中，在矩阵

中，按延拓前的数据点逐点移动滑动窗口，计算出每一个数据点处m₄₀、m₃₁、m₂₂、m₁₃、m₀₄的值。

相对于现有技术，本发明的有益效果是：本发明提出的基于四阶谱矩的位场边界识别方法能够均衡、收敛地识别弧形边界，且不产生干扰刻痕。同时，场源埋深的深浅对提取结果影响较小，即该方法对深度的敏感性较低，可用于提取位场中的弱曲率弧刻痕信号。

附图说明

图1为本发明实施例的流程图；

图2为本发明实施例的球体模型平面图；

图3为本发明实施例的不同埋深下球体模型的理论重力异常和四阶曲率弧刻痕系数Λ₄的提取结果：(a)埋深10km的理论重力异常；(b)埋深10km时Λ₄的提取结果；(c)埋深12km的理论重力异常；(d)埋深12km时Λ₄的提取结果；(e)埋深14km的理论重力异常；(f)埋深14km时Λ₄的提取结果；

图4为本发明实施例的常规边界识别方法对埋深10km的球体模型理论重力异常的提取结果：(a)总水平导数(THDR)；(b)总水平导数的倾斜角(TAHG)；(c)Theta map；(d)归一化标准差(NSTD)；

图5为本发明实施例的不同埋深下球体模型的理论磁异常和四阶曲率弧刻痕系数Λ₄的提取结果：(a)埋深10km的理论磁异常；(b)埋深10km时Λ₄的提取结果；(c)埋深12km的理论磁异常；(d)埋深12km时Λ₄的提取结果；(e)埋深14km的理论磁异常；(f)埋深14km时Λ₄的提取结果；

图6为本发明实施例的常规边界识别方法对埋深10km的球体模型理论磁异常的提取结果：(a)总水平导数(THDR)；(b)总水平导数的倾斜角(TAHG)；(c)Theta map；(d)归一化标准差(NSTD)。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

请见图1，本发明提供的一种基于四阶谱矩的位场边界识别方法，包括以下步骤：

步骤1：对研究区域的位场异常进行测量，在x、y方向上分别以间隔Δx、Δy等间距采样，采样点数分别为M和N，得到N×M的位场异常数据矩阵A；选取尺寸为N₁×M₁的滑动窗口提取网格点处的位场异常值，即沿x方向包含M₁个数据点，沿y方向包含N₁个数据点，M₁、N₁均为奇数；

步骤2：将N×M的矩阵A的列向外延拓(M₁-1)/2列，行向外延拓(N₁-1)/2行，得到大小为[N+(N₁-1)]×[M+(M₁-1)]的新矩阵

对矩阵

中每个数据点求二阶导数，得到

在矩阵

中，选定延拓前第一行第一列的数据点作为滑动窗口的中心点，利用窗口圈定中心点及其邻近数据点的异常值，形成大小为N₁×M₁的局部面元。令u(x_i,y_k)为滑动窗口中数据点(x_i,y_k)处的异常值，通过如下公式计算其四阶谱矩元的值：

式中，u_xx(x_i,y_k)、u_xy(x_i,y_k)、u_yy(x_i,y_k)分别为数据点(x_i,y_k)处位场异常的二阶导数；将上述各谱矩元的值赋给窗口中心的数据点；

步骤3：在矩阵

中，按延拓前的数据点逐点移动滑动窗口，计算出每一个数据点处m₄₀、m₃₁、m₂₂、m₁₃、m₀₄的值；

步骤4：根据矩阵A中每一个数据点处四阶谱矩元的值计算该点处统计不变量M₄、Δ₄的值：

M₄＝m₄₀+2m₂₂+m₀₄

步骤5：计算矩阵A中每一个数据点处的四阶曲率弧刻痕系数值：

得到矩阵B，即为四阶曲率弧刻痕系数对位场数据的提取结果，其高值指示场源的边界位置。

本发明提出的基于四阶谱矩的位场边界识别方法，是一种适用于提取弧形边界的方法，克服了常规传统方法和表面二阶谱矩不能有效识别弧形边界的缺陷。重磁异常的仿真试验结果表明，本发明提出的基于四阶谱矩的位场边界识别方法能够均衡收敛地识别弧形边界，且对深度的敏感性较低，可用于提取位场中的弱曲率弧刻痕信号。

为了验证是发明提供的一种基于四阶谱矩的边界识别方法，现将该方法应用于仿真试验数据。由于球体在地面的位场异常具有完整的曲率弧刻痕特征，是应用四阶谱矩来提取位场中曲率弧刻痕信息的理想模型，故设计平面图如图2所示的球体仿真模型，两球体半径均为8km，中心埋深均为10km，球心坐标为(-20km，-20km)的球体剩余密度为0.8g/cm³，磁化强度为0.8A/m，球心坐标为(20km，20km)的球体剩余密度0.2g/cm³，磁化强度为0.3A/m，x轴与y轴方向点距均为0.5km。

在理论试验部分，首先测试本发明提出的四阶曲率弧刻痕系数对深度的敏感性，然后与当前一些常规的位场边界识别方法进行比较，包括总水平导数(THDR)、总水平导数的倾斜角(TAHG)、Thetamap法和归一化标准差(NSTD)，以验证本发明提出的边界识别方法的有效性。以下给出两组试验分别检验四阶曲率弧刻痕系数在重、磁异常数据中的边界识别效果。

重、磁异常仿真试验中选取的滑动窗口大小分别为5×5和3×3。图3和图5分别为该球体模型在中心埋深10km、12km、14km时的理论重、磁异常及对应四阶曲率弧刻痕系数Λ₄的提取结果。分析其提取结果可知，Λ₄提取的曲率弧分布信息受场源强度影响较小，能稳定、有效地提取不同深度地质体的弧形边界信息，对深度的敏感性较低。图4和图6分别为几种常规边界识别方法对重、磁异常的提取结果，可见四种传统方法都能提取出两个球状场源的边界，但均存在边界发散的问题，除总水平导数外的另外三种方法的提取结果中存在干扰刻痕，这些缺点使得实际应用中难以准确分辨出地质体的边界位置。对比图3(b)和图5(b)易知，四阶曲率弧刻痕系数Λ₄相较于传统方法能更准确、收敛地刻画球状重磁异常体的边界，且没有干扰刻痕。从实验结果可以看出，本实施例的边界识别方法在提取场源边界时明显优于常规的传统方法。

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对较佳实施例的描述较为详细，并不能因此而认为是对本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

Claims (4)

1.一种基于四阶谱矩的位场边界识别方法，其特征在于，包括以下步骤：

步骤1：在x、y方向上分别以间隔Δx、Δy等间距采样，采样点数分别为M和N，得到研究区域的位场异常数据矩阵A；

步骤2：采用滑动窗口圈定窗口中心点及其邻近数据点的异常值，计算窗口面元的四阶谱矩，将结果定位于面元中心点；

步骤3：逐点移动滑动窗口，得到矩阵A上每个数据点对应的四阶谱矩元m₄₀、m₃₁、m₂₂、m₁₃、m₀₄的值；

步骤4：根据矩阵A中每一个数据点处四阶谱矩元的值计算该点处统计不变量M₄、Δ₄的值：

M₄＝m₄₀+2m₂₂+m₀₄

步骤5：计算矩阵A中每一个数据点处的四阶曲率弧刻痕系数值：

得到矩阵B，即为四阶曲率弧刻痕系数对位场数据的提取结果，其高值指示场源的边界位置。

2.根据权利要求1所述的基于四阶谱矩的位场边界识别方法，其特征在于：步骤1中，对研究区域的位场异常进行测量，在x、y方向上分别以间隔Δx、Δy等间距采样，采样点数分别为M和N，得到N×M的位场异常数据矩阵A；选取尺寸为N₁×M₁的滑动窗口，沿x方向包含M₁个数据点，沿y方向包含N₁个数据点，M₁、N₁均为奇数。

3.根据权利要求1所述的基于四阶谱矩的位场边界识别方法，其特征在于：步骤2中，将N×M的矩阵A的列向外延拓(M₁-1)/2列，行向外延拓(N₁-1)/2行，得到大小为[N+(N₁-1)]×[M+(M₁-1)]的新矩阵

对矩阵

中每个数据点求二阶导数，得到

在矩阵

中，选定延拓前第一行第一列的数据点作为滑动窗口的中心点，利用窗口圈定窗口中心点及其邻近数据点的异常值，形成大小为N₁×M₁的局部面元，通过如下公式计算其四阶谱矩元的值：

其中，z_xx(x_i,y_k)、z_xy(x_i,y_k)、z_yy(x_i,y_k)分别为数据点(x_i,y_k)处位场异常的二阶导数；将上述各谱矩元的值赋给窗口中心的数据点。

4.根据权利要求2所述的基于四阶谱矩的位场边界识别方法，其特征在于：步骤3中，在矩阵

中，按延拓前的数据点逐点移动滑动窗口，计算出每一个数据点处m₄₀、m₃₁、m₂₂、m₁₃、m₀₄的值。

Images