基于非结构化网格的重力梯度张量反演

2020-02-27黄天统彭新发朱自强

物探与化探 2020年1期

黄天统，彭新发，朱自强

(核工业二三〇研究所，湖南长沙 410007)

0 引言

作为一种高精度的地球物理勘探方法，在实际勘探过程中重力梯度张量测量的是重力位的二阶导数。与传统的重力勘探相比，重力梯度张量测量对地下介质密度变化反映更加灵敏，具有更好的抗干扰能力，能够更加直接地反映地下密度异常体的水平边界位置[1-3]。在国外，重力梯度张量勘探技术已经较为广泛地应用于海洋领域、航空领域以及卫星重力领域。

结构化网格是指网格区域内所有的内部点都具有相同的毗邻单元。由于非结构网格具有生成速度快、网格的数据结构简单、网格的质量好等优点[4]，广泛地应用于形状规则的物体剖分。但是随着计算机技术的飞速发展，结构化网格的弊端越来越明显，比如：待剖分物体的几何形状越来越复杂，结构化网格剖分带来了的剖分误差也越来越大，并且节点和网格单元在高质量情况下不能合理地呈梯度状分布在剖分的网格中[5]。

与结构化网格相比，非结构化网格弥补了其不能解决所有形状和任意连通区域网格剖分的不足；对于网格单元和节点的编号，非结构化网格并没有固定的规则，两个相邻单元格编号对应的网格单元不一定相邻，这样使得非结构化网格能够十分方便的控制单元和节点的分布，具有很强的适应性，能够更好地剖分复杂模型。

在地球物理数值模拟过程中，高效地拟合场源模型及地形分布将更加有利于模型的计算，非结构化网格与生俱来的优势使其广泛地应用于地球物理数值模拟中。汤井田[6]、任政勇[7]在基于非结构化网格的情况下实现了2.5维的自适应有限元直流电阻率的正演模拟；孙丽影[8]、刘长生[9]、童孝忠[10]通过非结构化网格进行局部加密完成了大地电磁二维、三维有限元数据模拟；赵晓博利用Delaunay三角化实现了中心回线瞬变电磁2.5维的正演模拟；曾思红[11]、曹书锦[12]利用有限单元法在非结构化网格下对重力梯度张量进行了二维、三维正演，Jahandari[13]则是运用有限体积法实现了基于非结构化网格的二维、三维的重力高精度正演模拟，取得了较好的效果；Okabe[14]推导了任意多边形和多面体网格下的重磁一阶、二阶导数的二维、三维解析表达式，具有很高的计算精度。

目前，对于重力梯度张量数据反演，大多都是基于长方形或者长方体等结构化网格进行的，但由于结构化网格的几何适应能力差，在较复杂地质体或地形的情况下会带来较大的剖分误差，从而使得反演效果有所下降。与结构化网格相比，非结构化网格能够减少模型的剖分误差，更好地拟合地下密度异常体的分布和地表地形的变化，因此，研究基于非结构化网格的重力梯度张量反演算法，可以提高重力梯度张量反演过程中的计算精度。现在，根据非结构化网格进行的地球物理反演主要还是基于二维三角形网格，elièvre[15]、吉日嘎拉图[16]利用地震初至波完成了二维慢度的反演，并完成了非结构化网格重震二维联合反演，使得反演的纵横向分辨率都很高，降低了反演结果的多解性问题。对于位场三维非结构化网格的反演，国内外文献中很少有涉及到。

1 正反演算法

1.1 基于非结构化网格重力梯度张量正演算法

重力梯度张量是重力位的二阶导数，在笛卡尔坐标系下，根据牛顿万有引力定律，一个剩余密度为ρ，观测点为P(p0,q0,r0)，位于Q(p,q,r)地下异常体ν在地下产生的重力位V有：

(1)

其中：G为万有引力常量，u=-(x2+y2+z2)-1/2,x=p-p0、y=q-q0、z=r-r0。则重力梯度张量表达式写成矩阵形式：

(2)

其中：Γ表示全重力梯度张量。重力梯度张量共有9个分量，由于在无源空间里重力的旋度和散度都为0，有Vxx+Vyy+Vzz=0，Vyx=Vxy，Vzx=Vxz，Vzy=Vyz，因此重力梯度张量中只有5个独立分量Vxx、Vxy、Vxz、Vyz、Vzz。

重力位在k方向上的一阶导数：

(3)

其中：k是k方向法向量：kT=(kx,ky,kz)或者kT=(cos (x,k),cos (y,k),cos (z,k))。

将式(1)带入式(3)得：

(4)

引入另一个方向向量l：lT=(lx,ly,lz)，对于不同的重力梯度张量分量，k、l的方向余弦取值如表1所示。对Vk求在l方向求导，得到重力位的二阶导数：

(5)

表1 方向余弦的取值与重力梯度张量各分量的关系

对式(5)利用散度定理得：

(6)

其中：S为地下异常体ν的表面，n为面S的法向量，指向异常体外部。

在四面体网格的形式下,式(6)的积分可以写成每个网格的积分总和的形式：

(7)

其中：

(8)

为了快速地求解出面Si上Ki的值，采用了坐标选择的方法。旋转分为两步：首先，以z轴为中心逆时针[15]旋转xOy平面(笛卡尔坐标系)，直到x方向与面的外法向量在xOy平面上的投影平行，旋转的角度为θ，得到最终的坐标系统(X,Y,Z)。然后，以Y轴为中心逆时针选择x′Oz轴，直到z方向与面的外法向量平行，旋转的角度为φ，得到最终的坐标系统(X,Y,Z)，如图1所示。

图1 面坐标系统旋转示意

坐标的旋转转换可以写成：

(9)

其中：0≤θ<2π,0≤φ≤π。通过旋转，面Si上的3个顶点的Z值相等，面Si的外法向量ni可以写成：

(10)

在(X,Y,Z)坐标系统下，式(8)可以写成：

(11)

其中：

(12)

通过对Si面上的每条边进行二维坐标旋转，可以将式(11)的面积分转换为线积分。以原点为中心，逆时针旋转X轴、Y轴，直到X方向与Si面上的j边重合，得到新的坐标系统(ξ,η)旋转的角度为ψ，0≤ψ<2π，如图2所示。坐标的旋转转换可以写成：

(13)

通过旋转，Si面上j边的两个端点的η值相等，j边的法向量sj为：

图2 边坐标系统旋转示意

(14)

那么，式(11)就可以写成：

(15)

从而得到：

(16)

对式(16)求积分得到：

ln[ξ+(ξ2+η2+Z2)1/2]·

(17)

最终得到重力位二阶导数的表达式:

(18)

如果Z或者cosψ的值为0，则式(17)的后面一项为0，但是此时式(17)的第一项的值不一定一直存在，例如η和Z的值都趋近于0，且ξ不为正数。当ξj和ξj+1都为负时，通过式(19)可以避免，对于其他情形可用一个较小的值代替。

(19)

1.2 反演的基本理论

由于重力梯度张量反问题是一个不适定问题，为了得到一个稳定的可行解，李耀国等人[16]在吉洪诺夫正则化理论的基础上提出了一个能够广泛应用于地球物理反问题的广义目标函数，其中数据拟合函数为:

其中：Wd=diag{1/ε1,…,1/εN}，εi为第i个观测数据的标准差。通常假设实际数据的误差是相互独立且满足高斯分布的，实测得到重力梯度张量数据dobs=F(m)+δ，δ表示的是观测点上的随机噪声。由于随机噪声δ是满足高斯分布N(0,ε2)的，ε为观测数据的标准差。根据最小二乘拟合，那么反演最终的数据拟合差为：

(21)

其中：N为测点数。

在最小化目标函数的过程中，在目标函数中引入的模型目标函数为：

(22)

φm(m)=αs‖Ws(m-mref)‖2+

αx‖Wx(m-mref)‖2+αy‖Wy(m-mref)‖2+

αz‖Wz(m-mref)‖2

(23)

其中：αs、αx、αy、αz为每项的加权系数，Ws为粗糙度矩阵，Wx、Wy、Wz分别为x、y、z三个方向上的差分算子。令：

(24)

然而，在非结构化网格下，差分矩阵Wx、Wy、Wz的求解不能像在结构化网格情况下那样简单的求解。为了解决差分矩阵的计算问题，在Lelièvre[15]提出的二维求解Wm的计算方法基础上进行衍生，得到三维情况下两个相邻网格Wm的计算方法。

(25)

其中：V1、V2为两个相邻四面体网格的体积，R为两个四面体网格的中心距离。

结合数据拟合函数式(1)和模型目标函数式(4)可以得到广义的反演目标函数：

1.3 深度加权函数

为了使得反演结果能够反演出深度上的结构，提高深度方向上的分辨率，李耀国等人[19]引入了一个深度加权函数来减小核矩阵在深度方向上的衰减：

(27)

其中：z是网格单元的中心埋深，z0是个常数，β是一个接近于3的常数，在已知网格剖分方式和观测面的高度后，通过调节z0和β的值，来尽可能地减少核函数在深度方向上的衰减，从而使得反演结果能够准确地分布在对应的深度上。式(24)结合深度加权函数得到

(28)

将深度加权函数结合式(22)得到带深度加权的模型目标函数：

φm(m)=‖Wm(m-mref)‖2

(29)

1.4 多分量联合反演

多分量的联合反演的目标函数只需要对式(26)中的数据拟合项进行适当的修改，即将各个分量的数据拟合项结合起来形成一个总的数据拟合项。结合模型目标函数，对于重力梯度张量多分量联合反演的目标函数可以写成：

a‖Wm(m-mref)‖2，

(30)

α,β=x,y,z

2 算法验算

图3所示为单个长方体模型，场源空间大小为1 000 m×1 000 m×500 m，其内埋藏着一个大小400 m×400 m×200 m的存在密度异常的长方体，长方体顶部埋深为150 m，底部埋深为350 m，水平中心位置在地表的投影坐标为(500 m,500 m)，长方体剩余密度为1.0 g/cm3，观测位置位于地表，x、y方向测线的长度1 000 m、点距50 m。根据已剖分的网格单元，应用前文推导的公式计算出来的重力梯度张量各分量的正演结果如图4所示。

图3 单个长方体模型非结构化网格剖分示意

图4 单个长方体模型重力梯度张量正演结果

从图4可以看出，重力梯度张量各个分量的正演结果都能够明显地反映地下密度异常体的边界位置。Vxy的4个极值点对应了长方体4个顶点的水平投影位置；Vxx、Vyy的极小值位置对应的是长方体水平中心位置，并且等值线变化最密集的部分正好是长方体的4条边界位置；Vxz、Vyz的极值位置对应的是长方体边界的水平中心位置；Vzz的等值线图类似一系列的同心矩形，等值线在地下异常体密度变化边界分布密集，能够很好地确定密度变化，极大值的中心位置对应的是长方体的中心位置。

将单个长方体模型计算出来的重力梯度张量正演结果，利用正则化方法并结合深度加权函数对重力梯度张量各分量进行反演，其中初始参考模型mref=0，深度加权矩阵中的β取值为 3。各分量反演结果如图5所示。

从图5可以看出，反演结果对异常顶部位置的分辨率比底部分辨率要高，这是由于位场本身的特征决定的。对比各个分量之间反演效果，各分量反演出来的剩余密度值与真实值比较接近，但是都出现了一定的负值，其中，Vxz、Vyy、Vyz反演出来的负值较大；各分量的反演结果在深度方向上延伸较大；Vzz反演出的剩余密度更接近真实值，反演的剩余密度极大值位置与真实位置非常吻合,Vxx、Vxy、Vxz、Vyy、Vyz、Vzz反演的密度值与理论模型密度的均方误差分别是0.214 8、0.205 3、0.201 5、0.217 2、0.213 1、0.192 7，说明了对于单分量的反演，基于非结构化网格的重力梯度张量反演方法是十分有效的，算法是可行的。

图5 单个长方体模型单分量反演结果

为了充分利用重力梯度张量各分量包含的密度、空间信息，对重力梯度张量的独立分量进行多分量的联合反演，结果如图6所示。对比图3与图6可知，多分量联合反演出来的密度值与理论模型的密度的均方误差为0.190 9，比任一单分量反演的均方误差都小，且利用多分量进行联合反演，其结果与理论模型极其吻合，这也证明了多分量联合反演的正确性。与单分量反演相比，联合反演出来的剩余密度值更加接近真实值，反演结果在水平方向和深度方向上都更加集中收敛于异常体的中心位置，利用多个独立分量进行联合反演能够更加准确反演出地下介质的剩余密度及埋藏位置，再次验证了算法的正确性。

图6 单个长方体模型多分量联合反演结果

3 Y型岩脉模型

场源空间大小为1 000 m×1 000 m×500 m，如图7所示，场源空间存在一个Y型岩脉，其剩余密度都为1.0 g/cm3，左侧板状体顶部埋深为100 m，顶部中心在地表的投影位置为(200 m,500 m)，顶部尺寸为200 m×400 m，底部埋深为300 m，底部中心在地表的投影位置为(300 m,500 m)，顶部尺寸为200 m×400 m。右侧板状体顶部埋深为100 m，顶部中心在地表的投影位置为(700 m,500 m)，顶部尺寸为200 m×400 m，底部埋深为400 m，底部中心在地表的投影位置为(500 m,500 m)，顶部尺寸为200 m×400 m。x、y方向点距为50 m，测线长度为1 000 m。利用非结构化网格剖分软件将场源空间剖分，网格数目为9 408。

加入5%的高斯噪声后，利用重力梯度张量的5个相互独立分量进行多分量的联合反演，多分量联合反演的密度值与理论模型密度的均方误差为0.2347，反演结果如图8所示。

对比图7与图8可知，利用多分量联合反演的结果与理论模型相似，反演的密度异常体能够很好地确定异常的位置、倾向。

4 结论

1)笔者利用吉洪诺夫正则化方法结合李耀国提出的广义目标函数将重力梯度张量各个分量正演结果分别进行了单分量反演；其次，为了充分利用重力梯度张量各分量数据使得反演结果更加准确，将5个相互独立的重力梯度张量分量进行了联合反演；最后，用长方体模型验证了算法的正确性，与长方体网格下的反演结果作对比说明非结构化网格下反演的优势，并通过Y型岩脉模型来说明方法的适应性。得出以下结论： 1) 对重力梯度张量各个分量的正演结果进行吉洪诺夫正则化反演，反演过程引入了深度加权函数，反演结果能够较好地确定地下异常体的水平、深度位置，反演剩余密度值也与理论值比较接近，验证了方法的正确性。