冷 潔 邱 峰

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

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

1 公式推導(dǎo)

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

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

(1)

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

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

(2)

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

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

(3)

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

Uxx=GρS·

(4a)

(4b)

(4c)

(4d)

(4e)

Uzz=GρS·

(4f)

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

由式(2)可得：

(5)

(6a)

(6b)

(6c)

(6d)

(6e)

(6f)

(6g)

(6h)

(6i)

(6j)

2 正演計算公式的正確性檢驗

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

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

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

黑色矩形為有限圓柱體模型的水平位置；白色實線為圖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 解析解的正確性檢驗

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

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

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

根據(jù)擾動重力位在場源外部空間應(yīng)當滿足拉普拉斯方程，重力位三階梯度張量的部分分量應(yīng)當滿足Uxxx+Uxyy+Uxzz=0、Uxxy+Uyyy+Uyzz=0、Uxxz+Uyyz+Uzzz=0[8]?；谠撛頇z驗所推公式的正確性，計算結(jié)果如圖5所示。可以看出，計算結(jié)果滿足上述3個方程，誤差均為0，這從另外一個方面證明了本文推導(dǎo)公式的正確性。

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

3 結(jié) 語

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