一种构建复杂形体引力梯度场的计算方法与流程

本发明属于重力及重力梯度测量技术领域，尤其是一种构建复杂形体引力梯度场的计算方法。

背景技术：

重力梯度仪自梯度、实验室基础梯度以及实验装置所产生的人工梯度激励的正演计算，需要计算物体周围引力梯度场。正演问题主要根据给定质量体的形状、产状、空间位置和物理特性等，通过理论或者数值计算来求得它在观测坐标系下产生的异常大小、特征和变化规律。只有求出不同形体的引力场分布，并总结出场的特征与几何参数以及物性之间互相联系的内在规律，才能运用这些规律对重力梯度异常做解释推断，因此，正演也是重力梯度反演及地质解释的基础。正演的基本原理如下：

地球重力的主体是万有引力，在空间r′处的密度体ρ(r′)在测点r处的重力位为：

其中，g是万有引力常数。重力场g是重力位函数u(r)的空间梯度，即重力位在直角坐标系三个方向上的一阶导数，

gx,gy,gz分别表示重力在x,y,z方向上的分量。重力梯度γ是重力位函数u(r)的二阶导数，其表示如下：

其中，ui,j(i,j＝x,y,z)为重力梯度张量各分量，其物理意义为gi在j方向上的空间变化率，单位为e(厄缶)，1e＝10^-9·1/s²，即在相距为1m的位置上重力的变化为10^-9m/s²。由于重力梯度张量矩阵是一个对称矩阵，且主对角元素之和为0，所以上式的9个分量中只有5个是独立的。

形状简单规则的质量体引起的重力梯度异常存在明确的理论解，如球体、棱柱体。forsberg给出了直角坐标系下单一矩形棱柱体在原点处六个重力梯度分量异常的解析公式，其中

式中，xi＝ξi-x，yj＝ηj-y，zk＝ζk-z，μijk＝(-1)ⁱ(-1)^j(-1)^k。

对于形状复杂质量体引力梯度的计算，主要采用对复杂形体进行分割，使之转化为一系列存在理论解的简单形体的组合，计算各简单形体的重力梯度再进行张量矩阵叠加，即得复杂形体的重梯分布。目前将复杂形体分割为简单规则体仍依靠手动完成，操作过程复杂且很难做到细致分割，导致分割后与实际形体拟合度较差，计算精度受到限制。因此只能实现复杂程度不高的形体周围引力梯度计算，效率低下，工程应用局限性大。

技术实现要素：

本发明的目的在于克服现有技术的不足，提供一种构建复杂形体引力梯度场的计算方法，解决目前手动将复杂形体有限分割成简单规则体所产生的实现复杂、效率地下且工程应用局限性大的问题。

本发明解决现有的技术问题是采取以下技术方案实现的：

一种构建复杂形体引力梯度场的计算方法，包括以下步骤：

步骤1、将复杂形体等效为带电绝缘体，将引力场等效为电场，对模型进行网格划分；

步骤2、根据库伦定律，计算两个物体之间的引力场，将密度等效为带电体的电荷密度并乘以转换系数，求得复杂形体周围的引力场；

步骤3、在直角坐标系中对各引力分量沿三个坐标轴进行求导，将引力各分量视为物体沿三个坐标方向的变形，求出在该变形下所产生的应变，其中线应变对应直线引力梯度，切应变对应交叉引力梯度。

所述步骤1对模型进行网格划分的具体方法为：对于规整的模型，将其划分为全六面体单元，对于形状复杂及含较多曲面的模型，将其划分为四面体单元，或者四面体与六面体单元的结合。

所述步骤2计算两个物体之间的引力场f的公式为：

式中，r为两者之间的距离，q1,q2为两个点电荷的电荷量，er为从q1到q2方向的矢径，k＝9×10⁹nm²/c²为库伦常数。

所述六面体单元在直角坐标系中受力所产生的变形及其应变之间的关系满足如下关系：

式中，u,v,w为微元体受力后分别沿x,y,z三个坐标方向的变形；εx,εy,εz分别为x,y,z三个方向的线应变，即微元体变形后沿该方向单位长度的改变量；γxy,γyz,γxz为绕三个坐标轴的切应变，即微元体两条相互垂直的棱边在变形后的直角改变量。

本发明的优点和积极效果是：

1、本发明利用万有引力与库伦力表达式的相似性，将质量体等效为带电绝缘体，将引力场等效为电场，实现了任意三维质量体周围引力梯度的计算。其充分利用电场分析模块的相关功能，方便地实现任意三维形体的自动精确划分，扩大了有限元法在的应用范围，提高模型拟合度高，计算误差小于2％。该方法可以大幅度提高计算效率，同时具有丰富的图像显示能力，可以方便地显示、考察任意观测面的仿真结果，便于后续分析。

2、本发明设计合理，适用于重力梯度仪自梯度、实验室基础梯度以及实验装置所产生的人工梯度激励的正演计算，为重力梯度反演及地质解释打下坚实基础。

附图说明

图1为直角坐标系下的矩形棱柱体；

图2为带电绝缘体周围电场分布示意图；

图3a为各阶梯度张量分量示意图(uxx分量)；

图3b为各阶梯度张量分量示意图(uxy分量)；

图3c为各阶梯度张量分量示意图(uxz分量)；

图3d为各阶梯度张量分量示意图(uyy分量)；

图3e为各阶梯度张量分量示意图(uyz分量)；

图3f为各阶梯度张量分量示意图(uzz分量)；

图4为相对误差云图。

具体实施方式

以下结合附图对本发明实施例做进一步详述：

一种构建复杂形体引力梯度场的计算方法，是基于以下原理实现的：

真空中两个静止的点电荷之间的相互作用力满足库伦定律，即与其电荷量的乘积成正比，与其距离的二次方成反比，作用力的方向在两者的连线上，同名电荷相斥，异名相吸，物体之间的引力场表达式为：

式中，r为两者之间的距离，q1,q2为两个点电荷的电荷量，er为从q1到q2方向的矢径，k＝9×10⁹nm²/c²为库伦常数。该公式与任意两个质点间的万有引力表达式形式相同。根据类比的思想，可以采用公式(5)去计算物体之间的引力场。进一步分析可知，一个带正电荷的绝缘体周围的电场分布与同样形状和体积的质量体周围的引力场分布是完全相同的，只需将质量体的密度等效为带电体的电荷密度，再乘以相应的转换系数。至此，可以成功求得复杂形体周围的引力场。

为了求得引力梯度场，还需在直角坐标系中对各引力分量沿三个坐标轴进行求导，然而计算公式(5)没有对电场力求方向导数的功能，这是仿真过程中需要解决的另一难题。联想到弹性力学的相关知识，正六面体弹性微元在直角坐标系中受力所产生的变形及其应变之间的关系满足如下关系。

式中，u,v,w为微元体受力后沿三个坐标方向的变形；εx,εy,εz分别为x,y,z三个方向的线应变，即微元体变形后沿该方向单位长度的改变量。x方向的线应变εx等于微元体受力后x向变形沿x向求偏导数；y方向的线应变εy等于微元体受力后y向变形沿y向求偏导数。z方向的线应变εz等于微元体受力后z向变形沿z向求偏导数。γxy,γyz,γxz为绕三个坐标轴的切应变，即微元体两条相互垂直的棱边在变形后的直角改变量。绕z方向的切应变γxy,γyx等于微元体受力后y向变形v沿x向求偏导数与x向变形u沿y向求偏导数之和，即x轴与y轴所夹直角的角度改变量；绕x方向的切应变γyz,γzy等于微元体受力后z向变形w沿y向求偏导数与y向变形v沿z向求偏导数之和，即y轴与z轴所夹直角的角度改变量；绕y方向的切应变γzx,γxz等于微元体受力后z向变形w沿x向求偏导数与x向变形u沿z向求偏导数之和，即x轴与z轴所夹直角的角度改变量。

从式中不难看出，由变形得到应变是一个沿三个坐标轴求导的过程，由此特点可以联想到，如果将之前求得的引力各分量视为物体沿三个坐标方向的变形，求出在该变形下所产生的应变，即实现了对引力分量沿坐标轴求导。其中，线应变对应直线引力梯度，切应变对应交叉引力梯度。

下面以图1所示测试模型为例对本发明做进一步说明。选择存在明确理论解的矩形棱柱体作为测试模型。该模型尺寸：x向、y向-200m～200m，z向150m～470m，模型密度为1kg/m³。观测域为x向、y向-500m～500m，z向-50m～500m的空间范围，求在该观测域内测试模型的引力梯度。

首先设置恰当的单元类型并对测试模型进行网格划分。该模型相对规整，可划分为全六面体单元，对于形状复杂及含较多曲面的模型，也可以划分为四面体单元或两者的结合，以达到良好的模型拟合效果。该过程可以通过软件自动完成，相比以往手动分割模型，可以明显提高工作效率和模型拟合度，从根本上减小计算误差。

然后采用有限元仿真计算引力梯度：对模型施加边界条件及载荷，进行仿真计算。以带电绝缘体来模拟相同形状和体积的矩形棱柱形质量体，对观测域设置相对介电常数，并将质量体所在单元设置相应的体电荷密度。利用电场分析模块进行仿真计算，得到如图2所示的带电绝缘体周围的电场强度分布图。将上述结果乘以转换系数，即为矩形棱柱质量体周围的引力场分布。

为了求得引力梯度场，还需在直角坐标系中对各引力分量沿三个坐标轴进行求导。本发明借助静力分析中基于变形求应变的方法以实现该过程。将之前计算得到的各节点引力分量视为位移载荷边界条件施加于仿真模型。例如，当需要计算直线引力梯度uxx和交叉引力梯度uxy时，其计算公式如下：

其对应的几何方程如下

此时，只需将x向引力分量gx作为位移载荷施加于仿真模型，并将模型各节点的y向和z向位移设置为零，消除了式(8)中第二式等号右边第一项对结果的影响，基于此位移载荷仿真计算εx和γxy，即可求得uxx和uxy。利用类似方法，可以求得观测域内各阶梯度张量。选取观测域内z＝-2平面作为观测平面，该面上各阶梯度张量如图3a、3b、3c、3d、3e、3f所示。

本发明可采用如下方法评定其准确度。由于各阶梯度张量均基于相同的仿真模型及求解域且在一次计算中全部求出，故仅以uxy为考察对象进行误差分析，即可代表其余各阶梯度张量。将上述仿真得到的z＝-2平面上共计57601个节点的uxy仿真值逐一输出并与理论值进行比对，按公式(9)计算各点相对误差，其中uxy理论值可利用公式(4)求出。

在本发明通过绘图软件，将考察平面上各点的相对误差以云图方式显示，如图4所示。观察各阶梯度张量分量示意图不难发现，uxy在观测面上坐标靠近x＝0及y＝0的区域真实值趋近于零，无法用公式(9)计算其相对误差，故在图4中将该区域进行白化处理，这并不影响对仿真准确度的考量。从图中可以看出，uxy的最大相对误差为0.016，仿真误差小于2％。同样方法考察其余各阶梯度张量，可以得到相同的结论。该仿真精度可以通过设置单元节点通过观测面以减少插值计算误差或局部加密观测面附近网格等方式进一步提高。

需要强调的是，本发明所述的实施例是说明性的，而不是限定性的，因此本发明包括并不限于具体实施方式中所述的实施例，凡是由本领域技术人员根据本发明的技术方案得出的其他实施方式，同样属于本发明保护的范围。