CN111580174A - 一种基于帕德近似的重磁数据向下延拓方法 - Google Patents

Abstract

本发明提供了一种基于帕德近似的重磁数据向下延拓方法。该方法包括：利用ISVD方法计算重磁数据的各阶垂向导数；利用各阶垂向导数计算帕德近似的系数方程，并求解系数方程的系数；利用求得的系数代入帕德近似，计算重磁数据的向下延拓结果。本发明提供的这种基于帕德近似的重磁数据向下延拓方法能够提供稳定且准确的重磁数据向下延拓。

Description

一种基于帕德近似的重磁数据向下延拓方法

技术领域

本发明涉及勘查技术与工程技术领域，特别是涉及一种基于帕德近似的重磁数据向下延拓方法。

背景技术

向下延拓是卫星、航空重磁数据处理中的经典和关键技术。重磁数据向下延拓可以增强数据信息，识别弱异常和叠加异常，这对近地表地球物理勘探中从强干扰中提取弱信息非常重要。泰勒级数展开法是一种简单快速的重磁数据向下延拓方法。泰勒级数展开法可以一定程度改善向下延拓的不稳定性问题，但是向下延拓的准确性也有待提高。

重磁场向下延拓的泰勒级数展开式可以表示为(Evjen，1936)

其中，u(x,y,z₀)为z₀高度的重磁场，h(h＞0)是向下延拓的高度，利用公式 (1)进行向下延拓的关键在于计算重磁场的不同阶垂向导数。Evjen(1936) 利用拉普拉斯方程计算计算高阶垂向导数，拉普拉斯方程如下：

计算重磁数据高阶垂向导数方法

但他利用的一阶垂向导数是用图解法计算出来的，其计算不便于现代形式的数字数据且所计算结果精度低。

自快速傅里叶变换(FFT)算法定义以来(Dean，1958；Cooley&Tukey， 1965)，Gunn(1975)实现了使用傅立叶变换来执行垂直导数。Blakely (1996)给出了波数域中重磁数据的垂直导数计算的滤波器响应的二维表达式，

其中，F和F^-1是傅里叶正反变换对，n是垂向导数的阶数。然而，由于波数域中高频部分的急剧放大作用，FFT方法中所使用的数据需要很好的去燥和平滑滤波，且由FFT方法计算的垂向导数也不稳定，这限制了重磁数据向下延拓的深度不超过2个数据点间隔。

Fedi和Florio(2001、2002)提出了一种计算垂直导数的积分二次垂向导数方法(ISVD方法)，该方法结果比FFT方法更稳定。为了计算重磁数据的一阶垂直导数，他们首先用波数域对重磁数据进行积分，得到重磁场的标量位

其次，根据拉普拉斯方程(3)，他们通过积分函数的二阶水平导数之和(通过有限差分计算)计算积分函数的二阶垂直导数。至于二阶水平导数，它们的计算可以在空间域中使用有限差分关系(例如Cordell和Grauch，1985； Trompat和等人，2003)以精确和稳定的方式完成，有限差分表达式

最后，他们重复第二步计算重磁数据的高阶垂向导数。利用ISVD法计算重磁数据的垂向导数时，采用了平滑滤波(垂直积分滤波)和有限差分法相结合的方法，比FFT法稳定得多。这种方法的优点是使重磁数据向下延拓的稳定性大大提高，但是向下延拓结果的精度仍然有待于提高。

基于泰勒级数展开，利用ISVD计算垂向导数的向下延拓称为 ISVD-Taylor向下延拓方法，上述方法提高了向下延拓的稳定性。然而，泰勒级数展开式不是数学逼近理论中最精确的近似(Ma et al，2002)。帕德近似(本文翻译自Padé approximation)对函数的逼近比Taylor级数展开的精度高，误差小(Zhang等人，2016)。

发明内容

本发明要解决的技术问题是提供一种基于帕德近似的重磁数据向下延拓方法及装置，能够提供稳定且准确的重磁数据向下延拓。

为解决上述技术问题，本发明提供了一种基于帕德近似的重磁数据向下延拓方法，所述方法包括：利用ISVD方法计算重磁数据的各阶垂向导数；利用各阶垂向导数计算帕德近似的系数方程，并求解系数方程的系数；利用求得的系数代入帕德近似，计算重磁数据的向下延拓结果。

在一些实施方式中，利用ISVD方法计算重磁数据的各阶垂向导数，包括：利用波数域的重磁数据换算方法计算重磁数据的标量位；利用有限差分法，根据标量位计算标量位水平方向的二阶导数；利用拉普拉斯方程，根据标量位水平方向的二阶导数得到标量位的垂直方向的二阶导数；重复上述计算过程，得到重磁数据的各阶垂向导数。

在一些实施方式中，利用波数域的重磁数据换算方法计算重磁数据的标量位，包括：根据如下公式计算重磁数据的标量位：

其中，F表示傅里叶变换，F^-1表示傅里叶逆变换，k_x,k_y表示与空间域坐标x,y对应的波数域坐标，z₀表示重磁数据在垂向坐标z的高度位置， u(x,y,z₀)是观测高度为z₀的重磁数据，

表示波数域垂向积分算子，v(x,y,z₀)表示换算得到的重磁数据的标量位。

在一些实施方式中，利用有限差分法，根据标量位计算标量位水平方向的二阶导数，包括：根据如下公式计算标量位水平方向的二阶导数：

其中，Δx与Δy均表示水平方向的采样间距，v_xx表示x轴方向的二阶导数，v_yy表示y轴方向的二阶导数。

在一些实施方式中，利用拉普拉斯方程，根据标量位水平方向的二阶导数得到标量位的垂直方向的二阶导数，包括：根据如下公式计算标量位的垂直方向的二阶导数：

u_z(x,y,z₀)＝v_zz(x,y,z₀)＝-[v_xx(x,y,z₀)+v_yy(x,y,z₀)]

其中，u_z表示重磁数据垂直方向的一阶导数，v_zz表示标量位垂直方向的二阶导数，v_xx及v_yy表示标量位的水平方向的二阶导数。

在一些实施方式中，利用各阶垂向导数计算帕德近似的系数方程，并求解系数方程的系数，包括：利用二阶垂向导数计算帕德近似的系数方程，并求解系数方程的系数。

在一些实施方式中，利用二阶垂向导数计算帕德近似的系数方程，并求解系数方程的系数，包括：根据如下公式计算系数方程的系数：

其中，p₀、q₁、q₂为帕德近似的系数方程中的系数，c₀、c₁、c₂为各阶垂向导数。

在一些实施方式中，利用求得的系数代入帕德近似，计算重磁数据的向下延拓结果，包括：根据如下公式计算重磁数据的向下延拓结果：

其中，h表示向下延拓的高度，u(x,y,z₀+h)表示重磁数据的向下延拓结果，u(x,y,z₀)表示重磁数据。

采用这样的设计后，本发明至少具有以下优点：

本发明提供的一种基于积分垂向导数的帕德近似重磁数据向下延拓方法，与传统方法相比，本发明利用帕德近似直接对重磁数据进行函数逼近，而不是间接逼近，计算简单、直观。本发明的帕德近似是2阶的有理逼近，在逼近理论上其精度比泰勒级数展开高，而且在泰勒级数展开不收敛的情况也可以保证收敛，提高了向下延拓结果的准确性。采用的积分二次垂向导数方法计算导数，有效增强了换算的稳定性和准确性。本发明中的帕德近似具有不同阶次的表达公式，本发明为了实现直观，选择计算误差更小、更稳定的二阶帕德近似，实现重磁数据的稳定且准确的向下延拓。

附图说明

上述仅是本发明技术方案的概述，为了能够更清楚了解本发明的技术手段，以下结合附图与具体实施方式对本发明作进一步的详细说明。

图1是本发明实施例提供的基于帕德近似的重磁数据向下延拓方法的流程图；

图2a是在z＝0米处场源体的无噪重力数据的示意图；

图2b是在z＝2米处、向下延拓前的无噪重力数据的示意图；

图2c是由z＝2米至z＝0米处的由FFT方法的向下延拓的示意图；

图2d是由z＝2米至z＝0米处的由FFT泰勒方法的向下延拓的示意图；

图2e是由z＝2米至z＝0米处的由ISVD泰勒方法的向下延拓的示意图；

图2f是由z＝2米至z＝0米处的由ISVD帕德方法的向下延拓的示意图；

图3a是在z＝0米处场源体的无噪重力数据的示意图；

图3b是在z＝5米处时向下延拓前的无噪重力数据的示意图；

图3c是由z＝5米至z＝0米处的由FFT方法的向下延拓的示意图；

图3d是由z＝5米至z＝0米处的由FFT泰勒方法的向下延拓的示意图；

图3e是由z＝5米至z＝0米处的由ISVD泰勒方法的向下延拓的示意图；

图3f是由z＝5米至z＝0米处的由ISVD帕德方法的向下延拓的示意图；

图4a是在z＝0米处场源体的无噪重力数据的示意图；

图4b是在z＝5米处时向下延拓前的无噪重力数据的示意图；

图4c是z＝5米至z＝0米处的由FFT方法的向下延拓的示意图；

图4d是z＝5米至z＝0米处的由FFT泰勒方法的向下延拓的示意图；

图4e是z＝5米至z＝0米处的由ISVD泰勒方法的向下延拓的示意图；

图4f是z＝5米至z＝0米处的由ISVD帕德方法的向下延拓的示意图；

图5a是在z＝0米处场源体的无噪重力数据的示意图；

图5b是在z＝15米向下延拓前的无噪重力数据的示意图；

图5c是z＝15米至z＝0米处的由FFT方法的向下延拓的示意图；

图5d是z＝15米至z＝0米处的由FFT泰勒方法的向下延拓的示意图；

图5e是z＝15米至z＝0米处的由ISVD泰勒方法的向下延拓的示意图；

图5f是z＝15米至z＝0米处的由ISVD帕德方法的向下延拓的示意图；

图6是不同方法自身的向下延拓值的均方根值；

图7a是在z＝0米处场源体的无噪重力数据的示意图；

图7b是在z＝9米处向下延拓前的无噪重力数据的示意图；

图7c是z＝9米至z＝0米处的由FFT方法的向下延拓的示意图；

图7d是z＝9米至z＝0米处的由FFT泰勒方法的向下延拓的示意图；

图7e是z＝9米至z＝0米处的由ISVD泰勒方法的向下延拓的示意图；

图7f是z＝9米至z＝0米处的由ISVD帕德方法的向下延拓的示意图；

图8a是Nechako盆地的布格重力异常；

图8b是通过FFT泰勒方法至z＝0.4千米处的向下延拓值；

图8c是通过ISVD泰勒方法至z＝0.4千米处的向下延拓值；

图8d是通过ISVD帕德方法至z＝0.4千米处的向下延拓值；

图9a是Nechako盆地的布格重力异常；

图9b是通过FFT泰勒方法至z＝2.0千米处的向下延拓值；

图9c是通过ISVD泰勒方法至z＝2.0千米处的向下延拓值；

图9d是通过ISVD帕德方法至z＝2.0千米处的向下延拓值。

具体实施方式

以下结合附图对本发明的优选实施例进行说明，应当理解，此处所描述的优选实施例仅用于说明和解释本发明，并不用于限定本发明。

本发明提出一种基于帕德近似的重磁数据向下延拓方法。图1示出了该方法的流程图。参见图1，技术方案的具体步骤如下：

S11，首先利用积分二次垂向导数方法(ISVD方法)计算重磁数据的各阶垂向导数{c_n}。

S12，然后利用步骤1获得各阶垂向导数{c_n}计算帕德近似的系数方程 P/Q，并求解其中的系数p_L,q_M。

S13，最后利用把步骤2求得的系数p_L,q_M代入帕德近似，计算重磁数据的向下延拓结果u(x,y,z₀+h)。

本发明属于勘查技术与工程领域的地球探测与信息技术，本技术实现了卫星、航空重磁数据的向下延拓，在保证稳定性的前提下，提高了向下延拓结果的精确性，对近地表重磁场勘探的精度改善具有十分重要意义。

上述步骤1中{c_n}是z₀(也是观测面的高度)高度上的重磁数据(实际观测的数据)。其一般表达式为：

本发明利用的是二阶帕德近似，因此所使用的{c_n}表达式为

公式(7)和(8)中的各阶垂向导数的计算是通过积分二次垂向导数方法(ISVD方法)的计算获得。

上述步骤1中积分二次垂向导数方法(ISVD方法)的计算过程为：

首先利用波数域的重磁数据换算方法计算u(x,y,z0+h)的标量位 v(x,y,z₀+h)：

其中，F表示傅里叶变换，F^-1表示傅里叶逆变换，k_x,k_y表示与空间域坐标x,y对应的波数域坐标，u(k_x,k_y,z₀)是u(x,y,z₀)的波数域表达，

表示波数域垂向积分算子。

上述涉及的波数域是通过傅里叶变换得到的，其中傅里叶变换及逆变换为：

然后利用有限差分法计算标量位

水平方向的二阶导数：

其中Δx,Δy表示水平方向的采样间距。

最后利用拉普拉斯方程，得到标量位v(x,y,z₀+h)的垂直方向的二阶导数，即重磁数据

的垂直方向的一阶导数u_z(x,y,z₀+h)：

u_z(x,y,z₀+h)＝v_zz(x,y,z₀+h)＝-[v_xx(x,y,z₀+h)+v_yy(x,y,z₀+h)](12)

在上式中，标量位的垂直方向的二阶导数，即是重磁数据的一阶垂向导数。再利用傅里叶逆变换得到：

u_z(x,y,z₀+h)＝F^-1[u_z(k_x,k_y,z₀+h)] (13)

上述步骤2中帕德近似的系数方程P/Q的一般形式为

本发明利用的是二阶帕德近似，因此所使用的系数方程P/Q实际表达式为：

代入{c_n}表达式(8)：

后得到：

上述步骤2中二阶帕德近似的系数方程P/Q的系数p_L,q_M求解为：

上述步骤3中代入系数p_L,q_M到帕德近似，计算重磁数据的向下延拓结果u(x,y,z₀+h)的一般形式为：

实际的二阶帕德近似重磁数据向下延拓表达式为：

在数学理论中泰勒级数展开并不是最为准确的近似。作为近似函数，帕德近似具有相较于泰勒级数展开更为优良的准确性和更小的误差。Zhang 等人使用帕德近似而不是泰勒级数展开来进行正演模拟及利用位场数据的地下界面反演。Zhou等人提出了一种使用帕德近似的向下延拓方法来在频域进行向下延拓因子的拓展，其中，参数k是波数模式。并且，为了避免高阶效应，Zhou等人仅仅示出了公式近似中3/2的项。

与Zhou等人提出的频域中的公式不同，我们提出使用帕德近似直接对位场进行拓展。

我们通过求解如下的线性方程来获得系数P_L和Q_M：

在本文中，我们称公式(18)为ISVD帕德方法，其中的系数通过方程 (14)及(20)计算。

由于帕德近似是多项式总和的有理函数，而多项式的可能阶数无限并且形式可变，我们将向下延拓的最高阶数限制为2，并将形式限制为0/2。

这样，对应的系数方程(14)成为：

将c_n带入(15)，我们得到：

通过求解方程(16)，我们得到基于二阶帕德近似的向下延拓：

3.综合结果

在本部分，我们为了测试本文提出方法(ISVD帕德方法)的性能而将其应用至综合重力数据。为了验证向下延拓方法的成就，第一模型是一个中心位于(55,45,13)米，范围是(40,10,20)米的长方体。图2a示出了位于地面0米处的长方体的理论无噪重力数据，其中地面位于距长方体顶面3米的位置。并且，数据的网格间距是1米，并且以后的数据也是这样。图2b示出了位于地面上2米处的长方体的理论无噪重力数据，其中地面位于由测量面距长方体顶面5米的位置。在图2c至图2f中的向下延拓的结果分别通过在图2b示出的理论无噪重力数据应用FFT方法、FFT泰勒方法、 ISVD泰勒方法及ISVD帕德方法至地面而运算得到。我们可以看到，当向下延拓的深度为2米(点间隔为2)时，全部四种方法能够完成向下延拓。采用FFT泰勒方法及ISVD泰勒方法的向下延拓的结果具有相较于FFT方法更少的边界失真。然而，本文提出的方法(ISVD帕德方法)是四种方法中最为准确的。

我们将向下延拓的深度增加至5米(点间隔为5)，结果被示出于图3 中。图4中的结果与图2中类似，除了FFT方法失效以外。

为了估计提出方法在向下延拓中的能力，应用三个长方体的第二模型被提出。在三个长方体中，第一个的中心位于(55,45,13)米，范围是 (40,10,20)米，第二个的中心位于(65,75,12)米，范围是(12,15,20) 米，第三个的中心位于(82,75,12)米，范围是(12,10,22)米。图4a示出了位于0米处的第二模型的理论无噪重力数据，其数据的点间隔与图2 中示出的相同。图4b示出了位于0米以上5米处的第二模型的理论无噪重力数据，其点间隔为5。通过FFT方法、FFT泰勒方法、ISVD泰勒方法及 ISVD帕德方法计算得到的向下延拓的结果在图2c至图2f中对应示出。当向下延拓的深度为5米(点间隔为5)并且模型复杂时，我们可以看到通过 ISVD帕德方法的向下延拓在全部方法中最为准确和稳定。它几乎覆盖了图 4a中的全部原始的异常，这包括一些小规模异常以及一些细节信息。

为了比较不同的方法，我们通过在1米至15米之间改变深度来完成向下延拓。15米深度处的结果被展示在图5中，这些方法的向下延拓结果相互之间有很大差别。我们可以看到即便向下延拓的深度是15米(点间隔为 15)，尽管对于向下延拓来说这是很大的深度，通过本文提出的方法计算得到的结果相较于理论值是稳定和准确的。由ISVD泰勒方法得到的结果也是稳定的，但是与本文提出方法相比则不是很准确。通过FFT方法和FFT 泰勒方法的向下延拓不能产生最终的结果。我们计算通过不同方法得到的向下延拓值自身的均方差(RMS)，并同时计算在0米深度处的综合数据与向下延拓值之间的均方误差(RMS误差)。所有的结果在图6中被示出。图 6a中的结果示出了FFT泰勒方法计算得到的向下延拓结果在7米范围内小于理论值，而在7米范围外大于理论值。通过ISVD泰勒方法计算得到的向下延拓结果始终小于理论值。通过ISVD帕德方法计算得到的向下延拓结果与其他两种方法相比较与理论值最为接近。在图6b中，我们可以看到ISVD 泰勒方法得到的理论值与向下延拓结果之间的RMS误差小于FFT方法，并且通过ISVD帕德方法得到的理论值与向下延拓结果之间的RMS误差为最小。

考虑到数据中的噪声，我们为数据中添加2.5％的高斯噪声，并将数据向下延拓至5米(点间隔为5)。并且，在进行向下延拓之前，使用2D中值滤波器及维纳滤波器来抑制噪声。可以看出，本文提出方法的结果相比于其他二者更为准确。但是，由于对于噪声污染的高密感性，本文提出方法并不稳定。

4.实际应用

在本部分中，一个实际的例子被执行，以便验证本文提出的方法。这个例子的空中重力数据在加拿大的Nechako盆地上测量，如Zhang等人在 2018年的应用。数据的点间隔为0.4千米×0.4千米，网格均布点的数量为179×179。向下延拓的深度是0.4千米(点间隔为1)，及2.0千米(点间隔为5)。向下延拓的全部其他参数均在TaPaDown.m的MATLAB代码中设置。图8a至图8d分别示出了重力异常、通过FFT泰勒方法的向下延拓的结果、通过ISVD泰勒方法的向下延拓的结果，以及ISVD帕德方法的向下延拓方法。我们可以看到，这三种方法可以得到0.4千米深度的向下延拓 (点间隔为1)。由FFT泰勒方法计算得到的结果示出了较少的异常小特征，而且甚至扭曲了异常的真实特征。通过ISVD泰勒方法得到的向下延拓相比于FFT泰勒方法更为稳定。而且，一些细节信息丢失了。本文提出的方法 (ISVD帕德方法)能够成功的向下延拓重力数据。本文提出方法的结果收敛且准确。并且它们比FFT泰勒方法及ISVD泰勒方法显示了更为细节的信息。我们将向下延拓的深度深化至2.0千米(点间隔为5)，结果显示通过FFT 泰勒方法及ISVD泰勒方法的向下延拓不再可行，然而通过本文提出方法的向下延拓仍然成功和稳定。我们能够看到ISVD帕德方法对于实际数据的向下延拓有效。

在本文中，我们提出了基于帕德近似展开的向下延拓方法，其导数通过在波数域及空域中的ISVD方法计算，这样能够改善向下延拓的稳定性和准确性，被称为ISVD帕德方法。同时，通过FFT泰勒、ISVD泰勒及ISVD 帕德方法，对向下延拓进行对比。我们使用综合模型来验证提出的方法。结果显示，本文提出方法，也就是ISVD帕德方法相比于FFT泰勒及ISVD 泰勒，能够得到更为稳定和准确的向下延拓。我们同时使用实际数据来验证提出的方法。结果显示，通过ISVD帕德方法的向下延拓能够有效的复现原始的小异常以及发现局部信息。

以上所述，仅是本发明的较佳实施例而已，并非对本发明作任何形式上的限制，本领域技术人员利用上述揭示的技术内容做出些许简单修改、等同变化或修饰，均落在本发明的保护范围内。

Claims (8)

1.一种基于帕德近似的重磁数据向下延拓方法，其特征在于，包括：

利用ISVD方法计算重磁数据的各阶垂向导数；

利用各阶垂向导数计算帕德近似的系数方程，并求解系数方程的系数；

利用求得的系数代入帕德近似，计算重磁数据的向下延拓结果。

2.根据权利要求1所述的基于帕德近似的重磁数据向下延拓方法，其特征在于，利用ISVD方法计算重磁数据的各阶垂向导数，包括：

利用波数域的重磁数据换算方法计算重磁数据的标量位；

利用有限差分法，根据标量位计算标量位水平方向的二阶导数；

利用拉普拉斯方程，根据标量位水平方向的二阶导数得到标量位的垂直方向的二阶导数；

重复上述计算过程，得到重磁数据的各阶垂向导数。

3.根据权利要求2所述的基于帕德近似的重磁数据向下延拓方法，其特征在于，利用波数域的重磁数据换算方法计算重磁数据的标量位，包括：

根据如下公式计算重磁数据的标量位：

其中，F表示傅里叶变换，F^-1表示傅里叶逆变换，k_x,k_y表示与空间域坐标x,y对应的波数域坐标，z₀表示重磁数据在垂向坐标z的高度位置，u(x,y,z₀)是观测高度为z₀的重磁数据，

表示波数域垂向积分算子，v(x,y,z₀)表示换算得到的重磁数据的标量位。

4.根据权利要求3所述的基于帕德近似的重磁数据向下延拓方法，其特征在于，利用有限差分法，根据标量位计算标量位水平方向的二阶导数，包括：

根据如下公式计算标量位水平方向的二阶导数：

其中，Δx与Δy均表示水平方向的采样间距，v_xx表示x轴方向的二阶导数，v_yy表示y轴方向的二阶导数。

5.根据权利要求4所述的基于帕德近似的重磁数据向下延拓方法，其特征在于，利用拉普拉斯方程，根据标量位水平方向的二阶导数得到标量位的垂直方向的二阶导数，包括：

根据如下公式计算标量位的垂直方向的二阶导数：

u_z(x,y,z₀)＝v_zz(x,y,z₀)＝-[v_xx(x,y,z₀)+v_yy(x,y,z₀)]

其中，u_z表示重磁数据垂直方向的一阶导数，v_zz表示标量位垂直方向的二阶导数，v_xx及v_yy表示标量位的水平方向的二阶导数。

6.根据权利要求1或2所述的基于帕德近似的重磁数据向下延拓方法，其特征在于，利用各阶垂向导数计算帕德近似的系数方程，并求解系数方程的系数，包括：

利用二阶垂向导数计算帕德近似的系数方程，并求解系数方程的系数。

7.根据权利要求6所述的基于帕德近似的重磁数据向下延拓方法，其特征在于，利用二阶垂向导数计算帕德近似的系数方程，并求解系数方程的系数，包括：

根据如下公式计算系数方程的系数：

其中，p₀、q₁、q₂为帕德近似的系数方程中的系数，c₀、c₁、c₂为各阶垂向导数。

8.根据权利要求1所述的基于帕德近似的重磁数据向下延拓方法，其特征在于，利用求得的系数代入帕德近似，计算重磁数据的向下延拓结果，包括：

根据如下公式计算重磁数据的向下延拓结果：

其中，h表示向下延拓的高度，u(x,y,z₀+h)表示重磁数据的向下延拓结果，u(x,y,z₀)表示重磁数据。

Images