具体实施方式

以下结合说明书附图和具体优选的实施例对本发明作进一步描述，但并不因此而限制本发明的保护范围。

如图1所示，本实施例基于地磁数据的垂直高度位场测量方法的步骤包括：

S1.频谱分析：获取观测平面的地磁位场数据并进行频谱分析，得到观测平面的频谱；

S2.自适应滤波算子构建：根据观测平面的频谱确定截止频率，并根据截止频率按照低通滤波器模型构建自适应滤波算子；

S3.向下延拓：以观测平面的频谱作为初值进行向下延拓的迭代更新，每次迭代时将当前延拓平面的频谱向上延拓到观测平面得到延拓值，根据观测平面的频谱与延拓值计算得到剩余谱，使用自适应滤波算子处理剩余谱并更新延拓值，直至得到最终的延拓平面的频谱输出；

S4.测量结果输出：将步骤S3最终得到的延拓平面的频谱进行傅里叶反变换，得到空间域位场的测量结果输出。

本实施例通过依据地磁位场数据的频谱分析构建自适应滤波算子，在向下延拓的迭代过程中利用自适应滤波算子处理剩余谱，可以达到低通滤波效果，自适应的调整滤波算子的低通频带，并尽可能利用数据中高频有用信息而抑制噪声部分，提高向下延拓迭代的收敛速度，能够快速、精确的得到地磁数据空间域位场，解决传统地磁数据库中在不同高度平面上的数据不完备性、精度低的问题。

上述在向下延拓的迭代过程中，还可以利用自适应滤波算子的低通特性延缓发散过程，提高迭代的稳定性，解决传统采用半收敛算子的迭代方法中会随着迭代次数增加而逐渐发散的问题，保证在较长迭代区间都可以收敛到最佳延拓效果。

本实施例中步骤S1的步骤包括：

S11.建立三维直角坐标系，其中地平面为Z＝0平面，垂直地平面向下为Z轴正方向，观测平面为T(z₀)，观测平面的地磁位场数据为

S12.将观测平面的地磁位场数据向下延拓到靠近场源的平面T(z₁)，得到延拓平面位场数据向下延拓高度为h＝z₀-z₁，z₀为观测平面的高度，z₁为靠近场源的平面的高度；

S13.分别将观测平面的地磁位场数据延拓平面位场数据以及空间域向上延拓核函数K(x,y)进行二维傅里叶变换，得到观测平面的频谱延拓平面的频谱以及向上延拓算子H_up(k_x,k_y)。

上述地磁位场数据延拓平面位场数据以及空间域向上延拓核函数K(x,y)之间关系可表示为：

其中，x,y为已知观测平面位场的位置坐标，ξ,γ为延拓平面位场的位置坐标。

将观测平面的地磁位场数据延拓平面位场数据进行二维傅里叶变换得到相应的频谱为：

应用积分变化表，得到频率域向上延拓算子H_up(k_x,k_y)即为：

通过傅里叶变换，把空间域的卷积运算转换到频率域的乘积运算，可得：

进一步可以简化为：

F＝H_up·G (8)

通过上述变换，可以将向下延拓问题转换为相当于已知H_up、F，求解G的解线性方程组的过程，从而使得可以与最小二乘问题联系在一起，即找到G^*使得：

后续经过迭代求解后即可得到上述G。

迭代求解前需先确定截止频率，本实施例步骤S2中确定截止频率的步骤包括：

S21.将观测平面的频谱零频移位到中心计算得到位场功率谱P(k_x,k_y)；

S22.以位场功率谱P(k_x,k_y)的中心作圆形成一系列环带，取环带内所有点的功率值作平均，得到径向平均功率谱；

S23.根据径向平均功率谱随频率的衰减程度变化确定出截止频率w_c。

由于功率谱的逆傅里叶变换是自相关谱，而观测数据中的高频噪声与位场信号是不相关的，但是由于混叠着小部分位场信号的高频成分，在功率谱的衰减应该是较为平缓，而位场信号中大部分的高低频信息会导致功率谱较大的衰减程度，即通过分析平均径向功率谱随频率的衰减程度，可以定性的找到一个区分位场信号与噪声的截止频率，同时能够尽可能多的利用到高频位场信号。

在具体应用实施例中，由式(6)计算频率域的向上延拓算子H_up(k_x,k_y)，以及由式(4)得到后，将频谱零频移位到中心计算其功率谱P(k_x,k_y)，即：

其中M，N为网格化数据的行数、列数。

本实施例进一步以位场功率谱P(k_x,k_y)的中心作圆，圆的半径w_r分别为基频的整数倍，从而形成一系列环带；取环带内所有点的功率值作平均，即得到径向平均功率谱然后分析观测平面位场的频率分布，以径向频率w_r为横坐标，对径向平均功率谱取对数为纵坐标，可得到对数径向平均功率谱曲线根据对数径向平均功率谱曲线的衰减程度差异确定截止频率w_c，具体将曲线衰减程度出现差异的拐点的横坐标作为截止频率w_c。

本实施例中自适应滤波算子通过使用低通滤波器模型，并使用截止频率对低通滤波器模型进行自适应调整构建得到。自适应滤波算子φ(w)实质上是低通滤波器，且可以根据位场数据的频率成分自适应的调整通带宽度，实现对高频噪声的压制，保证算法稳定性，且调整截止频率可以改变低通频带范围，使之与位场的频率成分相对应，建立基于数据的自适应滤波算子，充分利用高频有用信号。

在具体应用实施例中，先设计低通滤波算子如下，保证迭代算法稳定性。

利用截止频率w_c对低通滤波器φ进行自适应调整，建立基于数据-截止频率-低通滤波算子的对应关系。

综合截止频率w_c可调性与标准低通滤波算子φ(w)，得到最终的自适应滤波算子，

其中，w_c为截止频率，w为频率，n可取2或3。

采用上述自适应滤波算子，能够自适应的调整滤波算子的低通频带，尽可能利用数据中高频成分而压制噪声部分，提高向下延拓迭代收敛速度与实用性。

本实施例中步骤S3的步骤包括：

S31.把观测平面的位场频谱F(k_x,k_y)作为向下延拓的延拓平面频谱的初值G₀(k_x,k_y)，即：

G₀(k_x,k_y)＝F(k_x,k_y) (15)

S32.将当前延拓平面频谱G_i(k_x,k_y)向上延拓到观测平面得到延拓值，i＝0,1,…n，并将观测平面的位场频谱F(k_x,k_y)与延拓值作差得到剩余谱e，即：

e＝F(k_x,k_y)-H_up(k_x,k_y)·G_i(k_x,k_y) (16)

其中，H_up(k_x,k_y)为将空间域向上延拓核函数K(x,y)进行二维傅里叶变换得到的向上延拓算子；

S33.使用自适应滤波算子处理剩余谱e得到修正项，具体可按照下式计算得到：

ΔG＝φ·e (17)

S34.使用修正项更新当前延拓平面的频谱，即：

G_i+1(k_x,k_y)＝G_i(k_x,k_y)+ΔG (18)

S35.判断是否达到停止迭代条件，如果是则得到最终的延拓平面频谱G_n(k_x,k_y)输出，否则返回步骤S32。

剩余谱中包含位场数据中的中低频成分、高频成分与噪声部分，最佳效果是保留有用信息而滤除噪声，但是高频成分与噪声有混叠区域。本实施例利用自适应滤波算子处理剩余谱，能够达到低通滤波效果，在充分利用中低频信息的前提下，能够尽可能保留高频有用信号，抑制噪声对迭代过程的影响，提高稳定性，同时调整低通频带充分利用高频有用信息，提高收敛速度。

上述步骤S35中具体当迭代次数达到预设最大迭代次数或迭代结果误差小于预设阈值(|e|<ε(ε为很小的数)时，判定达到停止迭代条件。

本实施例中将迭代结束后得到的频谱后，具体按照下式通过傅里叶反变换得到待延拓平面位场数据：

磁源埋深体现在空间域为梯度变化大小，转化到频率域对应频率成分高低，磁源埋深浅的区域对应于空间中梯度较大，相当于频率域中的高频成分。本实施例通过利用径向平均功率谱分析磁数据的频率成分分布确定截止频率，基于截止频率构建自适应滤波算子，以通过对数据频率成分的分析，自适应的调整滤波算子的低通频带，使得能够尽可能利用数据中高频成分而压制噪声部分，提高迭代过程的收敛速度与稳定性，且能够保证在较长迭代区间都可以收敛到最佳延拓效果，从而基于数据自适应调整可以实现收敛速度快、稳定性好、延拓精度高的向下延拓迭代。

为验证本发明的有效性，在具体应用实施例中在向下延拓迭代过程中分别利用本发明自适应滤波算子与传统向上延拓算子并比较结果，不同截止频率下的观测平面位场数据的功率谱曲线，如图3所示，本发明自适应滤波算子与传统向上延拓算子虽然都是低通滤波器，但是本发明自适应滤波算子不仅具有低频通带，过滤带窄，且可以根据位场数据的频率成分自适应的调整通带宽度(截止频率越大，通带越长，可以使用的频率成分也就越多)，有效地利用频谱的中高频分量，提升迭代过程的性能，尤其适用于浅源位场的延拓。同时，低通特性还可以很好压制高频噪声的影响，提高算法稳定性。

为验证本发明上述方法的有效性，进一步在具体应用实例中使用本发明上述方法进行仿真实验，仿真条件具体为：

(1)采用的球体的磁异常表达式为：

其中μ₀表示真空磁导率；D表示磁化强度的偏角；I表示地磁场倾角，也表示磁化强度倾角；m表示磁矩，它与磁化率J之间的关系：

m＝J×πR³ (18)

其中，R表示球体半径。

(2)组合球体模型的球体磁源数量为5个，参数如表1。

表1 组合球体模型参数

接下来，通过设置两组不同埋深的仿真试验，以验证采用本发明方法可以根据数据频率成分的差异自适应调整，得到最佳延拓效果，并且保证收敛速度与稳定性。

本实施例定义三个精度指标来定量分析本发明的性能：

1)均方根误差Δg_rmse

2)误差矩阵Δg_me

Δg_me＝|g_e-g_t| (20)

3)相对误差变化率Δg_re

进一步本实施例对两组不同埋深模型的观测平面进行谱分析，不同埋深模型的观测平面位场数据的功率谱曲线如图4所示。

根据对数功率谱衰减程度对应于数据中不同频率成分，按照本发明方法可以确定截止频率为：

1)埋深为200m的浅源模型截止频率为0.0041，

2)埋深为1000m的深源模型的截止频率为0.0023

从上述可看出，浅源相对于深源，频率成分中高频占比增大，对数功率谱衰减程度的拐点越靠右。

进一步本实施例将球体磁源埋深统一设置为200m，得到的结果如图5所示，其中图5(a)是Z＝200m平面的位场，图5(b)是Z＝700m平面的位场。图5(c)是本发明延拓结果，图5(d)是延拓结果图5(c)与理论平面图5(a)作差得到的误差矩阵，图5(d)是迭代过程的均方根误差曲线，分析各结果图可得到：

1)由图5(c)的延拓结果与理论平面相比，可以直观的看出本发明取得极好的延拓效果，完全拟合到理论平面，仅仅存在一些边界的吉布斯效应，可以通过扩边方法进一步提升。

2)从图5(d)的误差矩阵结果中可以直观地看到延拓误差的分布，即本发明方法在大部分梯度变化小的区域以及磁源主体区域的延拓误差几乎为0，误差仅仅分布在磁源主体的中心部分(频率成分极高)以及由于磁源靠近边界导致的褶皱部分。

3)从图5(e)的均方根误差曲线中看得出，本发明不仅具有较快收敛速度，而且保证在较长迭代区间内保持最佳延拓效果，稳定性极好，可以抑制迭代算法由于本身的半收敛性质导致在迭代后期迅速发散。

进一步本实施例将球体磁源埋深统一设置为1000m，结果如图6所示，其中图6(a)是Z＝200m平面的位场，图6(b)是Z＝700m平面的位场。图6(c)是本发明延拓结果，图5(d)是延拓结果图6(c)与理论平面图6(a)作差得到的误差矩阵，图6(d)是迭代过程的均方根误差曲线，分析各结果图可得到：

1)埋深1000m的位场数据梯度变化较小，大部分为中低频信息，可得本发明(图6(c))几乎完全与理论平面相同，延拓效果相比埋深浅模型更好。

2)通过图6(d)的误差矩阵可以看出，误差分布情况与埋深浅模型相同，但是误差范围控制在0.5nT以内，具有极小的延拓误差，在低频数据延拓中展示出极高的延拓精度。

3)通过如图6(e)的均方根误差曲线可以看出，深源模型中表现出较好的延拓精度与稳定性，本发明在此基础上还进一步提高了收敛速度。

本实施例进一步分析不同埋深下的延拓效果，表2给出了延拓精度指标：最小均方根误差、最大误差矩阵值以及最小相对误差率。

表2 本发明延拓精度指标

本实施例通过两组不同埋深的模型仿真试验，其中埋深200m相当于频率成分中高频分量占比大，埋深1000m相当于频率分量主要以低频成分为主，从表中看出埋深代表的频率成分对延拓效果有很大影响，高频分量越多的位场数据，向下延拓的误差越大，即根据数据频率成分来自适应调整迭代算法是十分必要的，相对于传统迭代方法，本发明在对不同频率成分的磁数据都能够表现出收敛速度快、延拓精度高以及较强稳定性等优势。

本实施例进一步实测数据采用某地区航测的磁异常数据，对实测数据的延拓效果如图7所示，其中图7(a)是航测平面的位场数据，图7(b)是将图7(a)向上延拓500m后的位场数据。图7(c)是采用本发明方法将图7(b)向下延拓500m后的延拓结果，可以看到本发明在实际使用中也可以得到较好的延拓效果，延拓结果与理论平面比较误差较小，对位场中梯度变化复杂区域也有较好的延拓结果。

上述只是本发明的较佳实施例，并非对本发明作任何形式上的限制。虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明。因此，凡是未脱离本发明技术方案的内容，依据本发明技术实质对以上实施例所做的任何简单修改、等同变化及修饰，均应落在本发明技术方案保护的范围内。