有限圓柱體重力位三階梯度張量正演計(jì)算公式

冷 潔 邱 峰

1 海軍航空大學(xué)青島校區(qū)，青島市，2660412 大連測(cè)控技術(shù)研究所，大連市濱海街16號(hào)，116013

作為正反演計(jì)算中較為簡(jiǎn)單常用的模型構(gòu)建單元，有限圓柱體模型被廣泛應(yīng)用于地球物理學(xué)領(lǐng)域的正反演模擬計(jì)算、數(shù)據(jù)處理與解釋之中[1]。對(duì)于有限圓柱體模型，從重力位到重力二階梯度張量，其正演計(jì)算解析公式問(wèn)題已經(jīng)基本被學(xué)者解決[2-6]，而對(duì)于有限圓柱體模型的重力位三階梯度張量正演解析公式，還未見(jiàn)國(guó)內(nèi)外學(xué)者進(jìn)行報(bào)道或發(fā)表。鑒于重力曲率張量在未來(lái)多個(gè)領(lǐng)域的廣泛應(yīng)用前景，本文首先根據(jù)引力位計(jì)算公式推導(dǎo)得到有限圓柱體模型重力曲率張量的正演解析公式，而后通過(guò)理論分析與模型實(shí)驗(yàn)對(duì)比來(lái)分析驗(yàn)證所推公式的正確性。

1 公式推導(dǎo)

設(shè)有一長(zhǎng)為2L的剩余密度均勻的水平圓柱體，圓柱體的中心坐標(biāo)為(x0,y0,z0)，其截面積為S，剩余密度為ρ，坐標(biāo)系與水平圓柱體的關(guān)系如圖1所示，Y軸平行于水平圓柱體的走向，X軸垂直于柱軸。將其看作質(zhì)量集中于軸線Q(ξ,η,ζ)上的有限長(zhǎng)水平物質(zhì)線，故水平圓柱體任一體積元dv=dξdηdζ=S·dη。設(shè)空間中任一觀測(cè)點(diǎn)坐標(biāo)為P(x,y,z)，則根據(jù)牛頓萬(wàn)有引力定律可得有限長(zhǎng)水平圓柱體在P點(diǎn)產(chǎn)生的重力位為[7]：

圖1 圓柱體幾何示意Fig.1 Cylinder geometry schematic

(1)

式中，ρ為密度，G=6.672×10-11m3/kg·s2為萬(wàn)有引力常數(shù)，r=[(ξ-x)2+(η-y)2+(ζ-z)2]1/2為觀測(cè)點(diǎn)與體積元dξdηdζ之間的距離。

對(duì)式(1)積分要用到的積分公式為：

(2)

令ξ-x=a，η-y=b，ζ-z=c。

重力二階梯度是重力位各方向上的二階導(dǎo)數(shù)，其張量形式為：

(3)

根據(jù)眾多學(xué)者對(duì)重力梯度張量正演模擬和反演解釋的研究[2-6]可得，有限圓柱體的重力位梯度張量解析公式為：

Uxx=GρS·

(4a)

(4b)

(4c)

(4d)

(4e)

Uzz=GρS·

(4f)

重力曲率張量是重力位在各個(gè)方向上的三階導(dǎo)數(shù)，根據(jù)張量的對(duì)稱性可知，Uijk=Ujki=Ukij,i、j、k∈{x,y,z}。重力曲率張量的正演解析公式可由重力位梯度張量求導(dǎo)得到，也可通過(guò)對(duì)重力位進(jìn)行3次求偏導(dǎo)得到。

由式(2)可得：

(5)

(6a)

(6b)

(6c)

(6d)

(6e)

(6f)

(6g)

(6h)

(6i)

(6j)

2 正演計(jì)算公式的正確性檢驗(yàn)

為了檢驗(yàn)本文所推導(dǎo)公式的正確性，通過(guò)對(duì)理論模型重力位二階梯度張量異常結(jié)果進(jìn)行中心差分計(jì)算，得到所用模型的重力位三階梯度張量結(jié)果，而后與本文解析公式計(jì)算結(jié)果進(jìn)行對(duì)比。

2.1 理論模型的重力三階梯度張量分布

首先給出模型參數(shù)：圓柱體的橫截面半徑R=4 m，圓柱體長(zhǎng)2L=40 m，剩余密度為2.67 g/cm3，其中心坐標(biāo)為(0, 0, 40 m)，X、Y坐標(biāo)范圍均為-100～100 m。地面測(cè)網(wǎng)大小為200 m×200 m，網(wǎng)格間距為1 m×1 m。利用式(6)計(jì)算得到該有限圓柱體模型的重力三階梯度張量分布，如圖2所示。

黑色矩形為有限圓柱體模型的水平位置；白色實(shí)線為圖3剖面的水平位置；1 pMKS=10-12/ms2圖2 有限圓柱體的重力三階梯度張量分布Fig.2 Distribution of third-order gradient tensor of the gravitational potential caused by a cylinde

2.2 解析解的正確性檢驗(yàn)

由導(dǎo)數(shù)的定義可知，重力三階梯度張量可由其二階梯度張量通過(guò)中心差分得到，故可通過(guò)與二階梯度張量差分計(jì)算結(jié)果進(jìn)行對(duì)比，以檢驗(yàn)本文三階張量正演計(jì)算公式的正確性。圖3為圖2所示剖面的重力三階梯度張量曲線的差分結(jié)果，差分點(diǎn)距設(shè)置為1 m。由圖可知，本文所計(jì)算的差分結(jié)果都逼近于重力三階梯度張量計(jì)算結(jié)果。圖4為所計(jì)算的差分結(jié)果與理論結(jié)果之間的誤差曲線。由圖可知，整體誤差都較小，并且誤差在峰值處達(dá)到最大，證明了本文推導(dǎo)公式的正確性。

圖3 解析表達(dá)式與中心差分的計(jì)算結(jié)果對(duì)比Fig.3 Comparison between the calculating results by analytic expression and central difference

圖4 解析表達(dá)式與中心差分的計(jì)算結(jié)果誤差Fig.4 The calculating error results by analytic expression and central difference

根據(jù)擾動(dòng)重力位在場(chǎng)源外部空間應(yīng)當(dāng)滿足拉普拉斯方程，重力位三階梯度張量的部分分量應(yīng)當(dāng)滿足Uxxx+Uxyy+Uxzz=0、Uxxy+Uyyy+Uyzz=0、Uxxz+Uyyz+Uzzz=0[8]。基于該原理檢驗(yàn)所推公式的正確性，計(jì)算結(jié)果如圖5所示?？梢钥闯觯?jì)算結(jié)果滿足上述3個(gè)方程，誤差均為0，這從另外一個(gè)方面證明了本文推導(dǎo)公式的正確性。

圖5 重力位三階梯度張量部分元素相加結(jié)果Fig.5 The partial element addition results of third-order gradient tensor

3 結(jié) 語(yǔ)

鑒于有限圓柱體模型在地球物理學(xué)領(lǐng)域的正反演模擬與計(jì)算、數(shù)據(jù)處理與解釋中的重要性，本文根據(jù)引力位解析公式，經(jīng)過(guò)推導(dǎo)得到有限圓柱體重力曲率張量的正演解析公式，并基于理論模型實(shí)驗(yàn)，從對(duì)所推解析公式的正演結(jié)果與重力位梯度張量中心差分的計(jì)算結(jié)果進(jìn)行對(duì)比以及檢驗(yàn)正演計(jì)算結(jié)果是否滿足拉普拉斯方程2個(gè)方面驗(yàn)證了所推導(dǎo)公式的正確性。本文研究成果可為今后重力位三階梯度張量的仿真模擬、實(shí)測(cè)數(shù)據(jù)的處理與轉(zhuǎn)換計(jì)算以及反演解釋提供理論基礎(chǔ)。