重力及其梯度異常正演的Moving-footprint大尺度模型分解方法

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

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

0 引言

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

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

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

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

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

1 重力及其梯度異常正演

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

圖1 地下半空間單元劃分Fig.1 Underground half-space unit division

(1)

(2)

(3)

(4)

(5)

(6)

(7)

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

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

注：紅線包圍區(qū)域為子空間note：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)格點幾何格架函數(shù)技術(shù)的基礎(chǔ)上進行了改進。因此，首先在全空間內(nèi)選取子空間的大?。蝗缓笥嬎阕涌臻g內(nèi)部分測點的網(wǎng)格點的幾何格架函數(shù)，并存儲以備調(diào)用；隨后計算某觀測點異常時，先判斷該觀測點所在的子空間，再調(diào)用網(wǎng)格點格架函數(shù)和該子空間單元體的剩余密度，求和后獲得該觀測點異常；最后，利用Moving-footprint完成觀測區(qū)域所有點的重力及其梯度異常正演，主要計算步驟：

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

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

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

4)觀測點移動，子空間隨之移動，利用步驟3完成移動觀測點重力及其梯度異常的正演計算；最終完成所有觀測點的正演計算。

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

3 模型實驗

為了檢驗本文所提方法的計算效果，將地下半空間剖分為256×256×15個長方體單元，長方體單元大小為0.1 km×0.1 km×0.1 km。采用4個長方體組合模型進行檢驗，長方體組合模型的剩余密度都為0.5 g/cm3，長方體組合模型大小為4.0 km×4.0 km×0.5 km，其中心點坐標分別為(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所示，其中，正演計算點距為0.1 km×0.1 km。

注：紅線包圍區(qū)域為計算子空間，藍色立方體區(qū)域為異常體區(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

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

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

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

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

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

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

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

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

表1 256×256全空間不同子空間運行時間

圖8 全/子空間運算時間比值Fig.8 Full/subspace operation time ratio

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

(8)

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

均方差的公式為：

(9)

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

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

(10)

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

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

(11)

(12)

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

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

4 反演結(jié)果對比

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

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

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

當選擇更大模型剖分空間進行反演時，如256×256×9，嵌套已有的正演算法無法完成迭代過程。而嵌套本文所提Moving-footprint技術(shù)可以有效完成對大尺度模型的反演計算，迭代1次所需時間約為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é)論與建議

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

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