一种基于空间域的重磁位场解析延拓方法

摘要

本发明公开了一种基于空间域的重磁位场解析延拓方法，重磁位场指的是重力场或者磁场；解析延拓指的是通过观测面上的重力场或者磁场计算其它曲面(包括平面)上的重力场或磁场的方法。解析延拓的求解有空间域和频率域两种思路，本发明是针对空间域的一种解析延拓新方法，可以直接对航磁测量的磁异常数据向下延拓至地面增强磁异常信号，也可以将地面观测数据向上延拓至更高的平面上以突出深部场特征。本发明可以将任意曲面观测的异常数据延拓至另一个任意曲面，突破了频率域方法中延拓面一定是平面的局限，解决了空间域解析延拓方法中的奇异点问题。

权利要求

1.一种基于空间域的重磁位场解析延拓方法，其特征在于，该方法包括以下步骤：(a)观测面S(x,y)为平面，其高程为z0是常数，在S(x,y)上观测的重力或磁异常记为U(α,β,z0)，其中α,β分别为观测点在x方向和y方向上的坐标；

(b)将观测的重力或磁异常数据U(α,β,z0)向上延拓到延拓面S(x,y,z1)，z1(x,y)表示延拓面S(x,y,z1)的高程；延拓高度表示为Δz(x,y)＝z1(x,y)-z0，则U(α,β,z0)向上延拓结果U(x,y,z1)用积分表示为：其中， 表示观测点(α,β)与延拓点(x,y)之

间的距离，延拓点是观测点在延拓面上的垂直投影；

(c)对观测面S(x,y)上的位场数据进行网格化，得到网格化数据，总共有M行N列，点距分别为Δx和Δy，第i行第j列节点上的位场异常值表示为Ui,j、坐标为(αi,βi)；

(d)由于观测数据的有限性，结合步骤(c)，将公式(1)近似为观测数据范围内的分段积分的累加：(e)由于观测数据的分辨率限制，假设观测数据在一个点距范围内是常数，则将公式(2)表示为：其中Um,n表示延拓结果的第m行第n列的值；

(f)将公式(3)中的二重积分表示为公式(4)并求解得到一个反正切形式的精确解Pm,n；i,j：(g)结合公式(3)和(4)将向上延拓表示为离散观测数据的加权求和，第m行第n列延拓结果表示为：(h)将公式(5)用矩阵乘法表示为：

U＝KU0      (6)

其中U0表示位场观测数据，U为延拓结果，K是向上延拓的核矩阵，其元素由公式(4)计算得到，表示为：对于平面到平面的向上延拓：如果延拓面是平面，Δz(x,y)＝h为延拓高度，是一个常数；则矩阵K所有元素均用第一行元素计算得到，只需要存储第一行元素即可；计算公式(6)所示的矩阵乘法即可得到平面到平面的向上延拓结果；

对于平面到曲面的向上延拓：如果延拓面是曲面，其高程为z1(x,y)，则延拓高度Δz(x,y)由延拓点高程与对应的观测点高程相减得到，然后带入公式(4)即可计算得到矩阵K，计算公式(6)所示的矩阵乘法即可得到平面到曲面的向上延拓结果。

2.根据权利要求1所述的一种基于空间域的重磁位场解析延拓方法，其特征在于，该方法还用于实现向下延拓，向下延拓是向上延拓的反问题，对于平面到平面的向下延拓及曲面到平面的向下延拓：根据延拓高度和公式(4)和公式(7)计算核矩阵K，然后利用Landweber迭代法和光滑-拟合曲线求解式(6)得到向下延拓结果。

3.根据权利要求2所述的一种基于空间域的重磁位场解析延拓方法，其特征在于，所述利用Landweber迭代法和光滑-拟合曲线求解式(6)具体为：Landweber迭代法的最佳迭代次数通过光滑-拟合曲线的第一个极小值确定；光滑-拟合曲线中的“光滑”指的是向下延拓结果的水平梯度的二范数，作为纵坐标；“拟合”指的是向下延拓结果与观测数据之差的二范数与观测数据二范数的比值，作为横坐标；对不同迭代次数的结果均可计算出“光滑”和“拟合”值，然后绘制曲线，找出第一个极小值点对应的迭代次数就是最优的迭代次数，最优迭代次数对应的向下延拓结果就是所求解的向下延拓结果。

4.根据权利要求1所述的一种基于空间域的重磁位场解析延拓方法，其特征在于，该方法还用于实现曲面到曲面之间的解析延拓：曲面到曲面之间的解析延拓包括曲面到曲面的向下延拓及曲面到曲面的向上延拓；如果观测面是曲面，延拓面也是曲面，则首先需要求解观测曲面到最低的一个平面上的延拓场；如果是曲面到曲面的向上延拓，这个最低的平面取观测曲面的最低点所在平面；如果是曲面到曲面的向下延拓，这个最低的平面取为延拓曲面的最低点所在平面；将观测的位场向下延拓至最低平面得到等效场Ue；然后计算最低平面上的等效场Ue到延拓曲面的向上延拓，最终实现曲面到曲面的解析延拓。

说明书

一种基于空间域的重磁位场解析延拓方法

技术领域

[0001] 本发明属于重力、磁异常数据处理技术领域，涉及一种基于空间域的重磁位场解析延拓方法。

背景技术

[0002] 重磁勘探方法在地面、航空、海面及海底进行重力场或磁场观测，然后将原始数据经过一定校正得到重力异常或磁异常，反应的是观测面以下的密度和磁性不均匀性，从而应用于地质填图、矿藏勘探、管线探测等。重力勘探和磁法勘探这两种方法具有一定的相似性，故通常将重力异常和磁异常一起讨论，称之为重磁异常。重力场和磁场均是位场，故又称之为重磁位场。而重磁位场的解析延拓是对重磁异常数据处理的重要方法，通过解析延拓可以反映出重磁异常的多尺度信息。向上延拓可以突出深部场源信息，压制浅部场源的干扰；相反地，向下延拓则是增强浅部场源信息。

[0003] 传统的解析延拓方法是基于快速傅里叶变换在频率域进行计算，这种方法的优势是计算速度快，但是其缺点是只要求重磁异常的观测面是平面。空间域方法可以弥补此缺点，但传统的空间域延拓方法只适用于延拓高度大于一个点距的情况，而且不适用于计算曲面之间的延拓。传统空间域延拓方法在计算较小延拓高度的情况时，会出现奇异值。但是通常的观测面均是曲面，比如航磁观测的飞行器大多是定高飞行，换言之，航磁观测曲面近似于地形平行。这种情况下的解析延拓都会遇到延拓高度小于一个点距和曲面之间延拓的问题。

发明内容

[0004] 为了解决上述一般情况下传统解析延拓方法的不适用问题，就需要从解析延拓的理论公式中寻找解决方法。本发明通过对解析延拓的二重积分公式进行分段积分，求解得到了空间域解析延拓的一个新的核函数。通过新的核函数可以方便地实现平面到平面的向上延拓、平面到曲面的向上延拓、平面到平面的向下延拓、曲面到平面的向下延拓以及场源外的两个曲面之间的延拓。

[0005] 本发明通过以下技术方案得以实现：一种基于空间域的重磁位场解析延拓方法，该方法包括以下步骤：

[0006] (a)观测面S(x,y)为平面，其高程为z0是常数，在S(x,y)上观测的重力或磁异常记为U(α,β,z0)，其中α,β分别为观测点在x方向和y方向上的坐标。

[0007] (b)将观测的重力或磁异常数据U(α,β,z0)向上延拓到延拓面S(x,y,z1)，z1(x,y)表示延拓面S(x,y,z1)的高程；延拓高度表示为Δz(x,y)＝z1(x,y)-z0，则U(α,β,z0)向上延拓结果U(x,y,z1)用积分表示为：

[0008]

[0009] 其中， 表示观测点(α,β)与延拓点(x,y)之间的距离，延拓点是观测点在延拓面上的垂直投影。

[0010] (c)对观测面S(x,y)上的位场数据进行网格化，得到网格化数据，总共有M行N列，点距分别为Δx和Δy，第i行第j列节点上的位场异常值表示为Ui,j、坐标为(αi,βi)。

[0011] (d)由于观测数据的有限性，结合步骤(c)，将公式(1)近似为观测数据范围内的分段积分的累加：

[0012]

[0013] (e)由于观测数据的分辨率限制，假设观测数据在一个点距范围内是常数，则将公式(2)表示为：

[0014]

[0015] 其中Um,n表示延拓结果的第m行第n列的值；

[0016] (f)将公式(3)中的二重积分表示为公式(4)并求解得到一个反正切形式的精确解Pm,n；i,j：

[0017]

[0018]

[0019] (g)结合公式(3)和(4)将向上延拓表示为离散观测数据的加权求和，第m行第n列延拓结果表示为：

[0020]

[0021] (h)将公式(5)用矩阵乘法表示为：

[0022] U＝KU0    (6)

[0023] 其中U0表示位场观测数据，U为延拓结果，K是向上延拓的核矩阵，其元素由公式(4)计算得到，表示为：

[0024]

[0025] 对于平面到平面的向上延拓：如果延拓面是平面，Δz(x,y)＝h为延拓高度，是一个常数；则矩阵K所有元素均可以用第一行元素计算得到，只需要存储第一行元素即可；计算公式(6)所示的矩阵乘法即可得到平面到平面的向上延拓结果；

[0026] 对于平面到曲面的向上延拓：如果延拓面是曲面，其高程为z1(x,y)，则延拓高度Δz(x,y)可以由延拓点高程与对应的观测点高程相减得到，然后带入公式(4)即可计算得到矩阵K，计算公式(6)所示的矩阵乘法即可得到平面到曲面的向上延拓结果。

[0027] 进一步地，该方法还可以用于实现向下延拓，向下延拓是向上延拓的反问题，对于平面到平面的向下延拓及曲面到平面的向下延拓：

[0028] 根据延拓高度和公式(4)和公式(7)计算核矩阵K，然后利用Landweber迭代法和光滑-拟合曲线求解式(6)得到向下延拓结果。

[0029] 进一步地，所述利用Landweber迭代法和光滑-拟合曲线求解式(6)具体为：

[0030] Landweber迭代法的最佳迭代次数通过光滑-拟合曲线的第一个极小值确定；光滑-拟合曲线中的“光滑”指的是向下延拓结果的水平梯度的二范数，作为纵坐标；“拟合”指的是向下延拓结果与观测数据之差的二范数与观测数据二范数的比值，作为横坐标；对不同迭代次数的结果均可计算出“光滑”和“拟合”值，然后绘制曲线，找出第一个极小值点对应的迭代次数就是最优的迭代次数，最优迭代次数对应的向下延拓结果就是所求解的向下延拓结果。

[0031] 进一步地，该方法还可以用于实现曲面到曲面之间的解析延拓：

[0032] 曲面到曲面之间的解析延拓包括曲面到曲面的向下延拓及曲面到曲面的向上延拓；如果观测面是曲面，延拓面也是曲面，则首先需要求解观测曲面到最低的一个平面上的延拓场；如果是曲面到曲面的向上延拓，这个最低的平面取观测曲面的最低点所在平面；如果是曲面到曲面的向下延拓，这个最低的平面取为延拓曲面的最低点所在平面；将观测的位场向下延拓至最低平面得到等效场Ue；然后计算最低平面上的等效场Ue到延拓曲面的向上延拓，最终实现曲面到曲面的解析延拓，这个过程的理论依据是重磁位场的等效源原理。

[0033] 本发明的有益效果是：本发明方法可以直接对航磁测量的磁异常数据向下延拓至地面增强磁异常信号，也可以将地面观测数据向上延拓至更高的平面上以突出深部场特征。本发明可以将任意曲面观测的异常数据延拓至另一个任意曲面，突破了频率域方法中延拓面一定是平面的局限，解决了空间域解析延拓方法中的奇异点问题。

附图说明

[0034] 图1是本发明方法的流程图。

[0035] 图2是本发明方法的网格节点及积分示意图。

[0036] 图3是本发明方法模型实验所用的理论模拟数据；(a)为模拟的地形，黑色方框为立方体模型在平面上的投影；(b)为高程z＝0平面上的磁异常；(c)为图(a)所示的模拟地形上的磁异常；(d)为高程z＝8m平面上的磁异常；(e)为这三个磁异常数据和地形以及模型的三维视图。

[0037] 图4是本发明方法对磁异常从平面到平面的向上延拓实例结果；(a)～(c)是用频率域方法计算的结果；(d)～(f)是本发明方法的计算结果；(a)和(d)是对图3(b)所示的磁异常向上延拓至z＝8m的平面的结果；(b)和(e)分别为向上延拓结果(a)和(d)与图3(d)所示的理论数据之间的误差分布；(c)和(f)分别是对误差分布(b)和(e)进行统计的柱状图，误差均方根(RMSE)分别为2.13和3.40nT。

[0038] 图5是本发明方法对磁异常从平面向上延拓至曲面的实例结果；(a)是对图3(a)所示的磁异常向上延拓至图3(a)所示的模拟地形上的结果；(b)为向上延拓结果(a)与图3(c)所示的理论模拟数据之间的误差分布；(c)为对误差分布(b)进行统计的柱状图，误差均方根(RMSE)为1.76nT。

[0039] 图6是本发明方法对磁异常从平面向下延拓至平面的结果；(a)是对图3(d)所示的磁异常向下延拓至z＝0m平面的结果；(b)为向下延拓结果(a)与图3(b)所示的理论模拟磁异常之间的误差分布；(c)为对误差分布(b)进行统计的柱状图，误差均方根(RMSE)为3.4nT。

[0040] 图7是本发明方法对磁异常从曲面向下延拓至平面的结果；(a)是对图3(c)所示的曲面上的磁异常向下延拓至z＝0m平面的结果；(b)向下延拓结果与图3(b)所示的理论模拟数据之间的误差分布；(c)对误差分布进行统计的柱状图，均方根为0.69nT。

[0041] 图8是本发明方法对含有噪声的磁异常从曲面向下延拓至平面的结果；(a)对图3(c)所示的曲面上的磁异常加入高斯白噪声(噪声均值为0nT，标准差为1nT)后的异常分布；(b)白噪声的分布；(c)对此白噪声进行统计的柱状图；(d)将(a)所示的含有噪声的曲面上的磁异常向下延拓至z＝0m平面的结果；(e)向下延拓结果与图3(b)所示的理论模型数据之间的误差分布；(f)对误差分布进行统计的柱状图，均方差为5.29nT。

[0042] 图9是本发明在计算磁异常向下延拓(对应图8)的最优迭代次数选择曲线。

[0043] 图10是本发明方法对实测航磁异常进行解析延拓的结果。(a)航磁异常；(b)航磁观测的飞行高度；(c)对观测的航次异常向下延拓至高程为970m的平面上的结果；(d)对(c)所示的延拓延拓场向上延拓至高程z＝2200m的结果，表示对观测磁异常的一个曲化平(从曲面到平面的延拓)的结果；(e)将(c)所示延拓场向上延拓至观测曲面之上200m曲面(与观测面平行)上的结果，表示从曲面到曲面的向上延拓。

具体实施方式

[0044] 下面结合实施例对本发明的可行性做进一步的详细说明，以下实施例是对本发明的解释而本发明并不局限于以下实施例的结果。

[0045] 本发明提供的一种基于空间域的重磁位场解析延拓方法，可实现平面到平面的向上延拓、平面到曲面的向上延拓、平面到平面的向下延拓、曲面到平面的向下延拓以及场源外的两个曲面之间的延拓；

[0046] (1)向上延拓，包括如下步骤：

[0047] (a)观测面S(x,y)为平面，其高程为z0是常数，在S(x,y)上观测的重力或磁异常(以下简称为位场)记为U(α,β,z0)，其中α,β分别为观测点在x方向和y方向上的坐标；

[0048] (b)将观测的重力或磁异常数据U(α,β,z0)向上延拓到延拓面S(x,y,z1)，z1(x,y)表示延拓面S(x,y,z1)的高程。延拓高度可以表示为Δz(x,y)＝z1(x,y)-z0，则U(α,β,z0)向上延拓结果U(x,y,z1)可以用积分表示为：

[0049]

[0050] 其中， 表示观测点(α,β)与延拓点(x,y)之间的距离，延拓点是观测点在延拓面上的垂直投影；

[0051] (c)公式(1)为连续积分，但是实际观测数据是离散点，所以需要对观测面S(x,y)上的位场数据进行网格化，得到网格化数据，总共有M行N列，点距分别为Δx和Δy，第i行第j列节点上的位场异常值表示为Ui,j、坐标为(αi,βi)；

[0052] (d)由于观测数据的有限性，结合步骤(c)，则可以将公式(1)近似为观测数据范围内的分段积分的累加：

[0053]

[0054] (e)由于观测数据的分辨率限制，假设观测数据在一个点距范围内是常数，则公式(2)可以表示为：

[0055]

[0056] 其中Um,n表示延拓结果的第m行第n列的值；

[0057] (f)公式(3)中的二重积分是第i行第j列观测点和第m行第n列延拓点的坐标的函数，因此可以将此二重积分表示为公式(4)并求解得到一个反正切形式的精确解Pm,n；i,j：

[0058]

[0059] (g)结合公式(3)和(4)可以将向上延拓表示为离散观测数据的加权求和，第m行第n列延拓结果可以表示为：

[0060]

[0061] 其中Pm,n；i,j可视为加权系数；

[0062] (h)将公式(5)可以进一步地用矩阵乘法表示为：

[0063] U＝KU0    (6)

[0064] 其中U0表示位场观测数据，是一个M×N的列向量；U为延拓结果，也是一个M×N的列向量；K是向上延拓的核矩阵，是一个(M×N)×(M×N)的方阵，其元素由公式(4)计算得到，表示为：

[0065]

[0066] (1.1)平面到平面的向上延拓

[0067] 延拓面是平面的情况下，Δz(x,y)＝h为延拓高度，是一个常数。则矩阵K(公式7)所有元素都可以用第一行元素计算得到，所以只需要存储第一行元素即可，减少计算机内存消耗。计算公式(6)所示的矩阵乘法即可得到平面到平面的向上延拓结果；

[0068] (1.2)平面到曲面的向上延拓

[0069] 如果延拓面是曲面，其高程为z1(x,y)，则延拓高度Δz(x,y)可以由延拓点高程与对应的观测点高程相减得到，然后带入公式(4)即可计算得到矩阵K，然后计算公式(6)所示的矩阵乘法即可得到平面到曲面的向上延拓结果；

[0070] (2)向下延拓，包含如下步骤：

[0071] (A)向下延拓是向上延拓的反问题，上述的向上延拓原理是向下延拓的基础。如果U是观测数据，则求解式(6)表示的线性方程组，可以得到向下延拓结果U0；

[0072] (B)采用Landweber迭代法求解线性方程组(6)，最佳迭代次数可以通过光滑-拟合曲线的第一个极小值确定。光滑-拟合曲线中的“光滑”指的是向下延拓结果的水平梯度的二范数，作为纵坐标；“拟合”指的是向下延拓结果与观测数据之差的二范数与观测数据二范数的比值，作为横坐标。对不同迭代次数的结果均可计算出“光滑”和“拟合”值，然后绘制曲线，找出第一个极小值点对应的迭代次数就是最优的迭代次数，最优迭代次数对应的向下延拓结果就是所求解的向下延拓结果。

[0073] 平面到平面的向下延拓及曲面到平面的向下延拓：

[0074] 根据延拓高度和公式(4)和公式(7)计算核矩阵K，然后利用Landweber迭代法和光滑-拟合曲线求解线性方程组(6)得到最优向下延拓结果；

[0075] (3)曲面到曲面之间的解析延拓

[0076] 曲面到曲面之间的解析延拓包括曲面到曲面的向下延拓及曲面到曲面的向上延拓。如果观测面是曲面，延拓面也是曲面，则首先需要求解观测曲面到最低的一个平面上的延拓场。如果是曲面到曲面的向上延拓，这个最低的平面取观测曲面的最低点所在平面；如果是曲面到曲面的向下延拓，这个最低的平面取为延拓曲面的最低点所在平面。与“平面到平面的向下延拓及曲面到平面的向下延拓”相同的方法，将观测的位场向下延拓至最低平面得到等效场Ue；然后利用与(1.2)相同的方法计算最低平面上的等效场Ue到延拓曲面的向上延拓，最终实现曲面到曲面的解析延拓。

[0077] 实施示例1：对模拟计算的磁异常进行解析延拓

[0078] 利用立方体模型的磁异常正演公式分别计算得到一个立方体(图3e)在平面z＝0m，平面z＝8m以及一个曲面(图3a)上的磁异常，三个磁异常如图3(b)、(d)和(c)所示。利用上述的重磁位场解析延拓方法分别对模拟磁异常数据进行向上和向下延拓，然后与精确解进行对比，验证方法的可行性和准确性。

[0079] (1)磁异常从一个平面向上延拓至另一平面

[0080] 对图3b所示的z＝0m平面上的磁异常数据向上延拓8m，结果如图4a所示，将向上延拓结果与z＝8m平面上的真实磁异常作差计算其误差，如图4b所示；对误差进行统计，如图4c所示。均方差为2.13nT。图4d-f为频率域计算结果，对比可发现空间域的计算结果更精确。

[0081] (2)磁异常从平面向上延拓至曲面

[0082] 对图3b所示的z＝0m平面上的磁异常向上延拓至图3a所示的曲面上，结果如图5a所示，将向上延拓结果与曲面上的精确磁异常进行作差，计算误差分布如图5b所示。对误差统计得到均方差为1.76nT，如图5c所示。

[0083] (3)磁异常从平面向下延拓至平面

[0084] 对图3d所示的z＝8m平面上的磁异常向下延拓8m，结果如图6a所示。将向下延拓结果与z＝0m平面上的精确磁异常进行作差，计算其误差分布，如图6b所示。对误差统计的到均方差为2.4nT，如图6c所示。

[0085] (4)磁异常从曲面向下延拓至平面

[0086] 对图3c所示的曲面上的磁异常向下延拓至z＝0m平面上，延拓结果如图7a所示。将延拓结果与z＝0m平面精确磁异常作差，计算其误差分布，如图7b所示。对误差统计得到均方差为0.69nT，如图7c所示。

[0087] (5)对含有噪声的磁异常进行向下延拓

[0088] 在图3c所示的磁异常中加入均值为0nT，标准差为1nT的高斯白噪声(如图8b所示，噪声的统计如图8c所示，加入白噪声的磁异常结果如图8a所示。对图8a所示的含噪声磁异常向下延拓至z＝0m平面，结果如图8d所示。对结果和z＝0m平面精确磁异常进行作差计算其误差分布，如图8e所示。对误差进行统计得到均方差为5.29nT，如图8f所示。

[0089] (6)最优迭代次数的选择

[0090] 本发明在计算向下延拓时，需要进行迭代，最优迭代次数的选择依据提出的拟合-光滑曲线的第一个极小值点。以含噪声数据的向下延拓(图8)为例，其向下延拓最优迭代次数的选择所依据的拟合-光滑曲线如图9所示，最优迭代次数为277。

[0091] 实施示例2：实测航磁数据的解析延拓

[0092] 将本发明应用到航磁异常的解析延拓处理中，航磁异常如图10a所示，飞行曲面的高程如图10b所示。对于曲面上观测的磁异常，进一步处理之前，需要将其向下延拓至一个平面上。这个平面一般选择为飞行高度曲面的最低点所在平面，这里选择为z＝970m。将航磁异常延拓到z＝970m平面的结果如图10c所示，此向下延拓的高度并非一个常数，而是0～1430m变化的，传统的空间域方法无法实现。然后利用此中间结果进行其他平面之间的延拓处理，对图10c所示的磁异常向上延拓至z＝2200m平面得到航磁异常的一个曲化平的结果，如图10d所示。有些情况下，可能并不需要曲化平，而只需要将曲面上的磁异常延拓至另一个曲面上，比如延拓至于观测曲面平行的一个曲面。将图10c所示的结果向上延拓至与飞行曲面之上(平行)200m的结果如图10e所示。

[0093] 以上两个实施示例，证明了新的空间域延拓方法在实际应用中的可行性和准确性。对重磁位场异常进行解析延拓或者曲化平，对异常进一步解释和处理具有重要实用价值。