石澤玉，張志厚,2，劉鵬飛，范祥泰

(1.西南交通大學(xué) 地球科學(xué)與環(huán)境工程學(xué)院，四川 成都 611756; 2.西南交通大學(xué) 高速鐵路線路工程教育部重點(diǎn)實(shí)驗(yàn)室，四川 成都 610031)

0 引言

重力及其梯度異常是由于地球局部質(zhì)量分布不均勻而產(chǎn)生的[1]。重力及其梯度異常正演是根據(jù)已知地質(zhì)體的形狀、產(chǎn)狀和剩余密度等來(lái)計(jì)算異常的分布規(guī)律，是重力勘探定量解釋的基礎(chǔ)，其對(duì)重力反演[2-4]及輔助導(dǎo)航[5]都具有非常重要的意義。大規(guī)模重力及其梯度異常的正演速度決定了反演的可行性[2]，也成為海量重力數(shù)據(jù)處理迫切需要解決的現(xiàn)實(shí)問(wèn)題[6]。

為了提高重力及其梯度正演的計(jì)算效率，近年來(lái)眾多學(xué)者提出了多種重力或其梯度的正演方法與技術(shù)，如Li等[7]應(yīng)用小波變換和閾值小波系數(shù)組合壓縮重力異常正演靈敏度矩陣來(lái)減少計(jì)算量，但仍需要對(duì)靈敏度矩陣進(jìn)行計(jì)算；秦朋波等[8]發(fā)現(xiàn)了規(guī)則網(wǎng)絡(luò)情況下核函數(shù)的對(duì)稱(chēng)性，并提出了一種快速計(jì)算靈敏度矩陣的計(jì)算方法；熊光楚[9]推導(dǎo)了長(zhǎng)方體單元重力異常的傅里葉變換表達(dá)式；Shin等[10]提出了一種基于快速傅里葉變化的頻率域方法，該方法通過(guò)改進(jìn)已有算法加強(qiáng)了運(yùn)算中的周期性，減少了運(yùn)算中的數(shù)量級(jí)，在大型數(shù)據(jù)集進(jìn)行計(jì)算時(shí)可以有效提升計(jì)算效率；Wu等[11]受偏移抽樣技術(shù)的啟發(fā)，提出了重力異常正演的高斯快速傅里葉變換(Gauss-FFT)計(jì)算方法，相比標(biāo)準(zhǔn)的快速傅里葉(FFT)正演方法，Gauss-FFT計(jì)算精度更高、速度更快；Ren等[12]引入了自適應(yīng)快速多極子方法可對(duì)任意起伏地形的重力異常進(jìn)行快速正演。在空間域，姚長(zhǎng)利等[2]提出幾何格架函數(shù)的方法，首先將場(chǎng)源劃分成若干規(guī)則長(zhǎng)方體單元，然后計(jì)算部分測(cè)點(diǎn)的格架函數(shù)并存儲(chǔ)，其余測(cè)點(diǎn)格架函數(shù)可利用平移等效性和互換對(duì)稱(chēng)性直接調(diào)用，從而大大減少了計(jì)算量和存儲(chǔ)量；陳召曦等[3]在幾何格架函數(shù)方法的基礎(chǔ)上提出了多核CPU加速并行計(jì)算的方法，以此提升正反演速度。張志厚等[4]提出了網(wǎng)格點(diǎn)幾何格架函數(shù)的概念，該方法實(shí)質(zhì)上是對(duì)幾何格架函數(shù)方法進(jìn)行加速。兩者不同點(diǎn)在于幾何格架函數(shù)的概念是相對(duì)于長(zhǎng)方體單元，網(wǎng)格點(diǎn)幾何格架函數(shù)的概念是相對(duì)于長(zhǎng)方體的角點(diǎn)，其改進(jìn)在于避免了網(wǎng)格點(diǎn)幾何格架函數(shù)的多次重復(fù)計(jì)算。因此，當(dāng)剖分單元足夠大時(shí)，理論上網(wǎng)格點(diǎn)幾何格架函數(shù)的計(jì)算效率能夠提高近8倍。

但隨著航空地球物理的發(fā)展，大面積海量重力數(shù)據(jù)面臨著高精度快速處理的挑戰(zhàn)。而以上方法只能對(duì)較少的數(shù)據(jù)量進(jìn)行處理，如采用網(wǎng)格點(diǎn)幾何格架函數(shù)的策略[4]對(duì)256×256×15的單元體進(jìn)行正演至少需要3個(gè)多小時(shí)，當(dāng)剖分單元擴(kuò)大4倍到512×512×15時(shí)，其正演時(shí)間呈指數(shù)上漲，即在普通計(jì)算機(jī)上很難實(shí)現(xiàn)一次正演，難以滿(mǎn)足實(shí)際生產(chǎn)的需求。

受航空電磁Moving-footprint大尺度模型分解正演技術(shù)的啟發(fā)，本文提出了一種重力及其梯度異常正演的Moving-footprint大尺度模型分解計(jì)算方法。

文中將地下半空間規(guī)則剖分成若干長(zhǎng)方體單元，某觀測(cè)點(diǎn)重力及其梯度異常主要為其正下方一定范圍內(nèi)(子空間)的物性單元體產(chǎn)生，當(dāng)觀測(cè)點(diǎn)移動(dòng)時(shí)，子空間跟隨移動(dòng)，即為“Moving-footprint”。文中劃分了不同尺度的子空間進(jìn)行正演計(jì)算，并將計(jì)算結(jié)果與理論結(jié)果進(jìn)行對(duì)比，以此來(lái)驗(yàn)證方法的適用性。

1 重力及其梯度異常正演

重力及其梯度正演是將地下半空間剖分成若干個(gè)長(zhǎng)方體單元(如圖1所示)，然后計(jì)算每一個(gè)長(zhǎng)方體單元在觀測(cè)點(diǎn)觀測(cè)到的異常值，再將每一個(gè)長(zhǎng)方體的異常值求和，得到的結(jié)果即為地下半空間內(nèi)的異常[13]。單個(gè)長(zhǎng)方體單元在觀測(cè)點(diǎn)產(chǎn)生的重力異常及重力梯度異常理論計(jì)算表達(dá)為[14]式(1)～(7)。

圖1 地下半空間單元?jiǎng)澐諪ig.1 Underground half-space unit division

(1)

(2)

(3)

(4)

(5)

(6)

(7)

2 基于Moving-footprint技術(shù)的重力及其梯度異常正演方法

Moving-footprint即移動(dòng)腳印技術(shù)，該技術(shù)廣泛應(yīng)用于航空電磁探測(cè)領(lǐng)域，實(shí)現(xiàn)了大規(guī)模航空電磁數(shù)據(jù)的快速正反演[15-16]。Yin等[16]將航空電磁系統(tǒng)的Moving-footprint定義為地下有限導(dǎo)電半空間中的某一子空間，認(rèn)為該子空間整體電磁響應(yīng)約為地下半空間總電磁響應(yīng)的90%，當(dāng)發(fā)射接收系統(tǒng)移動(dòng)，該子空間隨之移動(dòng)，即為航空電磁系統(tǒng)的Moving-footprint大尺度模型分解技術(shù)。在重力異常正演計(jì)算中應(yīng)用Moving-footprint技術(shù)，即只考慮在距離觀測(cè)點(diǎn)一定的范圍內(nèi)的網(wǎng)格點(diǎn)(子空間內(nèi))在觀測(cè)點(diǎn)產(chǎn)生的異常影響(如圖2所示)，而超出選定的子空間范圍的網(wǎng)格點(diǎn)由于距離太遠(yuǎn)，在觀測(cè)點(diǎn)產(chǎn)生的異常影響十分微弱，在計(jì)算中可以忽略不計(jì)。文中方法與航空電磁領(lǐng)域的Moving-footprint不同之處是選取的子空間在深度上與地下半空間整體深度一致，而不是航空電磁的淺層部分深度。

注：紅線包圍區(qū)域?yàn)樽涌臻gnote：the area enclosed by the red line is the subspace圖2 地下網(wǎng)格體子空間劃分示意Fig.2 Underground grid subspace division schematic diagram

本文所提Moving-footprint的重力及其梯度異常正演方法是在網(wǎng)格點(diǎn)幾何格架函數(shù)技術(shù)的基礎(chǔ)上進(jìn)行了改進(jìn)。因此，首先在全空間內(nèi)選取子空間的大?。蝗缓笥?jì)算子空間內(nèi)部分測(cè)點(diǎn)的網(wǎng)格點(diǎn)的幾何格架函數(shù)，并存儲(chǔ)以備調(diào)用；隨后計(jì)算某觀測(cè)點(diǎn)異常時(shí)，先判斷該觀測(cè)點(diǎn)所在的子空間，再調(diào)用網(wǎng)格點(diǎn)格架函數(shù)和該子空間單元體的剩余密度，求和后獲得該觀測(cè)點(diǎn)異常；最后，利用Moving-footprint完成觀測(cè)區(qū)域所有點(diǎn)的重力及其梯度異常正演，主要計(jì)算步驟：

1)將地下半空間剖分為規(guī)則的長(zhǎng)方體單元，并對(duì)單元體的剩余密度進(jìn)行賦值；

2)確定子空間大小，計(jì)算子空間的網(wǎng)格點(diǎn)格架函數(shù)并進(jìn)行存儲(chǔ)以備調(diào)用；

3)根據(jù)觀測(cè)點(diǎn)確定子空間的相對(duì)位置，利用平移等效性和對(duì)稱(chēng)互換性調(diào)用網(wǎng)格點(diǎn)幾何格架函數(shù)，并利用式(1)～(7)完成該觀測(cè)點(diǎn)重力及其梯度異常的正演計(jì)算；

4)觀測(cè)點(diǎn)移動(dòng)，子空間隨之移動(dòng)，利用步驟3完成移動(dòng)觀測(cè)點(diǎn)重力及其梯度異常的正演計(jì)算；最終完成所有觀測(cè)點(diǎn)的正演計(jì)算。

基于Moving-footprint技術(shù)，只計(jì)算對(duì)觀測(cè)點(diǎn)起主要貢獻(xiàn)的長(zhǎng)方體單元產(chǎn)生的異常，從而大大減少了總運(yùn)算量，提高了計(jì)算效率。

3 模型實(shí)驗(yàn)

為了檢驗(yàn)本文所提方法的計(jì)算效果，將地下半空間剖分為256×256×15個(gè)長(zhǎng)方體單元，長(zhǎng)方體單元大小為0.1 km×0.1 km×0.1 km。采用4個(gè)長(zhǎng)方體組合模型進(jìn)行檢驗(yàn)，長(zhǎng)方體組合模型的剩余密度都為0.5 g/cm3，長(zhǎng)方體組合模型大小為4.0 km×4.0 km×0.5 km，其中心點(diǎn)坐標(biāo)分別為(13.0 km, 13.0 km, 0.75 km)、(20.0 km, 6.0 km, 0.75 km)、(6.0 km, 6.0 km, 0.75 km)及(20.0 km, 20.0 km, 0.75 km)，模型示意如圖3所示，其中，正演計(jì)算點(diǎn)距為0.1 km×0.1 km。

注：紅線包圍區(qū)域?yàn)橛?jì)算子空間，藍(lán)色立方體區(qū)域?yàn)楫惓ｓw區(qū)域note：the area surrounded by the red line is the calculation subspace, and the blue cube is the gravity anomaly area圖3 模型示意Fig.3 Model diagram

通過(guò)式(1)～(7)利用網(wǎng)格點(diǎn)幾何格架函數(shù)方法[4]計(jì)算獲得重力及其梯度異常。圖4所示為選取256×256全空間計(jì)算所得準(zhǔn)確的地下異常體正演結(jié)果，其與地下異常體一一對(duì)應(yīng)。采用32×32、24×24及16×16的子空間分別計(jì)算，獲得重力及其梯度異常的結(jié)果分別如圖5～圖7所示，各子空間的計(jì)算時(shí)間如表1所示，以及全空間與各子空間的計(jì)算時(shí)間比如圖8所示。

由圖4～圖8可以看出：①隨著子空間范圍的縮小，運(yùn)算時(shí)間隨之縮短，大大提高了計(jì)算效率；②子空間計(jì)算范圍縮小，計(jì)算精度下降。

為了評(píng)價(jià)本文所提方法的精度，將子空間為32×32的計(jì)算結(jié)果(圖5)與理論值(圖4)相減，其結(jié)果如圖9所示?？梢钥闯稣`差值基本上在零值附近。

為了定量評(píng)價(jià)誤差的大小，文中統(tǒng)計(jì)了重力及其梯度異常的最大、最小值，以及計(jì)算結(jié)果與理論值誤差的均值和均方差(如表2所示)。均值公式為：

圖4 256×256全空間運(yùn)行結(jié)果Fig.4 256×256 full space operation result

圖5 32×32子空間運(yùn)行結(jié)果Fig.5 32×32 subspace operation result

圖6 24×24子空間運(yùn)行結(jié)果Fig.6 24×24 subspace operation result

圖7 16×16子空間運(yùn)行結(jié)果Fig.7 16×16 subspace operation result

表1 256×256全空間不同子空間運(yùn)行時(shí)間

圖8 全/子空間運(yùn)算時(shí)間比值Fig.8 Full/subspace operation time ratio

圖9 256×256全空間與選取32×32子空間計(jì)算誤差Fig.9 256×256 full space and selected calculated 32×32 subspace error

(8)

式中：a1,a2,…,an為計(jì)算結(jié)果與理論值的誤差矩陣元素；n為矩陣所包含的元素的數(shù)量。

均方差的公式為：

(9)

由表2可得，重力異常計(jì)算值與理論值誤差的均值與均方差分別為1.326 5 g.u.、0.717 5 g.u.，相比其理論最大值(72.127 g.u.)、最小值(0.089 6 g.u.)的范圍，誤差相對(duì)較小。

為了進(jìn)一步定量衡量計(jì)算結(jié)果的精度，文中同時(shí)也統(tǒng)計(jì)了子空間32×32計(jì)算結(jié)果的均方差和平均相對(duì)誤差(如表3所示)，計(jì)算結(jié)果的均方差公式為：

(10)

式中：xi為子空間各觀測(cè)點(diǎn)計(jì)算結(jié)果；(x*)i為全空間各觀測(cè)點(diǎn)計(jì)算結(jié)果；n為觀測(cè)點(diǎn)數(shù)量。理論值與其平均值的均方偏差為：

表2 全空間重力及其梯度異常最大、最小值及計(jì)算值與理論值誤差的均值和均方差

(11)

(12)

表3 32×32子空間計(jì)算結(jié)果的均方差和平均相對(duì)誤差

由表3可得：32×32子空間計(jì)算所得結(jié)果中重力異常和Uxz、Uyz兩個(gè)梯度異常所得結(jié)果的精度較高。重力異常的平均相對(duì)誤差為9.01%，Uxz、Uyz的平均相對(duì)誤差數(shù)值在5%以下，另Uzz、Uxx、Uyy的平均相對(duì)誤差數(shù)值在20%以下，Uxy的平均相對(duì)誤差數(shù)值為29.09%。通過(guò)表3，并結(jié)合表2中理論最大值和最小值，可以看出，重力及部分重力梯度異常計(jì)算結(jié)果誤差較小。

4 反演結(jié)果對(duì)比

有限內(nèi)存擬牛頓方法已被證明在解決重力及其梯度正演中具有一定的優(yōu)勢(shì)性[8]。因此應(yīng)用有限內(nèi)存擬牛頓方法對(duì)Moving-footprint方法的計(jì)算效果進(jìn)行檢驗(yàn)。為驗(yàn)證本文所提正演方法在計(jì)算中的優(yōu)勢(shì)，采用文中方法與文獻(xiàn)[4]所提正演方法在反演中的效果進(jìn)行對(duì)比。首先選擇較小空間的數(shù)據(jù)集進(jìn)行運(yùn)算。全空間網(wǎng)格剖分為16×16×9，子空間為8×8×9。模型為一長(zhǎng)方體模型，大小為500 m×500 m×500 m，頂面埋深為500 m，剩余密度為1 g/cm3。分別應(yīng)用現(xiàn)有計(jì)算方法與本文所提Moving-footprint方法得出的計(jì)算結(jié)果如圖10所示。

正演方法和本文方法這兩種計(jì)算方式所需的時(shí)間分別為：29.222 s、26.602 s。由此可知應(yīng)用Moving-footprint方法可以降低運(yùn)算時(shí)間。當(dāng)選擇32×32×9的全空間進(jìn)行計(jì)算時(shí)，所得結(jié)果如圖11所示，現(xiàn)有方法計(jì)算所需時(shí)間為1 543.53 s，而在應(yīng)用Moving-footprint方法后，計(jì)算時(shí)間為703.25 s，因此應(yīng)用Moving-footprint技術(shù)對(duì)計(jì)算效率進(jìn)行了提升，且隨著模型空間剖分?jǐn)?shù)量的增加，反演的計(jì)算效率有顯著提升。

a—原始方法；b—Moving-footprint方法a—the result obtained by the original method；b—the result obtained by applying the Moving-footprint method圖10 全空間16×16×9不同方法反演結(jié)果Fig.10 Inversion results of different methods in full space 16×16×9

當(dāng)選擇更大模型剖分空間進(jìn)行反演時(shí)，如256×256×9，嵌套已有的正演算法無(wú)法完成迭代過(guò)程。而嵌套本文所提Moving-footprint技術(shù)可以有效完成對(duì)大尺度模型的反演計(jì)算，迭代1次所需時(shí)間約為30 min。

a—原始方法；b—Moving-footprint方法a—the result obtained by the original method；b—the result obtained by applying the Moving-footprint method圖11 全空間32×32×9不同方法反演結(jié)果Fig.11 Inversion results of different methods in full space 32×32×9

5 結(jié)論與建議

本文借鑒航空電磁正演計(jì)算的“Moving-footprint”技術(shù)，提出了基于“Moving-footprint”重力及其梯度異常的正演計(jì)算方法。該方法選擇全空間范圍內(nèi)的一定子空間，只計(jì)算存儲(chǔ)在子空間范圍內(nèi)的網(wǎng)格點(diǎn)的格架函數(shù)，即只考慮子空間范圍內(nèi)的計(jì)算點(diǎn)在觀測(cè)點(diǎn)產(chǎn)生的重力及其梯度異常。在保證一定的計(jì)算精度的同時(shí)，縮短了運(yùn)算時(shí)間，提升了計(jì)算效率。

海量重力數(shù)據(jù)線性迭代反演時(shí)，初步迭代應(yīng)用“Moving-footprint”技術(shù)選取合適的子空間大小，提升了計(jì)算效率。從而使得大范圍全空間重力及其梯度異常正反演計(jì)算成為可能。