CN107291659B - 平面地磁异常场一步向上延拓平面模量梯度场的递归余弦变换法 - Google Patents

Abstract

本发明提供的是一种平面地磁异常场一步向上延拓平面模量梯度场的递归余弦变换法。先进行平面地磁异常场测量数据预处理；利用径向基函数插值法将观测平面上不规则离散分布的地磁异常场测量数据进行网格化，获得网格化的地磁异常测量数据。再利用递归余弦变换将网格化的地磁异常测量数据转换为余弦变换波谱；利用一步波数域转换因子，由地磁异常余弦变换波谱计算地磁模量梯度场各分量的波谱。最后分别对水平x分量、水平y分量、垂直z分量的波谱进行行向移相π/2的递归反余弦变换得到延拓平面上的地磁模量梯度场的分量。本发明的一步法的计算量约为两步法计算量的一半；余弦变换法减小Gibbs边界效应，提高了算法精度；可应用于任意长度的数据序列。

Description

平面地磁异常场一步向上延拓平面模量梯度场的递归余弦变 换法

技术领域

本发明涉及的是一种地磁位场数据处理方法，具体涉及到一种平面地磁异常场一步向上延拓平面模量梯度场的递归余弦变换法。

背景技术

地磁延拓是一种重要的地磁位场数据处理方法，在地磁匹配定位、地磁观测及其资料解释等方面具有重要应用。地磁基准图是实现地磁匹配定位的数据基础。为进行地磁匹配定位，一般需要获得航行深度面上的地磁基准图。不同航行器的航行高度各不相同，通过直接测量方式得到各航行高度上的地磁数据是不现实的。大面积的地磁测量值多数为地磁异常场数据，地磁模量梯度测量虽已开始，但匹配定位需要大范围导航区域的基准图，现有的地磁模量梯度测量数据难以满足这一要求。根据已有的航测地磁数据或者海测地磁数据，通过延拓和导数的同时计算一步获得地磁模量梯度场，是获得航行面地磁模量梯度场的有效途径之一。向上延拓还可融合不同高度的航磁数据，使之归化到同一高度面。地磁观测及其资料解释方面，利用向上延拓可以进行区域场与局部场、深源场与浅源场的分离。

地磁模量梯度场具有以下优点：第一，地磁模量梯度与不同磁传感器间的姿态无关，降低了多个传感器测量轴配置精度的要求；第二，便于根据载体磁场分布特性，选择载体磁场模量相等的地点合理布置磁传感器，通过求取模量梯度消除载体磁场等的影响。第三，地磁模量梯度场有三个分量，可利用的信息源较地磁异常场多；第四，模量梯度场能有效地突出浅源异常。

Cordell L和Grauch V指出因有限截断和离散采样，非周期函数的离散傅里叶变换会导致频谱的“零频”过小而高频过高；为此提出先将观测数据拓展成周期函数，再进行离散傅里叶变换。在拓展之前采用等效源方法对观测数据进行扩边(Cordell L,GrauchV.Reconciliation of the discrete and integral Fourier transforms,Geophysics,1982,47(2):337-243)。Ricard Y和Blakely R J对边界效应进行了较深入分析，提出了一种经验方法提高频谱的计算精度(A method to minimize edge effects in two-dimensional discrete Fourier transforms,Geophysics,1988,53(8):1113-1117)。郭华等人利用理论公式推导证明实测梯度数据进行延拓转换处理的可行性，研究了地磁梯度延拓后的数据解释规律，阐明了梯度数据向上延拓的理论意义(郭华，王平，朱春华，杜颖.向上延拓对航磁梯度数据的影响及其规律研究，地球物理学进展，2015，30(3):1214-1223)。王小兵等人利用基于傅里叶变换的频率域迭代法，测试了实际基准图的延拓精度和对误差的方法效果(王小兵，张金生，王哲，王仕成.地磁匹配中位场频率域延拓方法研究，探测与控制学报，2009,31(增刊)：29-33)。

目前的航空或海洋测量数据主要是地磁异常值，其向上延拓方法也多针对地磁异常场，如何由平面地磁异常场一步向上延拓获得平面地磁场模量梯度场的问题还较少涉及。基于傅里叶变换的平面地磁异常向上延拓平面地磁模量梯度场存在两种两步法，一种是先由观测平面上的地磁异常场转换模量梯度场，再由模量梯度场向上延拓得到延拓平面上的模量梯度场；另一种先由观测平面上的地磁异常场向上延拓得到延拓面上的地磁异常场，再由延拓面上的地磁异常场转换得到模量梯度场。但这两种两步法都会使用两步基于傅里叶变换的波谱变换与反变换，增加了这种算法的计算时间，不利于向上延拓与位场转换的实时性。另外，傅里叶变换通常要求数据序列的空间周期性，而实际上数据序列两个边界值往往不相等。因此，使用傅里叶变换进行位场延拓或转换都会在边界形成一个跳跃，增加了高频分量，形成Gibbs边界效应。余弦变换可看作实偶序列的傅里叶变换，变换的输入和输出都为实数。与傅里叶变换相比，余弦变换具有更高的能量压缩性能，在一阶马尔科夫过程中依据最小均方误差原则是最接近Karhunen-Loeve变换性能的，可减小Gibbs边界效应。

发明内容

本发明的目的在于提供一种能简化位场延拓与转换的步骤，提高一步位场延拓与转换的精度的平面地磁异常场一步向上延拓平面模量梯度场的递归余弦变换法。

本发明的目的是这样实现的：

步骤1、选取网格化参数，对平面地磁异常场测量数据进行预处理，剔除超出合理范围的粗大测量误差；

步骤2、利用径向基函数插值法，将观测平面上的离散地磁异常场测量数据进行网格化，获得观测平面上的网格化地磁异常测量数据ΔT(n_x,n_y,z₀)，n_x＝1,2,…,M_x，n_y＝1,2，…,M_y，观测平面方程为z＝z₀，向上延拓平面方程为z＝z₀-Δz，z轴垂直向下；

步骤3、先后对观测平面上的地磁异常场进行行向和列向的一维递归余弦变换，这两个方向的计算次序可以交换，得到地磁异常场的二维余弦变换波谱ΔT^C(u_x,u_y,z₀)；

步骤4、按式

由观测平面上的地磁异常场波谱ΔT^C(u_x,u_y,z₀)计算延拓平面上的地磁模量梯度场的波谱

和

式中，R_1WC(u_x,u_y,Δz)为余弦变换下的由观测平面地磁异常场计算延拓平面地磁模量梯度场的一步波数域转换因子，其表达式为

步骤5、先后对

进行行向的移相π/2的一维递归反余弦变换和列向的一维递归反余弦变换，定义为行向正弦、列向余弦反变换，这两个方向的计算次序可以交换，得到延拓平面上的地磁模量梯度场的水平x分量ΔT_x(n_x,n_y,z₀-Δz)；

步骤6、先后对

进行行向的一维递归反余弦变换和列向的移相π/2的一维递归反正弦变换，定义为行向余弦、列向正弦反变换，这两个方向的计算次序可以交换，得到延拓平面上的地磁模量梯度场ΔT_y(n_x,n_y,z₀-Δz)的水平y分量；

步骤7、先后对

进行行向和列向的一维递归反余弦变换，这两个方向的计算次序可以交换，得到延拓平面上的地磁模量梯度场ΔT_z(n_x,n_y,z₀-Δz)的垂直z分量。

本发明的一种平面地磁异常场一步向上延拓平面模量梯度场的递归余弦变换法。先进行平面地磁异常场测量数据预处理，剔除粗大误差；利用径向基函数插值法将观测平面上不规则离散分布的地磁异常场测量数据进行网格化，获得网格化的地磁异常测量数据。再利用递归余弦变换将网格化的地磁异常测量数据转换为余弦变换波谱；利用一步波数域转换因子，由地磁异常余弦变换波谱计算地磁模量梯度场各分量的波谱。最后对水平x分量的波谱进行行向移相π/2的递归反余弦变换，得到延拓平面上的地磁模量梯度场的水平x分量；对水平y分量的波谱进行列向移相π/2的递归反余弦变换，得到延拓平面上的地磁模量梯度场的水平y分量；对垂直z分量的波谱进行递归反余弦变换，得到延拓平面上的地磁模量梯度场的垂直z分量。其中对于并行计算而言，步骤5、6和7是可以进行的。

本发明提出了一种基于递归余弦变换与反变换(含移相π/2的递归余弦反变换)的平面地磁异常场一步向上延拓平面地磁模量梯度场的方法，推导出余弦变换下的由观测平面地磁异常场计算延拓平面地磁模量梯度场的一步波数域转换因子。利用递归余弦变换与反变换(含移相π/2的递归余弦反变换)，由观测平面上的地磁异常场一步计算延拓平面上的地磁模量梯度场。由于每一步位场转换都有计算时间的消耗和计算误差的存在，因此相比于两步法，一步法的计算时间和计算误差都得以减少，也简化了计算流程。递归余弦变换与反变换可以适应任意长度的数据序列，快速傅里叶变换与反变换通常要求数据序列的长度是基数的整数次幂，补零有时会增大不少的计算量。

本发明的方法能由平面地磁异常场直接地一步向上延拓，同时得到延拓平面上的地磁模量梯度场的三个分量；解决现有技术由观测平面上的地磁异常场获得向上延拓面上的地磁模量梯度场，需进行向上延拓和地磁异常场转换地磁模量梯度场的两步转换问题；两步位场转换通常需要两步波谱变换和反变换，本方法仅需一步波谱变换和反变换，简化了位场延拓与转换的步骤，减少了计算时间；而且本方法采用递归余弦变换与反变换，在保证计算量小的同时，相比于傅里叶变换与反变换减小Gibbs边界效应，提高了这种一步位场延拓与转换的算法精度。

本发明与现有技术比较具有以下优点：提出的一种平面地磁异常场一步向上延拓垂向模量梯度场的递归余弦变换法，具有计算量小及算法精度较高等优点，解决现有技术由观测平面上的地磁异常场获得向上延拓平面上的地磁模量梯度场的三分量，需进行向上延拓和地磁异常场转换模量梯度场的两步运算问题；同时本发明仅需一步波谱变换和反变换就可同时得到延拓平面上的地磁模量梯度场的三分量，简化了位场处理的步骤，减少了位场处理的计算量；而且在进行波谱变换和反变换过程中采用递归余弦变换与反变换，在保证算法计算量小的同时，减小Gibbs边界效应，提高了由平面地磁异常场一步向上延拓平面地磁模量梯度场的转换精度。递归余弦变换与反变换可以适应任意长度的数据序列，快速傅里叶变换与反变换通常要求数据序列的长度是基数的整数次幂，补零有时会增大不少的计算量。

附图说明

图1是平面地磁异常场一步向上延拓平面地磁模量梯度场的余弦变换法的方框图。

图2是平面地磁异常场一步向上延拓平面地磁模量梯度场的余弦变换法的流程图。

图3(a)-图3(i)是球体模型数据对应的傅里叶变换与余弦变换一步向上延拓模量梯度图。其中：图3(a)ΔT_x(z₀-Δz)；图3(b)

图3(c)

图3(d)ΔT_y(z₀-Δz)；图3(e)

图3(f)

图3(g)ΔT_z(z₀-Δz)；图3(h)

图3(i)

图4(a)-图4(i)是球体和长方体混合模型数据对应的傅里叶变换与余弦变换一步向上延拓模量梯度图。其中：图4(a)ΔT_x(z₀-Δz)；图4(b)

图4(c)

图4(d)ΔT_y(z₀-Δz)；图4(e)

图4(f)

图4(g)ΔT_z(z₀-Δz)；图4(h)

图4(i)

具体实施方式

下面结合附图举例对本发明进行详细描述：

步骤1、平面地磁异常场测量数据的预处理。设定测量误差及地磁异常值的合理范围，对观测平面上的地磁异常场测量数据超出这一合理范围的粗大误差数据予以剔除。

步骤2、平面地磁异常场测量数据的网格化处理。选取网格化参数，利用径向基函数插值法，将观测平面上离散的地磁异常场测量数据进行网格化，获得网格化的地磁异常测量数据。

径向基函数插值法是一种精确度高的插值法，适用于对大量点数据进行插值计算，具有较高的预测精度，能较好地反映数据的变化；它是利用基函数确定已知数据点到内插网格节点的最佳权重。待插值点的径向基函数平面插值表达式为

式中，(x,y)为插值点的坐标，f₀(x,y)＝c_1xx+c_1yy+c₀，(x_j,y_j)为采样点的坐标，λ_j为对应于每一个采样点的权值，

为径向基函数，||·||为欧几里德范数。选取薄板样条函数作为径向基函数，令d_·j＝||x-x_j,y-y_j||，则有

由式(1)得到，采样点的已知值所给定的方程约束条件为

设f(x,y)具有二次连续导数，则能量函数表示为：

最小化能量函数，得到所需满足的正交条件为

为了求解式(1)中的系数c₀、c_1x、c_1y及权值λ_j，联立式(3)和式(5)得到关于系数c₀、c_1x、c_1y及权值λ_j的线性方程组为

线性方程组的系数矩阵是一个对称正定矩阵。求解得到系数c₀、c_1x、c_1y及权值λ_j，将其代入式(1)得到观测面上的网格化地磁异常场数据，n_x＝1,2,…,M_x，n_y＝1,2,…,M_y，观测面方程为z＝z₀，z轴垂直向下。

步骤3、利用递归余弦变换由ΔT(n_x,n_y,z₀)计算平面地磁异常场的余弦变换波谱

定义

步1)按式(8)，由ΔT(u_x,u_y,z₀)计算

式中，u_x＝n_x，θ_ux＝(u_x-1)π/M_x，

与

由递归公式(9)确定。

式中，m_x＝1,2,…,K_x，K_x＝[(M_x+1)/2]，符号[·]表示取整。

步2)按式(11)，由

计算ΔT^C(u_x,u_y,z₀)，

式中，u_y＝n_y，θ_uy＝(u_y-1)π/M_y，

与

由递归公式(12)确定。

式中，m_y＝1,2,…,K_y，K_y＝[(M_y+1)/2]，符号[·]表示取整。

步骤4、分别按式(14)、(15)和(16)，由观测面上的地磁异常场的余弦变换波谱ΔT^C(u_x,u_y,z₀)计算延拓平面上的地磁模量梯度场三分量的波谱

和

步骤5、对

进行行向的移相π/2的一维递归反余弦变换和列向的一维递归反余弦变换，得到延拓平面上的地磁模量梯度场的x分量ΔT_x(n_x,n_y,z₀-Δz)。

步1)按式(17)，由

计算

式中，θ_nx＝(n_x-1/2)π/M_x，

由递归公式(18)确定。

步2)按式(19)，由

计算ΔT_x(n_x,n_y,z₀-Δz)，

式中，θ_ny＝(n_y-1)π/M_y，

与

由递归公式(20)确定。

步骤6、对

进行行向的一维递归反余弦变换和列向的移相π/2的一维递归反正弦变换，得到延拓平面上的地磁模量梯度场的y分量ΔT_y(n_x,n_y,z₀-Δz)；

步1)按式(21)，由

计算

式中，θ_nx＝(n_x-1/2)π/M_x，

和

由递归公式(22)确定。

步2)按式(23)，由

计算ΔT_y(n_x,n_y,z₀-Δz)，

式中，θ_ny＝(n_y-1)π/M_y，

由递归公式(24)确定。

步骤7、对

进行行向和列向的一维递归反余弦变换，得到延拓平面上的地磁模量梯度场的z分量ΔT_z(n_x,n_y,z₀-Δz)。

步1)按式(25)，由

计算

式中，

与

由递归公式(26)确定。

步2)按式(27)，由

计算ΔT_z(n_x,n_y,z₀-Δz)，

式中，

与

由递归公式(28)确定。

利用递归余弦变换与反变换由平面地磁异常场一步向上延拓，得到平面地磁模量梯度场三分量，平面地磁异常场一步向上延拓平面地磁模量梯度场的余弦变换法的方框图如图1所示。平面地磁异常场一步向上延拓平面地磁模量梯度场的余弦变换法的算法流程图如图2所示。

为直接反映平面地磁异常场一步向上延拓平面地磁模量梯度场的三分量的位场转换效果，定义一个无量纲的相对误差指标

和

傅里叶变换与反变换下平面地磁模量梯度场j分量(j＝x,y,z)的位场转换相对误差

为

式中，ΔT_j(z₀-Δz)为平面地磁模量梯度场j分量的真实值，

为使用傅里叶变换与反变换得到的平面地磁模量梯度场j分量的延拓值，符号表示求平均值，符号| |表示求绝对值。

余弦变换与反变换下平面地磁模量梯度场j分量的位场转换相对误差

为

式中，

为使用余弦变换与反变换得到的平面地磁模量梯度场j分量的延拓值。

实验一：观测平面为z＝z₀＝0，观测区域的x和y坐标范围都为[-5110m,5110m]，网格数为256×256，向上延拓的高度为Δz＝-300m。六个球体磁源的位置分别为[0,0,500m]、[2000m,800m,400m]、[1800m,-1600m,500m]、[2000m,-1800m,600m]、[-1000m,2000m,500m]和[1500m,0,1000m]，磁化强度分别为1.6×10⁵A/m、1.4×10⁵A/m、1.3×10⁵A/m、1.5×10⁵A/m、1.2×10⁵A/m和1.8×10⁵A/m，球半径分别为10m、12m、8m、15m、20m和18m，磁倾角分别45°、60°、75°、45°、60°和25°，磁偏角分别为45°、35°、15°、30°、45°和45°。地磁场的磁倾角和磁偏角分别为60°和30°。没有磁测噪声时，延拓平面上的地磁模量梯度三分量的模型值、基于傅里叶变换与反变换的延拓值和基于余弦变换与反变换的延拓值如图3所示，计算得到两种变换下的相对误差分别是

和

这在理论上表明基于余弦变换与反变换的延拓地磁模量梯度精度高于基于傅里叶变换与反变换的精度。

观测面的磁测噪声为零均值、方差为10nT的高斯白噪声时进行50次蒙特卡洛仿真实验得到延拓的平均相对误差

和

球体模型数据的一步向上延拓地磁模量梯度的仿真实验表明，基于余弦变换与反变换的延拓地磁模量梯度精度高于基于傅里叶变换与反变换的精度。

实验二：在上述仿真实验区域添加两个长方体磁源。第一个长方体的中心位置为[0,400m,800m]，长、宽和高分别为30m、20m和15m，磁化强度为0.02A/m，磁倾角和磁偏角分别为45°和55°。第二个长方体的中心位置为[-1500m,1000m,600m]，长、宽和高分别为25m、18m和12m，磁化强度为0.012A/m，磁倾角和磁偏角分别为60°和40°。没有磁测噪声时，延拓平面上的地磁模量梯度三分量的模型值、基于傅里叶变换与反变换的延拓值和基于余弦变换与反变换的延拓值如图4所示，计算得到两种变换下的相对误差

和

这在理论上表明基于余弦变换与反变换的延拓地磁模量梯度精度高于基于傅里叶变换与反变换的精度。

观测面的磁测噪声为零均值、方差为10nT的高斯白噪声时进行50次蒙特卡洛仿真实验得到延拓的平均相对误差

和

球体和长方体模型数据的一步向上延拓地磁模量梯度的仿真实验同样表明，基于余弦变换与反变换的延拓地磁模量梯度精度高于基于傅里叶变换与反变换的精度。

本发明有益效果说明如下：

由于步骤5、步骤6和步骤7对应的三次波谱反变换在使用并行计算时可以同时进行，因此可以认为平面地磁异常场转换向上延拓平面地磁模量梯度场的一步法共需要一次波谱变换和一次波谱反变换。而通常的两种两步法，需要两次波谱变换和两次波谱反变换。分析比较表明了使用同种波谱变换条件下，一步法的计算量约为两步法计算量的一半。球体模型数据和球体与长方体混合模型数据仿真比较了基于傅里叶变换和余弦变换的一步向上延拓法的计算结果，结果表明余弦变换法减小Gibbs边界效应，提高了算法精度。递归余弦变换与反变换对数据序列的长度没有要求，可应用于任意长度的数据序列。

Claims (1)

1.一种平面地磁异常场一步向上延拓平面模量梯度场的递归余弦变换法，其特征是：

步骤1、选取网格化参数，对平面地磁异常场测量数据进行预处理；

步骤2、利用径向基函数插值法，将观测平面上的离散地磁异常场测量数据进行网格化，获得观测平面上的网格化地磁异常测量数据；

步骤3、先后对观测平面上的地磁异常场进行行向和列向的一维递归余弦变换，得到地磁异常场的二维余弦变换波谱；

步骤4、由观测平面上的地磁异常场波谱计算延拓平面上的地磁模量梯度场的水平x、水平y和垂直z波谱，具体包括：

按式

由观测平面上的地磁异常场的二维余弦变换波谱ΔT^C(u_x,u_y,z₀)计算延拓平面上的地磁模量梯度场的波谱

和

式中，R_1WC(u_x,u_y,Δz)为余弦变换下的由观测平面地磁异常场计算延拓平面地磁模量梯度场的一步波数域转换因子，其表达式为

步骤5、先后对延拓平面上的地磁模量梯度场的水平x波谱进行行向的移相π/2的一维递归反余弦变换和列向的一维递归反余弦变换，定义为行向正弦、列向余弦反变换，得到延拓平面上的地磁模量梯度场的水平x分量；

步骤6、先后对延拓平面上的地磁模量梯度场的水平y和波谱进行行向的一维递归反余弦变换和列向的移相π/2的一维递归反正弦变换，定义为行向余弦、列向正弦反变换，得到延拓平面上的地磁模量梯度场的水平y分量；

步骤7、先后对延拓平面上的地磁模量梯度场的垂直z波谱进行行向和列向的一维递归反余弦变换，得到延拓平面上的地磁模量梯度场的垂直z分量。

Images