邱 峰

(江西省水利科學(xué)院 建材與巖土研究所,江西 南昌 330029)

0 引言

重力位三階梯度張量是重力位在坐標(biāo)軸3個(gè)方向的三階導(dǎo)數(shù),表示重力梯度張量各個(gè)元素在3個(gè)方向上的空間變化率,相比傳統(tǒng)的重力異常數(shù)據(jù)以及重力梯度異常數(shù)據(jù),其對(duì)于微弱的密度差異引起的變化更為敏感。在重力勘探及大地測(cè)量領(lǐng)域,測(cè)量地球的重力位場(chǎng)的三階導(dǎo)數(shù)及其時(shí)變場(chǎng)最早是由Balakin等[1]于1997年在Dulkyn項(xiàng)目中提出。之后基于此思想與概念,DiFrancesco等[2]在2009年提出重力三階梯度張量的航空測(cè)量概念,Brieden等[3]又在2010年提出利用衛(wèi)星測(cè)量地球時(shí)變重力場(chǎng)的計(jì)劃。同時(shí),Fantino等[4]和Casotto等[5]以及Fukushima[6]從2009～2013年相繼發(fā)表了闡述重力位一階至三階梯度的球諧表達(dá)式及其計(jì)算問(wèn)題的相關(guān)系列論文。而Rosi等[7]于2015年在實(shí)驗(yàn)室內(nèi)利用冷原子干涉技術(shù)測(cè)量重力位三階垂直梯度的成功實(shí)現(xiàn),使得重力位三階梯度張量研究取得突破性的進(jìn)展。此外,由于重力位三階垂向梯度對(duì)局部質(zhì)量源具有高度敏感性,也被用于實(shí)驗(yàn)室測(cè)量萬(wàn)有引力常數(shù)[8];并且,由于地球重力位三階垂直梯度具有不為零的特點(diǎn),使得重力垂直梯度隨高度變化存在非線性特征[9]。目前,重力位三階梯度張量測(cè)量技術(shù)已經(jīng)初見(jiàn)雛形,隨著其將來(lái)的進(jìn)一步發(fā)展,勢(shì)必在重力勘探以及大地測(cè)量等領(lǐng)域引起新一輪的技術(shù)革命。

DEXP(depth from extreme point)變換成像法是近幾年剛被提出并且發(fā)展較快的一種方法,由Fedi[10]提出,由于在計(jì)算過(guò)程中深度加權(quán)函數(shù)加入了異常體對(duì)應(yīng)的構(gòu)造指數(shù),使得相比于其他成像結(jié)果更為準(zhǔn)確。該方法的適用性比較廣,不僅可以用于重力異常不同階次導(dǎo)數(shù)來(lái)估計(jì)異常體的物性參數(shù)(如場(chǎng)源深度、剩余密度及構(gòu)造指數(shù)等),也可以應(yīng)用于自然電位數(shù)據(jù)處理,并有較好的結(jié)果[11]。由于重力位三階梯度張量具有7個(gè)獨(dú)立元素,因此其空間分布形態(tài)豐富多樣,能夠獲取比傳統(tǒng)重力異常以及重力梯度張量更多的場(chǎng)源信息[12],并且采用重力位三階梯度張量異常進(jìn)行DEXP定量解釋具有更強(qiáng)的抗噪能力、受背景場(chǎng)影響也更小[13]。

由于重力位三階梯度張量勘探方法還沒(méi)有提出完善的測(cè)量體系,尚處于概念提出與理論基礎(chǔ)研究階段,因此在實(shí)際數(shù)據(jù)資料處理中,只能通過(guò)基于重力異?；蛑亓μ荻葟埩繑?shù)據(jù)進(jìn)行數(shù)據(jù)轉(zhuǎn)換得到所需要的重力位三階梯度張量數(shù)據(jù)。波數(shù)域數(shù)據(jù)轉(zhuǎn)換計(jì)算方法是其中較為常用且有效的方法之一[14-15]。重力位三階梯度張量波數(shù)域轉(zhuǎn)換計(jì)算方法還未見(jiàn)國(guó)內(nèi)外學(xué)者進(jìn)行報(bào)道或發(fā)表,因此本文首先推導(dǎo)重力位三階梯度張量的波數(shù)域轉(zhuǎn)換計(jì)算公式,然后利用理論模型進(jìn)行試驗(yàn)對(duì)比驗(yàn)證,最后將該方法應(yīng)用到Vinton鹽丘區(qū)DEXP數(shù)據(jù)處理中,取得較好的結(jié)果。

1 方法原理

1.1 重力位三階梯度張量異常的波數(shù)域轉(zhuǎn)換計(jì)算方法

1.1.1 基于重力異常的波數(shù)域轉(zhuǎn)換計(jì)算方法

對(duì)于自由空間的重力矢量g(x,y,z) = (gx,gy,gz),其3個(gè)分量是連續(xù)可微的。Mickus等[14]根據(jù)拉普拉斯方程得到了重力矢量各分量的傅里葉變換式:

(1)

重力位三階梯度張量是重力三分量在x、y以及z方向上的二階導(dǎo)數(shù),即Wαβγ=?2gα/(?β?γ),?{α,β,γ}∈{x,y,z},其僅有7個(gè)獨(dú)立元素。因此根據(jù)式(1),由重力異常數(shù)據(jù)得到重力位三階梯度張量數(shù)據(jù)的表達(dá)式為:

(2)

1.1.2 基于重力梯度張量的波數(shù)域轉(zhuǎn)換計(jì)算方法

根據(jù)式(1)同樣可得重力梯度張量的傅里葉轉(zhuǎn)換表達(dá)式,Mickus等[14]已經(jīng)得到基于重力異常數(shù)據(jù)波數(shù)域轉(zhuǎn)換重力梯度張量的表達(dá)式:

Tij= FFT-1{[K(k)]Gz(k)},

(3)

其中:

。

(4)

由于重力三階梯度張量是重力梯度張量在坐標(biāo)軸3個(gè)方向上的導(dǎo)數(shù),所以結(jié)合式(4)即可得基于重力梯度張量數(shù)據(jù)波數(shù)域轉(zhuǎn)換重力三階梯度張量的表達(dá)式:

(5)

其中:Гij為重力梯度張量各分量的二維離散快速傅里葉變換。

1.2 重力位三階梯度張量異常的比值DEXP定量解釋方法

根據(jù)Fedi[10]提出的DEXP變換方法理論可知,重力位異常在坐標(biāo)軸3個(gè)方向上的p階導(dǎo)數(shù)的DEXP變換函數(shù)可表示為:

Ωp=zNp/2fp,

(6)

其中:zNp/2為p階DEXP變換的深度加權(quán)函數(shù);fp為重力位異常p階導(dǎo)數(shù)。

由于重力位三階梯度張量異常為重力位在坐標(biāo)軸x、y、z方向得到三階導(dǎo)數(shù),所以式(6)中Np可以寫(xiě)成N3=N+3,因此重力位三階梯度張量各分量的DEXP變換函數(shù)為:

(7)

2 模型試驗(yàn)

首先令均勻直立棱柱體模型(如圖1所示)的長(zhǎng)(a)、寬(b)及高(c)分別為150 m、100 m以及60 m,剩余密度為0.5 g/cm3,中心埋深為100 m,計(jì)算點(diǎn)距為2 m,觀測(cè)面為地面即z=0 m。直立棱柱體重力位一階至三階重力梯度正演公式眾多文獻(xiàn)已給出[16],因此可得模型重力位一階至三階重力梯度張量異常,如圖2所示。

圖1 直立長(zhǎng)方體模型幾何示意Fig.1 Geometric schematic diagram of upright cuboid model

圖2 模型重力位一階至三階張量異常空間分布Fig.2 The gravity tensor anomalies from first-order to third-order of upright cuboid model

2.1 基于重力異常的重力位三階梯度張量異常轉(zhuǎn)換計(jì)算

根據(jù)式(2)可知,為得到重力位三階梯度張量各分量結(jié)果,首先須將重力異常模擬數(shù)據(jù)進(jìn)行二維離散快速傅里葉變換(FFT),再乘以各分量對(duì)應(yīng)的轉(zhuǎn)化系數(shù)計(jì)算各自的頻譜,最后再對(duì)所得頻譜分別進(jìn)行二維離散快速傅里葉反變換(IFFT),即可得到模擬數(shù)據(jù)的轉(zhuǎn)換結(jié)果。本文提出的重力位三階梯度張量計(jì)算算法使用到了二維離散快速傅里葉變換方法,因此FFT轉(zhuǎn)換計(jì)算結(jié)果會(huì)存在一定的誤差。為了直觀地說(shuō)明重力位三階梯度張量FFT轉(zhuǎn)換計(jì)算結(jié)果與模型正演計(jì)算結(jié)果的差異,本文計(jì)算并可視化了模型對(duì)應(yīng)的每個(gè)分量之間的誤差,如圖3所示。由圖可知,Wxxy、Wxyy、Wxyz分量的誤差較小,誤差絕對(duì)值最大不超過(guò)1.6 pMKS;但剩余的7個(gè)分量(Wxxx、Wxxz、Wxzz、Wxyy、Wyyz、Wyzz、Wzzz),由于FFT計(jì)算過(guò)程中受邊界效應(yīng)影響較大,使得分布在計(jì)算邊界的數(shù)據(jù)結(jié)果誤差也比較大,其誤差較大的數(shù)據(jù)值分布在x、y剖面上-400～ -380 m和380～ 400 m范圍之間,如果去除這些受邊界效應(yīng)影響較大的數(shù)據(jù),則重力位三階梯度張量FFT轉(zhuǎn)換與正演計(jì)算之間的差異較小。

圖3 基于重力異常的重力位三階梯度張量FFT轉(zhuǎn)換與正演計(jì)算各分量之間的差異Fig.3 The deviations between the FFT-derived and model-calculated third-order gradient tensor of gravitational potential from gravity value

此外,通過(guò)均方根誤差(RMS)來(lái)量化本文所使用方法存在的誤差,其計(jì)算公式如下:

?{α,β,γ}∈{x,y,z} ,

(8)

表1為正演模型重力位三階梯度張量的均方根誤差,其中去掉部分邊界畸變值計(jì)算的是-390～390 m范圍之間7個(gè)分量(Wxxx、Wxxz、Wxzz、Wxyy、Wyyz、Wyzz、Wzzz)的均方根誤差。由表1可知,均方根誤差比較大的是受邊界影響較大的7個(gè)分量(Wxxx、Wxxz、Wxzz、Wxyy、Wyyz、Wyzz、Wzzz),其均方根誤差均超過(guò)7 pMKS,RMS最大的為Wzzz分量;剩余的3個(gè)分量計(jì)算的均方根誤差較小,均不超過(guò)0.2 pMKS,其中RMS最小的為Wxyz分量。對(duì)于去掉部分邊界畸變值的計(jì)算結(jié)果,每個(gè)分量的均方根誤差不超過(guò)2.55 pMKS,其中RMS最大的為Wxxx分量,RMS最小的為Wxyz分量。

表1 基于重力異常的重力三階梯度張量異常轉(zhuǎn)換計(jì)算結(jié)果與理論值間的均方根誤差Table 1 RMS error between model and FFT third-order gradient tensor of gravitational potential

2.2 基于重力梯度張量異常的重力位三階梯度張量異常轉(zhuǎn)換計(jì)算

同理,根據(jù)式(5)可得出基于重力梯度張量波數(shù)域轉(zhuǎn)換計(jì)算的重力位三階梯度張量異常。為直觀說(shuō)明重力位三階梯度張量FFT轉(zhuǎn)換計(jì)算結(jié)果與模型正演計(jì)算結(jié)果的差異,同樣計(jì)算并可視化了模型對(duì)應(yīng)的每個(gè)分量之間的誤差,如圖4所示。由圖可知,每個(gè)分量轉(zhuǎn)換計(jì)算結(jié)果與模型正演計(jì)算結(jié)果的誤差非常小,誤差絕對(duì)值最大不超過(guò)2.2 pMKS,其中FFT轉(zhuǎn)換計(jì)算誤差最大的為Wzzz分量,誤差最小的為Wxyz分量。

圖4 基于重力梯度張量數(shù)據(jù)的重力位三階梯度張量FFT轉(zhuǎn)換與正演計(jì)算各分量之間的差異Fig.4 The deviations between the FFT-derived and model-calculated third-order gradient tensor of gravitational potential from gravity gradient tensor

通過(guò)式(8),可以計(jì)算出基于重力梯度張量數(shù)據(jù)的重力位三階梯度張量FFT轉(zhuǎn)換計(jì)算結(jié)果與理論值之間的均方根誤差,如表2所示。由表可知,每個(gè)分量所計(jì)算的均方根誤差都比較小,RMS最大不超過(guò)1 pMKS,對(duì)應(yīng)為Wzzz分量,RMS最小為0.04 pMKS,對(duì)應(yīng)為Wxyz分量。并且與重力異常數(shù)據(jù)的FFT轉(zhuǎn)換計(jì)算結(jié)果相比,基于重力梯度張量的換算結(jié)果受邊界效應(yīng)的影響非常小,使得每個(gè)分量與理論值之間的誤差也非常小,對(duì)應(yīng)的均方根誤差也很小,更接近理論值。

表2 基于重力梯度張量的重力三階梯度張量異常轉(zhuǎn)換計(jì)算結(jié)果與理論值間的均方根誤差Table 2 RMS error between model and FFT third-order gradient tensor of gravitational potential

3 實(shí)際應(yīng)用

3.1 數(shù)據(jù)及其預(yù)處理

本文所采用的文頓鹽丘重力梯度張量(FTG)實(shí)測(cè)數(shù)據(jù)由Bell Geospace公司在2008年利用Air-FTG全張量重力梯度儀測(cè)量得到。其中,該測(cè)區(qū)范圍為北緯30.07°～30.23°、西經(jīng)93.53°～93.66°,總測(cè)線長(zhǎng)度為1 087.5 km,測(cè)區(qū)面積為192.6 km2。在實(shí)際航空測(cè)量中,其飛行測(cè)線方向?yàn)镾N向,其中,中間兩側(cè)53條測(cè)線間距為250 m,中間部分測(cè)線間距為125 m,航測(cè)飛行高度設(shè)定為80 m。由于數(shù)據(jù)量較大,為方便數(shù)據(jù)資料處理分析,本文只選取了面積較小的文頓鹽丘測(cè)區(qū)。在實(shí)際異常數(shù)據(jù)資料的處理中,噪聲以及區(qū)域背景場(chǎng)的存在往往會(huì)對(duì)數(shù)據(jù)處理結(jié)果產(chǎn)生較大的誤差,尤其是背景場(chǎng)。因此,為提高后續(xù)數(shù)據(jù)處理結(jié)果的精度,先對(duì)實(shí)測(cè)重力異常以及重力梯度張量數(shù)據(jù)進(jìn)行去噪處理,然后對(duì)去噪后的數(shù)據(jù)進(jìn)行濾波,移除背景場(chǎng)。圖5a～f為處理后的重力梯度張量異常數(shù)據(jù)。

a—Txx異常;b—Txy異常;c—Txz異常;d—Tyy異常;e—Tyz異常;f—Tzz異常

3.2 重力位三階梯度張量異常的轉(zhuǎn)換計(jì)算結(jié)果

根據(jù)基于模型重力異常與重力梯度張量數(shù)據(jù)轉(zhuǎn)換的重力位三階梯度張量對(duì)比結(jié)果可知,由重力梯度張量轉(zhuǎn)換的重力三階梯度張量受邊界效應(yīng)的影響較低,計(jì)算結(jié)果精度更高。因此在本章實(shí)際數(shù)據(jù)資料處理中,結(jié)合前文數(shù)據(jù)波數(shù)域轉(zhuǎn)換理論,采用如圖5所示的經(jīng)過(guò)處理后的Vinton鹽丘區(qū)重力梯度張量數(shù)據(jù)進(jìn)行波數(shù)域轉(zhuǎn)換,即可得到該研究區(qū)較高精度的重力位三階梯度張量異常數(shù)據(jù),結(jié)果如圖6所示。圖中黑色實(shí)線為后續(xù)數(shù)據(jù)DEXP變換處理所選取的剖面位置。

a—Wxxx;b—Wxxy;c—Wxxz;d—Wxyy;e—Wxyz;f—Wxzz;g—Wyyy;h—Wyyz;i—Wyzz;j—Wzzz

3.3 重力位三階梯度張量異常的比值DEXP定量解釋結(jié)果

根據(jù)DEXP理論,在對(duì)Vinton鹽丘區(qū)實(shí)測(cè)數(shù)據(jù)進(jìn)行DEXP變換處理時(shí),需要預(yù)先知道該地區(qū)地下異常體的構(gòu)造指數(shù)。根據(jù)Li[17]的研究可知,該鹽丘體重力異常的構(gòu)造指數(shù)為1.01,從而可以確定該鹽丘體的重力位三階梯度張量數(shù)據(jù)的構(gòu)造指數(shù)為3.01。結(jié)合式(7),將波數(shù)域轉(zhuǎn)換所得的重力位三階梯度張量異常數(shù)據(jù)進(jìn)行DEXP變換,即可得各分量所選剖面上的深度成像結(jié)果,如圖7所示。

a—Ωxxx;b—Ωxxy;c—Ωxxz;d—Ωxyy;e—Ωxzz;f—Ωyyy;g—Ωyyz;h—Ωyzz;i—Ωzzz

由圖7可知,能反映鹽丘區(qū)深度信息的有Ωxxz結(jié)果(圖7c)、Ωyyz結(jié)果(圖7g)及Ωzzz結(jié)果(圖7i),其各自的深度成像結(jié)果分別為580 m、660 m及600 m,減去航測(cè)飛行高度80 m,其反映的鹽丘區(qū)平均埋藏深度約為533 m;Ωxxx結(jié)果(圖7a)反映的是鹽丘體EW向長(zhǎng)度,其結(jié)果為1 350 m,Ωyyy結(jié)果(圖7f)可以反映鹽丘體SN向長(zhǎng)度,其對(duì)應(yīng)結(jié)果為1 125 m。Oliveria等[18]給出的鹽丘區(qū)底面深度結(jié)果是460 m,Geng等[19]給出的底面深度結(jié)果是650 m,Ennen[20]給出Vinton鹽丘區(qū)EW向長(zhǎng)度為1 520 m,SN向長(zhǎng)度為1 280 m。結(jié)合已知資料,認(rèn)為本文所得結(jié)果是合理的。

4 結(jié)論

本文經(jīng)過(guò)推導(dǎo)得到基于重力異常及重力梯度張量模擬數(shù)據(jù)的重力位三階梯度張量波數(shù)域轉(zhuǎn)換公式,并基于直立棱柱體模型進(jìn)行試驗(yàn)驗(yàn)證,結(jié)果表明所采用的波數(shù)域轉(zhuǎn)換方法是可行的,并且相對(duì)于重力異常數(shù)據(jù),基于重力梯度張量數(shù)據(jù)的波數(shù)域轉(zhuǎn)換結(jié)果精度更高。同時(shí),將Vinton鹽丘區(qū)實(shí)測(cè)重力梯度張量數(shù)據(jù)用于波數(shù)域轉(zhuǎn)換方法中,得到該地區(qū)重力位三階梯度張量數(shù)據(jù),并將所得結(jié)果用于DEXP方法的數(shù)據(jù)解釋中,得到該鹽丘區(qū)的底面深度為533 m,EW向長(zhǎng)度為1 350 m,SN向長(zhǎng)度為1 125 m,所得結(jié)果與前人的研究結(jié)果相符合。

致謝：感謝美國(guó)Bell Geospace公司提供全張量重力梯度數(shù)據(jù)。