一种重力场水平分量的获取方法

摘要

针对传统基于快速傅里叶变换算法的波数域方法因波数域算子在零波数存在奇异而导致计算精度低的问题，本发明公开了一种重力场水平分量的获取方法。本发明通过设定应用场景、加权系数计算，和一种二维离散卷积快速算法，实现重力场垂直分量到水平分量的快速、高精度计算，实现了大规模重力测量数据条件下分量转换在计算效率和计算精度上的统一。本发明解决了传统波数域方法因波数域转换算子在零波数存在奇异而导致水平分量计算精度低的问题，为满足大规模重力数据处理、反演和解释等提供了方法支撑。

权利要求

1.一种重力场水平分量的获取方法，其特征在于，包括以下步骤：步骤一：设定应用场景；

建立直角坐标系，其中z轴垂直向下；

给定规则网格重力垂直分量图，规则网格重力垂直分量图其x方向上等间距分布有M个网格节点，其y方向上等间距分布有N个网格节点，这样规则网格重力垂直分量图上有M×N个网格节点，M×N个网格节点对应M×N个重力垂直分量数据；

规则网格重力垂直分量图上M×N个网格节点对应的重力垂直分量数据记为gz(ξi,ηj)，i＝1,2,…,M，j＝1,2,…,N；(ξi,ηj)表示观测点所在的网格节点坐标，gz(ξi,ηj)表示网格节点坐标(ξi,ηj)处重力垂直分量值；M表示规则网格重力垂直分量图其x方向上的网格节点数，N表示规则网格重力垂直分量图其y方向上的网格节点数；

源点所在的网格节点坐标(xm,yn)处对应的重力场水平x分量数据记为gx(xm,yn)，源点所在的网格节点坐标(xm,yn)处对应的重力场水平y分量数据记为gy(xm,yn)，m＝1,2,…,M，n＝1,2,…,N；gx(xm,yn)和gy(xm,yn)即为待求解的重力场水平分量；

步骤二：根据公式(1)和公式(2)，分别计算加权系数；

式(1)和式(2)中，ωx(xm,yn；ξi,ηj)表示计算重力场水平x分量用的加权系数；ωy(xm,yn；ξi,ηj)表示计算重力场水平y分量用的加权系数；(xm,yn)和(ξi,ηj)都表示网格节点坐标，其中(ξi,ηj)表示观测点所在的网格节点坐标，(xm,yn)表示源点所在的网格节点坐标，当下角标m＝i时，xm与ξi两者数值相同，下角标不同时，两者数值不同；同样，当下角标n＝j时，yn与ηj两者数值相同，下角标不同时，两者数值不同；μpq，Xp，Yq，Rpq表示辅助计算变量，其计算式为：X1＝xm-ξi+0.5Δx，X2＝xm-ξi-0.5ΔxY1＝yn-ηj+0.5Δy，Y2＝yn-ηj-0.5Δyμpq＝(-1)p(-1)q

式中，p＝1,2；q＝1,2；Δx表示x方向上相邻网格节点间的间距，Δy表示y方向上相邻网格节点间的间距；

步骤三：根据步骤二计算得到的加权系数计算重力场水平分量；

2.根据权利要求1所述的重力场水平分量的获取方法，其特征在于，步骤三的实现方法为：(1)将步骤二计算得到的加权系数ωx(xm,yn；ξi,ηj)、ωy(xm,yn；ξi,ηj)分别排列成矩阵t的形式，其中矩阵t的形式如下：式(3)中，矩阵t中的各矩阵元素用ti,j表示，i＝1,2,…,M，j＝1,2,…,N；

矩阵t中的矩阵元素ti,j与对应的加权系数ω(x1,y1；ξi,ηj)存在如下关系：ti,j＝ω(x1,y1；ξi,ηj)            (4)加权系数ω(x1,y1；ξi,ηj)是一般表示，代表ωx(xm,yn；ξi,ηj)或者ωy(xm,yn；ξi,ηj)；

(2)将矩阵t补零扩展成矩阵

将矩阵 分成四个块矩阵，记为：

将块矩阵互换位置，得到矩阵cext：

(3)gz(ξi,ηj)，其中i＝1,2,…,M，j＝1,2,…,N；将gz(ξi,ηj)排列成矩阵g，矩阵g行数为M，列数为N；矩阵g中的矩阵元素gi,j与gz(ξi,ηj)存在如下关系：gi,j＝gz(ξi,ηj)        (8)将矩阵g补零扩展成行数为2M，列数为2N的矩阵gext；

式中， 表示Nx×Ny零矩阵，Nx×Ny即M×N；

(4)计算

式中，fft2()表示二维快速傅里叶变换；

(5)计算

式中，“.*”表示对应元素相乘运算；

(6)计算

式中，ifft2()表示二维快速傅里叶反变换；

(7)提取矩阵fext的前M行前N列，构成矩阵f，即为二维离散卷积计算结果，也就是重力场水平分量。

说明书

一种重力场水平分量的获取方法

技术领域

[0001] 本发明涉及一种应用于重力勘探领域的重力场分量获取方法，特别是一种重力场水平分量的获取方法。

背景技术

[0002] 重力勘探中常用的重力仪，测量的是重力场的垂直分量。重力垂直分量数据对于地质异常体的深部信息具有较强的分辨能力。受制于仪器水平，目前测量重力场水平分量的重力仪没有得到广泛使用。重力场水平分量对异常体的水平边界有很好的分辨能力。综合利用重力场三分量信息，能够提高重力勘探资料解释效果。

[0003] 根据经典位场理论，重力场水平分量与垂直分量之间存在定量关系，两者之间可以相互转换。根据实测的重力场垂直分量获取重力场水平分量，一般采用基于快速傅里叶变换算法(FFT)的波数域转换方法。波数域方法借助快速傅里叶变换算法的高效性，能够实现由大规模实测重力场垂直分量数据到重力场水平分量数据的快速计算。但是，由于波数域转换算子在零波数点存在奇异，导致波数域方法的计算精度低。不同于波数域方法，文献(Bingzhu Wang,Edward S.Krebes,Dhananjay Ravat.High-precision potential-field and gradient-component transformations and derivative computations using cubic B-splines.Geophysics,2008,73(5):I35-I42.)公开了一种空间域三次B样条函数重力场分量转换方法，该方法有效改善了重力场分量转换的计算精度，但是计算效率低。

[0004] 目前已有的获取重力场水平分量的方法存在的最大问题是不能同时保证计算效率和计算精度，无法满足大规模数据情况下快速、高精度重力场分量转换的要求。因此，寻找一种计算效率高、同时能保证计算精度的重力场水平分量的获取方法，具有重要的现实意义。

发明内容

[0005] 本发明针对现有技术存在的缺陷，提出了一种重力场水平分量的获取方法，其能够实现大规模重力测量数据情况下重力场水平分量的快速、高精度获取。

[0006] 为解决上述技术问题，本发明采用以下技术方案：

[0007] 一种重力场水平分量的获取方法，其步骤为：

[0008] 步骤一：设定应用场景；

[0009] 建立直角坐标系，其中z轴垂直向下。

[0010] 给定规则网格重力垂直分量图，规则网格重力垂直分量图其x方向上等间距分布有M个网格节点，其y方向上等间距分布有N个网格节点，这样规则网格重力垂直分量图上有M×N个网格节点，M×N个网格节点对应M×N个重力垂直分量数据。

[0011] 规则网格重力垂直分量图上M×N个网格节点对应的重力垂直分量数据记为gz(ξi,ηj)，i＝1,2,…,M，j＝1,2,…,N；(ξi,ηj)表示观测点所在的网格节点坐标，gz(ξi,ηj)表示网格节点坐标(ξi,ηj)处重力垂直分量值；M表示规则网格重力垂直分量图其x方向上的网格节点数，N表示规则网格重力垂直分量图其y方向上的网格节点数。

[0012] 源点所在的网格节点坐标(xm,yn)处对应的重力场水平x分量数据记为gx(xm,yn)，源点所在的网格节点坐标(xm,yn)处对应的重力场水平y分量数据记为gy(xm,yn)，m＝1,2,…,M，n＝1,2,…,N。gx(xm,yn)和gy(xm,yn)即为待求解的重力场水平分量。

[0013] 步骤二：根据公式(1)和公式(2)，分别计算加权系数；

[0014]

[0015]

[0016] 式(1)和式(2)中，ωx(xm,yn；ξi,ηj)表示计算重力场水平x分量用的加权系数；ωy(xm,yn；ξi,ηj)表示计算重力场水平y分量用的加权系数；(xm,yn)和(ξi,ηj)都表示网格节点坐标，其中(ξi,ηj)表示观测点所在的网格节点坐标，(xm,yn)表示源点所在的网格节点坐标，当下角标m＝i时，xm与ξi两者数值相同，下角标不同时，两者数值不同；同样，当下角标n＝j时，yn与ηj两者数值相同，下角标不同时，两者数值不同；μpq，Xp，Yq，Rpq(p＝1,2；q＝1,2)表示辅助计算变量，其计算式为

[0017] X1＝xm-ξi+0.5Δx，X2＝xm-ξi-0.5Δx

[0018] Y1＝yn-ηj+0.5Δy，Y2＝yn-ηj-0.5Δy

[0019]

[0020] 式中，Δx表示x方向上相邻网格节点间的间距，Δy表示y方向上相邻网格节点间的间距；

[0021] 步骤三：根据步骤二计算得到的加权系数计算重力场水平分量；

[0022] 根据公式(1)和公式(2)计算得到的加权系数，和一种二维离散卷积快速算法，实现公式(3)和公式(4)给出的重力场水平分量gx(xm,yn)和gy(xm,yn)的快速、高精度计算[0023]

[0024]

[0025] 与现有重力场垂直分量计算水平分量的转换方法相比，采用公式(1)和(2)给出的方法计算加权系数，能够提高转换结果的计算精度，并且保证分量转换的数值计算稳定性。

[0026] 步骤三中，根据公式(1)和公式(2)给出的加权系数和一种二维离散卷积快速算法，实现公式(3)和公式(4)给出的重力场水平分量快速、高精度计算，其步骤为：

[0027] (1)将步骤二计算得到的加权系数ωx(xm,yn；ξi,ηj)、ωy(xm,yn；ξi,ηj)均分别排列成矩阵t的形式，其中矩阵t的形式如下：

[0028]

[0029] 式(3)中，矩阵t中的各矩阵元素用ti,j表示，i＝1,2,…,M，j＝1,2,…,N。

[0030] 矩阵t中的矩阵元素ti,j与对应的加权系数ω(x1,y1；ξi,ηj)存在关系

[0031] ti,j＝ω(x1,y1；ξi,ηj)   (4)

[0032] 加权系数ω(x1,y1；ξi,ηj)是一般表示，代表ωx(xm,yn；ξi,ηj)或者ωy(xm,yn；ξi,ηj)。

[0033] (2)将矩阵t补零扩展成矩阵

[0034]

[0035] 将矩阵 分成四个块矩阵，记为

[0036]

[0037] 将块矩阵互换位置，得到矩阵cext

[0038]

[0039] (3)将gz(ξi,ηj)(i＝1,2,…,M，j＝1,2,…,N)排列成矩阵g，矩阵g中的矩阵元素gi,j与gz(ξi,ηj)存在关系

[0040] gi,j＝gz(ξi,ηj)   (8)

[0041] 将矩阵g补零扩展成矩阵gext

[0042]

[0043] 式中， 表示Nx×Ny零矩阵；

[0044] (4)计算

[0045] 式中，fft2()表示二维快速傅里叶变换；

[0046] (5)计算

[0047] 式中，“.*”表示对应元素相乘运算；

[0048] (6)计算

[0049] 式中，ifft2()表示二维快速傅里叶反变换；

[0050] (7)提取矩阵fext的前M行前N列，构成矩阵f，即为二维离散卷积计算结果，也就是重力场水平分量。

[0051] 与现有重力场垂直分量计算水平分量方法相比，采用上述二维离散卷积快速算法是本专利发明点之一，它保证了在大规模重力场垂直分量数据情况下水平分量计算的高效性。

[0052] 本发明是一个有机整体，根据步骤二提供的一种特殊的加权系数计算公式，采用步骤三运用的二维离散卷积快速算法，实现了大规模数据情况下重力场垂直分量计算水平分量在效率和精度上的统一。本发明解决了传统波数域方法因波数域算子在零波数存在奇异而导致水平分量计算精度低的问题，为满足大规模重力测量数据水平分量转换提供了方法支撑。

[0053] 与现有重力场分量转换方法相比，本发明具有以下优点：

[0054] (1)能够实现大规模重力测量数据情况下重力场水平分量快速、高精度获取；

[0055] (2)能够保证重力场分量转换数值计算的稳定性；

[0056] (3)大规模重力测量数据分量转换时，不但计算效率和计算精度高，并且所需计算机内存小。附图说明：

[0057] 图1为本发明的流程图；

[0058] 图2为重力观测数据网格节点坐标和重力水平分量网格节点坐标示意图；

[0059] 图3为棱柱体组合模型；

[0060] 图4为重力场垂直分量理论值；

[0061] 图5(a)为重力场水平x分量新方法计算结果；

[0062] 图5(b)为重力场水平x分量理论值；

[0063] 图6(a)为重力场水平y分量新方法计算结果；

[0064] 图6(b)为重力场水平y分量理论值。

[0065] 图中符号说明如下：

[0066] g/cm3表示克每立方厘米；mGal表示毫伽。具体实施方式：

[0067] 下面结合附图对本发明中的方法作进一步详细描述。

[0068] 参照图1，本发明提供一种重力场水平分量的获取方法，包括以下步骤：

[0069] 步骤一：设定应用场景

[0070] 参照图2，建立直角坐标系，z轴垂直向下为正；

[0071] 给定规则网格重力垂直分量图，规则网格重力垂直分量图其x方向上等间距分布有M个网格节点，其y方向上等间距分布有N个网格节点，这样规则网格重力垂直分量图上有M×N个网格节点，M×N个网格节点对应M×N个重力垂直分量数据。

[0072] 规则网格重力垂直分量图上M×N个网格节点对应的重力垂直分量数据记为gz(ξi,ηj)，i＝1,2,…,M，j＝1,2,…,N；(ξi,ηj)表示观测点所在的网格节点坐标，gz(ξi,ηj)表示网格节点坐标(ξi,ηj)处重力垂直分量值；M表示规则网格重力垂直分量图其x方向上的网格节点数，N表示规则网格重力垂直分量图其y方向上的网格节点数。

[0073] 源点所在的网格节点坐标(xm,yn)处对应的重力场水平x分量数据记为gx(xm,yn)，源点所在的网格节点坐标(xm,yn)处对应的重力场水平y分量数据记为gy(xm,yn)，m＝1,2,…,M，n＝1,2,…,N；gx(xm,yn)和gy(xm,yn)即为待求解的重力场水平分量。

[0074] 步骤二：根据公式(1)和公式(2)，分别计算加权系数。

[0075]

[0076]

[0077] 式(1)和式(2)中，ωx(xm,yn；ξi,ηj)表示计算重力场水平x分量用的加权系数；ωy(xm,yn；ξi,ηj)表示计算重力场水平y分量用的加权系数；(xm,yn)和(ξi,ηj)都表示网格节点坐标，其中(ξi,ηj)表示观测点所在的网格节点坐标，(xm,yn)表示源点所在的网格节点坐标，当下角标m＝i时，xm与ξi两者数值相同，下角标不同时，两者数值不同；同样，当下角标n＝j时，yn与ηj两者数值相同，下角标不同时，两者数值不同。μpq，Xp，Yq，Rpq表示辅助计算变量，其计算式为：

[0078] X1＝xm-ξi+0.5Δx，X2＝xm-ξi-0.5Δx

[0079] Y1＝yn-ηj+0.5Δy，Y2＝yn-ηj-0.5Δy

[0080]

[0081] 式中，p＝1,2；q＝1,2；Δx表示x方向上相邻网格节点间的间距，Δy表示y方向上相邻网格节点间的间距。

[0082] 步骤三：计算重力场水平分量。

[0083] 根据公式(1)和公式(2)给出的加权系数和一种二维离散卷积快速算法，实现公式(3)和公式(4)给出的重力场水平分量快速、高精度计算

[0084]

[0085]

[0086] 步骤三其实现步骤为：

[0087] (1)将加权系数ωx(xm,yn；ξi,ηj)、ωy(xm,yn；ξi,ηj)分别排列成矩阵t的形式，其中矩阵t的形式如下：

[0088]

[0089] 式(3)中，矩阵t中的各矩阵元素用ti,j表示；i＝1,2,…,M，j＝1,2,…,N。

[0090] 矩阵t中的矩阵元素ti,j与对应的加权系数ω(x1,y1；ξi,ηj)存在关系

[0091] ti,j＝ω(x1,y1；ξi,ηj)   (4)

[0092] 加权系数ω(x1,y1；ξi,ηj)是一般表示，代表ωx(xm,yn；ξi,ηj)或者ωy(xm,yn；ξi,ηj)。

[0093] (2)将矩阵t补零扩展成矩阵

[0094]

[0095] 将矩阵 分成四个块矩阵，记为

[0096]

[0097] 将块矩阵互换位置，得到矩阵cext

[0098]

[0099] (3)将重力场垂直分量数据gz(ξi,ηj)(i＝1,2,…,M，j＝1,2,…,N)排列成矩阵g，矩阵元素gi,j与重力场垂直分量数据存在关系

[0100] gi,j＝gz(ξi,ηj)   (8)

[0101] 将矩阵g补零扩展成矩阵gext

[0102]

[0103] 式中， 表示Nx×Ny零矩阵；

[0104] (4)计算

[0105] 式中，fft2()表示二维快速傅里叶变换；

[0106] (5)计算

[0107] 式中，“.*”表示对应元素相乘运算；

[0108] (6)计算

[0109] 式中，ifft2()表示二维快速傅里叶反变换；

[0110] (7)提取矩阵fext的前M行前N列，构成矩阵f，即为二维离散卷积计算结果，也就是重力场水平分量。

[0111] 下面对本发明方法的效果进行检验。

[0112] 为了说明本发明所提出的方法用于重力场垂直分量计算水平分量时的计算效率和计算精度，设计了如图3所示的棱柱体组合模型，模型参数由表1给出。观测平面范围为：x方向从-99800m到99800m，y方向从-99800m到99800m。模拟1:40000比例尺重力测量，网格间距取为Δx＝Δy＝400m，观测点数为500×500。

[0113] 将本发明提供的重力场水平分量的获取方法和传统基于快速傅里叶变换算法的波数域方法进行比较，其中本发明方法和传统波数域方法的程序代码均由Matlab语言编写，程序运行计算机为个人笔记本电脑，CPU主频为2.8GHz，内存为32GB。由表2可知，本发明提供的方法整个运行所需时间约为0.09秒，而传统波数域方法约为0.26秒，由此可见本发明提供的新方法计算效率高，且优于传统波数域方法。

[0114] 参见图4至图6，图4为重力场垂直分量理论值；图5(a)为重力场水平x分量新方法计算结果；图5(b)为重力场水平x分量理论值；图6(a)为重力场水平y分量新方法计算结果；图6(b)为重力场水平y分量理论值。对比可知，从形态上看，本发明提供的方法的计算结果与理论值吻合很好。以均方根误差衡量计算结果精度，由表2可知，本发明方法其重力场水平分量计算结果均方根误差为0.28nT，而传统波数域方法计算结果均方误差为2.98nT，可见本发明提供的方法计算精度高，且优于传统波数域方法。

[0115] 表1棱柱异常体模型参数

[0116] 3  x方向展布范围 y方向展布范围 z方向展布范围 密度(g/cm)
异常体1 -30km至-10km -10km至10km 0.1km至1.1km 1
异常体2 10km至30km -10km至10km 0.1km至1.2km 2
异常体3 -29km至29km -9.5km至9.5km 1.3km至2.2km 1

[0117] 表2计算时间和计算精度

[0118]  x分量计算精度(mGal) y分量计算精度(mGal) 计算时间(s)
波数域方法 0.31 0.54 0.26
空间域新方法 0.17 0.19 0.09

[0119] 以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。