任意密度分布复杂地质体重力场快速、高精度正演方法

摘要

本发明公开了任意密度分布复杂地质体重力场快速、高精度正演方法。本发明通过复杂地质体模型表示、棱柱体组合模型重力场计算(包括加权系数计算、二维离散卷积计算、重力场值合成)的步骤，实现了重力场正演计算在效率和精度上的统一。本发明解决了现有重力场正演方法不能同时保证计算效率和计算精度，无法满足大规模重力场三维密度反演、人机交互建模和解释的需求的问题。

说明书

技术领域：

本发明涉及一种重力场正演方法，特别是任意密度分布复杂地质体重力场快速、高精度 正演方法。

背景技术：

重力场正演是指根据密度分布计算重力场，反演是指根据观测重力值计算密度分布。正 演是反演的基础，正演计算的效率直接影响反演计算的效率，而正演计算精度直接影响反演 结果的质量。在重力勘探领域，伴随测绘技术和仪器的发展，重力测量在测量精度、空间分 辨率和测量范围上都有了显著提高，为重力反演提供了大规模高精度、高分辨率数据，重力 反演发展到三维密度反演阶段，成为国内外学者研究的热点。

随着计算机软硬件水平的不断提高，人机交互建模、解释也日益得到人们的重视，在地 球物理勘探中发挥着越来越重要作用。人机交互建模能够使人们通过直观的方式对地质体进 行建模，更容易融合地下结构的先验信息。反演方法与人机交互建模、解释方法相辅相成， 将极大提高人们对地球内部结构的认识。交互建模过程中，首先对地下结构进行剖分，根据 先验信息设计地质体分布。在勾勒出地质体分布后，进行正演计算，将正演结果与观测数据 进行比对，再对模型进行调整，如此反复，实现建模。

现阶段，由于重力观测数据的增加，正演计算量急剧增大，导致在一般计算机上难以实 现，成为制约三维密度反演发展的计算瓶颈。同时，人机交互建模对实时性有很高的要求， 而正演计算是交互建模中最大计算量所在，正演计算的效率直接影响人机交互建模的效果。

针对正演计算，众多国内外学者进行了研究。正演计算首先对地质体剖分，然后根据剖 分方式，采用某种方法计算重力场。文献(Zhdanov,M.S.,R.Ellis,S.Mukherjee. Three-dimensionalregularizedfocusinginversionofgravitygradienttensorcomponentdata. Geophysics,2004.69(4):925-937.)公开了一种结构化剖分方式，利用小棱柱体逼近复杂地质模 型，在计算过程中，采用等体积的球体重力场公式近似计算小棱柱体重力场，提高了正演计 算效率，但正演计算精度有所降低；文献(姚长利,郝天珧,管志宁,张聿文.重磁遗传算法 三维反演中高速计算及有效存储方法技术.地球物理学报,2003.46(2):252-258.)采用结构化 剖分方式，根据离散化后数学问题的特点，提出了“格架分离”技术和“格架等效计算方案”， 较好解决了计算效率和计算精度问题，但对于大规模剖分情形，该文献所给出的正演方法的 计算效率仍然比较低；文献(Tontini,F.C.,L.Cocchi,C.Carmisciano.Rapid3-Dforwardmodelofpotential fieldswithapplicationtothePalinuroSeamountmagneticanomaly(southernTyrrhenianSea,Italy).Journalof GeophysicalResearch,2009.114.)采用结构化剖分方法，采用三维傅里叶变换，给出了任意密度 分布情形下重力异常正演的波数域表达式，借助三维快速傅里叶变换算法，实现了快速正演 计算，该方法效率极高，但为克服截断效应，使用该方法前需要对剖分区域进行扩边，影响 了正演计算精度；此外，还有学者采用非结构化剖分方式，采用有限元方法计算重力场，这 种方法能够较精确刻画复杂地质体，计算精度高但计算效率很低。

剖分方式和计算方法共同决定了正演计算的效率和精度。正演计算的效率和精度是一对 矛盾体，目前已有的正演方法存在的最大问题是不能同时保证计算效率和精度，无法满足大 规模重力三维密度反演、人机交互建模和解释的需求。因此，寻找一种计算效率高、同时能 保证计算精度的正演计算方法具有重要的现实意义。

发明内容：

本发明针对目前现有重力场正演方法不能同时保证计算效率和计算精度，无法满足大规 模重力三维密度反演、人机交互建模和解释的需求问题，提出了任意密度分布复杂地质体重 力场快速、高精度正演方法。

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

任意密度分布复杂地质体重力场快速、高精度正演方法，其步骤为：

步骤一：复杂地质模型表示：建立包含所有目标区域的规则棱柱体模型，使得目标区域 (包含起伏地形)完全嵌入在该棱柱体模型中；将该棱柱体划分成许多小棱柱体，每个小棱 柱体密度为常值，不同棱柱体密度取值不同，以此刻画任意密度分布情形下复杂地质体；将 位于空气部分的小棱柱体的密度值设为零，以此刻画起伏地形；

步骤二：棱柱体组合模型重力场计算：步骤一中给出的棱柱体组合模型重力场其计算公 式为

$g_z(x_m,y_n,z_0)=\underset{r=1}{\overset{L}{\Sigma }}\underset{p=1}{\overset{M}{\Sigma }}\underset{q=1}{\overset{N}{\Sigma }}\rho ({\xi }_p,{\eta }_q,{\zeta }_r)h(x_m-{\xi }_p,y_n-{\eta }_q,z_0-{\zeta }_r)---\left(1\right)$

式(1)中，(x_m,y_n,z₀)表示观测点坐标，z₀为常值；L表示z方向棱柱体剖分个数；M表示 x方向棱柱体剖分个数；N表示y方向棱柱体剖分个数；(ξ_p,η_q,ζ_r)表示编号为(p,q,r)的小棱 柱体几何中心坐标；ρ(ξ_p,η_q,ζ_r)表示该棱柱体的密度值；h(x_m-ξ_p,y_n-η_q,z₀-ζ_r)表示加权系数；

实现上式的计算，分为三个环节：

首先，计算加权系数h(x_m-ξ_p,y_n-η_q,z₀-ζ_r)，其计算公式为

$h\left(x_m-{\xi }_p,y_n-{\eta }_q,z_0-{\zeta }_r\right)=-\gamma \underset{i=1}{\overset{2}{\Sigma }}\underset{j=1}{\overset{2}{\Sigma }}\underset{k=1}{\overset{2}{\Sigma }}{\mu }_{ijk}[z_k\arctan \frac{x_i{y}_j}{z_k{R}_{ijk}}-x_i\log (R_{ijk}+y_j)-y_j\log (R_{ijk}+x_i\left)\right]---\left(2\right)$

式(2)中，γ表示万有引力常数，Δx，Δy，Δz表示小棱柱体几何尺寸，arctan()表示反余 切函数运算符，log()表示自然对数运算符；其它符号含义如下

x₁＝ξ_p-0.5Δx-x_m，x₂＝ξ_p+0.5Δx-x_m，y₁＝η_q-0.5Δy-y_n，y₂＝η_q+0.5Δy-y_n，z₁＝ζ_r-0.5Δz-z₀， z₂＝ζ_r+0.5Δz-z₀，μ_ijk＝(-1)ⁱ(-1)^j(-1)^k，i＝1,2，j＝1,2，k＝1,2

其次，采用二维离散卷积快速计算方法来计算一层(相对z方向而言)棱柱体组合模型 重力场，其计算公式为

$g_z^r(x_m,y_n,z_0)=\underset{p=1}{\overset{M}{\Sigma }}\underset{q=1}{\overset{N}{\Sigma }}\rho ({\xi }_p,{\eta }_q,{\zeta }_r)h(x_m-{\xi }_p,y_n-{\eta }_q,z_0-{\zeta }_r)---\left(3\right)$

式(3)中，表示第r层(r＝1,2,…,L)棱柱体组合模型在高度面z₀产生的重力 场；(x_m,y_n,z₀)表示离散观测点坐标；

最后，将各层棱柱体组合模型重力场进行累加，得到整个组合模型的重力场， 即

$g_z(x_m,y_n,z_0)=\underset{r=1}{\overset{L}{\Sigma }}{g}_z^r(x_m,y_n,z_0)---\left(4\right)$

步骤二中所述的二维离散卷积快速计算方法，其步骤为：

(1)将加权系数h(x₁-ξ_p,y₁-η_q,z₀-ζ_r)排列成矩阵t，记为

$t=\left(\begin{array}{ccccc}{t}_{1-M,1-N}& \mathrm{...}& t_{1-M,0}& \mathrm{...}& t_{1-M,N-1}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ t_{0,1-N}& \mathrm{...}& t_{0,0}& \mathrm{....}& t_{0,N-1}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ t_{M-1,1-N}& \mathrm{...}& t_{M-1,0}& \mathrm{...}& t_{M-1,N-1}\end{array}\right)---\left(5\right)$

式(5)中，矩阵元素t_i,j与加权系数h(x₁-ξ_p,y₁-η_q,z₀-ζ_r)存在关系

t_i,j＝h(x₁-ξ_i+1,y_j+1-η_q,z₀-ζ_r)(6)

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

$\tilde{t}=\left(\begin{array}{cccccc}0& 0& \mathrm{...}& 0& \mathrm{...}& 0\\ 0& t_{1-M,1-N}& \mathrm{...}& t_{1-M,0}& \mathrm{...}& t_{1-M,N-1}\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ 0& t_{0,1-N}& \mathrm{...}& t_{0,0}& \mathrm{...}& t_{0,N-1}\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ 0& t_{M-1,1-N}& \mathrm{...}& t_{M-1,0}& \mathrm{...}& t_{M-1,N-1}\end{array}\right)---\left(7\right)$

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

$\tilde{t}=\left(\begin{array}{cc}{\tilde{t}}_{11}& {\tilde{t}}_{12}\\ {\tilde{t}}_{21}& {\tilde{t}}_{22}\end{array}\right)---\left(8\right)$

将块矩阵互换位置，得到矩阵c^ext

$c^{ext}=\left(\begin{array}{cc}{\tilde{t}}_{22}& {\tilde{t}}_{21}\\ {\tilde{t}}_{12}& {\tilde{t}}_{11}\end{array}\right)---\left(9\right)$

(3)将第r层密度值ρ(ξ_p,η_q,ζ_r)(p＝1,2,…,M，q＝1,2,…,N)排列成矩阵g，矩阵元素 g_i,j与密度值存在关系

g_i,j＝ρ(ξ_i,η_j,ζ_r)(10)

将矩阵g补零扩展成矩阵g^ext

$g^{ext}=\left(\begin{array}{cc}g& 0_{M\times N}\\ 0_{M\times N}& 0_{M\times N}\end{array}\right)---\left(11\right)$

式(11)中，0_M×N表示M×N零矩阵；

(4)计算${\hat{c}}^{ext}:=fft2\left(c^{ext}\right),{\hat{g}}^{ext}:=fft2\left(g^{ext}\right)$

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

(5)计算${\hat{f}}^{ext}:={\hat{c}}^{ext}.*{\hat{g}}^{ext}$

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

(6)计算$f^{ext}:=ifft2\left({\hat{f}}^{ext}\right)$

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

(7)提取矩阵f^ext的前M行前N列，构成矩阵f，即为二维离散卷积计算结果。

本发明是一个有机整体，即在特定的模型表示方式条件下，建立棱柱体重力场叠加模型， 根据一种特殊的加权系数计算公式，采用二维离散卷积快速计算方法，实现了重力场正演计 算在效率和精度上的统一。

与现有技术相比，本发明具有以下优点：

(1)模型表示方法简单、灵活，很容易刻画任意密度分布复杂地质体以及起伏地形；

(2)能够实现任意密度分布情况下复杂地质体重力场的快速、高精度计算，可以满足大 规模重力三维密度反演、人机交互建模和解释的需求；

(3)大规模正演计算时，算法不但计算效率和计算精度高，并且所需计算机内存小。

附图说明：

1、图1为重力场快速、高精度正演流程图；

2、图2为复杂地质模型表示；

3、图3为组合模型剖面图；

4、图4为组合模型重力场正演计算值；

5、图5为组合模型重力场理论值；

6、图6重力场理论值与计算值的差值；

图中符号说明如下：

L：表示z方向剖分小棱柱体个数；

M：表示x方向剖分小棱柱体个数；

N：表示y方向剖分小棱柱体个数；

ρ：表示密度；

具体实施方式：

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

1、复杂地质模型表示：

首先，建立包含所有目标区域的规则棱柱体模型，确定棱柱体在x，y，z方向的起始位 置，使得目标区域(包含起伏地形)完全嵌入在该棱柱体模型中；

其次，根据实际问题需求，将棱柱体划分成许多规则小棱柱体(如图2所示)，确定小棱 柱体的几何尺寸Δx，Δy，Δz；

最后，根据目标区域的密度分布，对每个小棱柱体密度进行赋值，位于空气部分的小棱 柱体，其密度值设为零；

2、棱柱体组合模型重力场计算：

步骤一中给出的棱柱体组合模型，其重力场计算公式为

$g_z(x_m,y_n,z_0)=\underset{r=1}{\overset{L}{\Sigma }}\underset{p=1}{\overset{M}{\Sigma }}\underset{q=1}{\overset{N}{\Sigma }}\rho ({\xi }_p,{\eta }_q,{\zeta }_r)h(x_m-{\xi }_p,y_n-{\eta }_q,z_0-{\zeta }_r)---\left(12\right)$

式(12)中，(x_m,y_n,z₀)表示观测点坐标，z₀为常值；L表示z方向棱柱体剖分个数；M表 示x方向棱柱体剖分个数；N表示y方向棱柱体剖分个数；(ξ_p,η_q,ζ_r)表示编号为(p,q,r)的小 棱柱体几何中心坐标；ρ(ξ_p,η_q,ζ_r)表示该棱柱体的密度值；h(x_m-ξ_p,y_n-η_q,z₀-ζ_r)表示加权系 数；

实现式(12)的计算，分为三个环节：

首先，根据观测点坐标(x_m,y_n,z₀)和小棱柱体几何中心坐标(ξ_p,η_q,ζ_r)，计算加权系数 h(x_m-ξ_p,y_n-η_q,z₀-ζ_r)，其计算公式为

$h\left(x_m-{\xi }_p,y_n-{\eta }_q,z_0-{\zeta }_r\right)=-\gamma \underset{i=1}{\overset{2}{\Sigma }}\underset{j=1}{\overset{2}{\Sigma }}\underset{k=1}{\overset{2}{\Sigma }}{\mu }_{ijk}[z_k\arctan \frac{x_i{y}_j}{z_k{R}_{ijk}}-x_i\log (R_{ijk}+y_j)-y_j\log (R_{ijk}+x_i\left)\right]---\left(13\right)$

式(13)中，γ表示万有引力常数，Δx，Δy，Δz表示小棱柱体几何尺寸，arctan()表示反余 切函数运算符，log()表示自然对数运算符；其它符号含义如下

x₁＝ξ_p-0.5Δx-x_m，x₂＝ξ_p+0.5Δx-x_m，y₁＝η_q-0.5Δy-y_n，y₂＝η_q+0.5Δy-y_n，z₁＝ζ_r-0.5Δz-z₀， z₂＝ζ_r+0.5Δz-z₀，μ_ijk＝(-1)ⁱ(-1)^j(-1)^k，i＝1,2，j＝1,2，k＝1,2

其次，采用二维离散卷积快速计算方法来计算一层(相对z方向而言)棱柱体组合模型 重力场，其计算公式为

$g_z^r(x_m,y_n,z_0)=\underset{p=1}{\overset{M}{\Sigma }}\underset{q=1}{\overset{N}{\Sigma }}\rho ({\xi }_p,{\eta }_q,{\zeta }_r)h(x_m-{\xi }_p,y_n-{\eta }_q,z_0-{\zeta }_r)---\left(14\right)$

式(14)中，表示第r层(r＝1,2,…,L)棱柱体组合模型在高度面z₀产生的重 力场；(x_m,y_n,z₀)表示离散观测点坐标；

最后，将各层棱柱体组合模型重力场(r＝1,2,…,L)进行累加，得到整个组 合模型的重力场，即

$g_z(x_m,y_n,z_0)=\underset{r=1}{\overset{L}{\Sigma }}{g}_z^r(x_m,y_n,z_0)---\left(15\right)$

步骤二中所述的二维离散卷积快速计算方法，其步骤为：

(1)将加权系数h(x₁-ξ_p,y₁-η_q,z₀-ζ_r)排列成矩阵t，记为

$t=\left(\begin{array}{ccccc}{t}_{1-M,1-N}& \mathrm{...}& t_{1-M,0}& \mathrm{...}& t_{1-M,N-1}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ t_{0,1-N}& \mathrm{...}& t_{0,0}& \mathrm{....}& t_{0,N-1}\\ .& .& .& .& .\\ .& .& .& .& .\\ .& .& .& .& .\\ t_{M-1,1-N}& \mathrm{...}& t_{M-1,0}& \mathrm{...}& t_{M-1,N-1}\end{array}\right)---\left(16\right)$

式(16)中，矩阵元素t_i,j与加权系数h(x₁-ξ_p,y₁-η_q,z₀-ζ_r)存在关系

t_i,j＝h(x₁-ξ_i+1,y_j+1-η_q,z₀-ζ_r)(17)

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

$\tilde{t}=\left(\begin{array}{cccccc}0& 0& \mathrm{...}& 0& \mathrm{...}& 0\\ 0& t_{1-M,1-N}& \mathrm{...}& t_{1-M,0}& \mathrm{...}& t_{1-M,N-1}\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ 0& t_{0,1-N}& \mathrm{...}& t_{0,0}& \mathrm{...}& t_{0,N-1}\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ .& .& .& .& .& .\\ 0& t_{M-1,1-N}& \mathrm{...}& t_{M-1,0}& \mathrm{...}& t_{M-1,N-1}\end{array}\right)---\left(18\right)$

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

$\tilde{t}=\left(\begin{array}{cc}{\tilde{t}}_{11}& {\tilde{t}}_{12}\\ {\tilde{t}}_{21}& {\tilde{t}}_{22}\end{array}\right)---\left(19\right)$

将块矩阵互换位置，得到矩阵c^ext

$c^{ext}=\left(\begin{array}{cc}{\tilde{t}}_{22}& {\tilde{t}}_{21}\\ {\tilde{t}}_{12}& {\tilde{t}}_{11}\end{array}\right)---\left(20\right)$

(3)将第r层密度值ρ(ξ_p,η_q,ζ_r)(p＝1,2,…,M，q＝1,2,…,N)排列成矩阵g，矩阵元素 g_i,j与密度值存在关系

g_i,j＝ρ(ξ_i,η_j,ζ_r)(21)

将矩阵g补零扩展成矩阵g^ext

$g^{ext}=\left(\begin{array}{cc}g& 0_{M\times N}\\ 0_{M\times N}& 0_{M\times N}\end{array}\right)---\left(22\right)$

式(22)中，0_M×N表示M×N零矩阵；

(4)计算${\hat{c}}^{ext}:=fft2\left(c^{ext}\right),{\hat{g}}^{ext}:=fft2\left(g^{ext}\right)$

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

(5)计算${\hat{f}}^{ext}:={\hat{c}}^{ext}.*{\hat{g}}^{ext}$

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

(6)计算$f^{ext}:=ifft2\left({\hat{f}}^{ext}\right)$

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

(7)提取矩阵f^ext的前M行前N列，构成矩阵f，即为二维离散卷积计算结果。

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

为了说明本发明所提出的方法用于计算任意密度分布情况下复杂地质构造重力场时的效 率和精度，设计了如下组合模型(图3所示)：

密度均匀的棱柱体内嵌一个密度均匀的球体，球心与棱柱体中心重合。棱柱体范围为：x 方向从-10000m到10000m，y方向从-10000m到10000m，z方向从0m到3000m(z轴向下为正)； 球体半径为1000m。棱柱体密度为1g/cm³，球体的密度为5g/cm³。将棱柱体剖分成1000× 1000×500个大小相同的小棱柱体，计算高度为-200m平面(图3中虚线所示)上的重力场，计 算点个数为1000×1000。

正演算法利用Fortran语言编程实现，运行程序所用的个人台式机配置为：CPU为i7-2620， 主频为2.7GHz，内存为32GB，四核八线程。运行所需时间约为60秒，由此可见正演算法效率 很高。组合模型重力场正演算法计算值和理论值分别如图4和图5所示，从形态上看，两者是 一致的。理论值减去计算值得到差值(图6所示)，对差值进行统计，统计结果由表1给出，可 知正演算法精度很高。

表1组合模型重力场理论值和正演计算误差统计量(单位：mGal)

最大值 最小值 均值 均方值 理论值 145.54444 29.475854 91.824949 16.393163 误差 0.001254 -0.000299 0.000107 0.000219

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