基于起伏地層的曲面重力場快速高精度正演

鄭翾宇， 柳建新， 陳 波， 趙廣東， 陳龍偉， 郭榮文

(1.中南大學 地球科學與信息物理學院，長沙 410083;2.有色資源與地質災害探查湖南省重點實驗室，長沙 410083)

0 引言

地形模型的正演是重力數據定量解釋的基礎。正演問題中的地形改正和等效補償，反演問題中所涉及的重要密度界面(如沉積層底部界面、地殼與地幔的分界面等)，都需要一個快速而準確的算法來計算出地形模型產生的重力場。地形模型重力場的正演方法可分為空間域和頻率域兩類[1]。空間域方法是將地形模型剖分成多個便于計算的規則單元體，分別計算每個單元體在觀測點處產生的重力異常，最后累加求和得到該地形模型總的重力異常。這些規則單元體產生的重力場具有解析表達，通過解析式可以精確地計算出空間任意一點處的重力異常場。許多學者在空間域正演方法方面做了大量工作，Hubbert[2]用線積分的形式來計算任意二度體的重力場，同時也采用了圖形法來計算重力場；Nagy[3]給出了長方體模型的重力場解析式，同時指出復雜形體可以用長方體作為基本單元疊加來模擬，為重力場的地形校正提供了基本的思路；Paul[4]給出了以三角形為面元的任意三度體的重力場；Plouff[5]將Talwani[6]的薄板公式推廣到邊界為任意多邊形的直立棱柱體的異常解析式，同時利用該解析式做地形校正；Li[7]回顧和整理了直立長方體、上下底面為任意多邊形的直立棱柱體以及任意多面體模型的重力場及其高階導數在全三維空間的正演公式，并初步研究了重力場及磁場在場源邊界處的連續性問題；Nagy[8]詳細討論了直立長方體模型的重力位及其三階以下導數在全空間的連續性問題，詳細推導了各奇點處的解析積分式，并分析了奇點處的極限性質；Garcia-Abd eslem[9]推導了密度隨深度呈三次多項式變化的直立長方體的重力場解析公式，并用之前的解析和數值相結合的方法[10]對結果進行了驗證。這些研究極大地豐富和完善了空間域重力正演方法，但是空間域方法的解析式往往較為復雜，推導過程比較繁瑣，尤其當場源體模型較為復雜，觀測點數較多時，由于計算次數會隨著模型剖分個數以及觀測點數急劇地增加，計算過程將會特別耗時，計算效率成為制約空間域方法運用的瓶頸。

相比于空間域方法，頻率域方法引入了快速傅里葉變換(FFT)，計算速度得到顯著提升。同時，在異常頻譜的表達式中，場源體的幾何特征表現為簡單的乘積因子，異常頻譜的表達式更為簡潔和緊湊，頻率域方法也因其具有的獨特優勢而被廣泛采用。相鵬等[11]采用頻率域方法研究了磁性界面；馮娟等[12]用頻率域方法研究了三維密度界面和重力位場之間的關系；姜永濤等[13]采用了頻率域算法研究了中國西部地區的莫霍面深度；盧鵬羽等[14]基于重力全張量數據采用頻率域算法研究了密度界面。

頻率域算法的發展從20世紀60年代開始，Bhattacharyya[15]給出了長方體磁化的磁場頻率域表達式，并指出在頻譜表達式中場源形體和物性等參數是彼此分離的乘積因子，從而場源體的物性及尺寸、埋深等各參數可以直接利用異常場的頻譜來分析；Schouten[16-17]采用FFT算法計算了層狀體的磁場異常?；谶@一思路，Parker[18]將快速傅里葉變換引入重力位場正演，提出了用于地形校正的Parker公式，該方法是從頻率域計算均一密度的長方體疊加地形模型的位場。在利用傳統的Parker公式正演起伏地層產生的重力異常時，往往存在一定的誤差。一是傅里葉變換時產生的；其次Parker公式本身存在一定的缺陷，這一缺陷影響正演結果的收斂速度，導致正演誤差的產生。Forsberg[19]在詳細研究了長方體疊加的場源體模型的重力場之后，將Parker公式中級數展開的思想改進為空間域級數展開法，再利用FFT算法計算一系列褶積的和；馮銳[20]將Parker公式的密度由原來的橫向變化推廣至沿深度方向按指數或線性變化的情況；Xia等[21]提出了通過增加上下界面平均值來提高Parker公式的數值穩定性加速級數收斂的思想；Chai[22]詳細介紹了偏移抽樣方法，該方法大幅提高了頻率域重力正演精度；柴玉璞[23]提出了變衰減系數指數和變系數多項式密度模型的Parker公式，為了提高傅里葉算法的數值精度文中采用乘子法[24]和偏移抽樣方法；柴玉璞[25]將偏移抽樣方法進行了重新的疏導和整理并發展為一套完整的理論，將離散傅里葉變換的精度提高數十倍，減小了基于傅里葉變換的頻率域位場正演誤差；Wu等[26]指出，采用標準FFT方法正演產生不同程度畸變的原因是離散傅里葉變換這種算法(梯形積分)不能很好地逼近傅里葉這一震蕩積分，從這個角度出發引入偏移參數并提出Gauss-FFT算法，該算法不僅大幅減小了頻率域計算時“邊界效應”、“混疊效應”、“詭源效應”等問題，同時還極大地提高了頻率域方法的計算精度。這些頻率域方法的提出以及傅里葉算法的改進，為頻率域方法提供了廣闊的應用空間。

筆者從重力位場出發對Parker公式進行了重新的推導，參考Xia等[21]的思路，引入上下界面平均值，將界面的絕對起伏值轉化為相對平均界面的起伏量，加快了正演的收斂速度，提高了正演公式的數值穩定性?？紤]到實際情況下的觀測面并非水平，甚至有時我們需要一個曲面上的重力異常值，對此提出了一個曲面觀測算法，將傳統的水平觀測面拓展到了起伏觀測面，實現了上下界面起伏地層在起伏觀測面上的快速高精度重力正演計算，并采用Gauss-FFT算法進一步保障正演的數值精度。

1 方法原理

1.1 正演理論方法

自Parker將傅里葉變換引入到重力場正演計算并推導出著名的Parker公式以來，這種以快速高效著稱的頻率域算法已被廣泛應用于地形校正和其他正演計算，但是傳統的Parker公式存在一定的缺陷，當界面起伏較大時不能保證正演結果的數值穩定。我們對Parker公式進行了重新的推導，并引入上下界面平均值，改進后的Parker公式具有更高的數值穩定性。

如圖1所示的直角坐標系中，x-y平面置于水平，z軸豎直向上，與重力方向一致。重力異常的場源體是由上下起伏面h1(ξ,η)和h2(ξ,η)所圍成的地層模型，場源體的坐標為(ξ,η,ζ)。觀測面為z=z0的平面， 觀測點坐標為(x,y,z0)。

圖1 起伏地層模型和觀測面示意圖Fig.1 Complicated undulating stratigraphic model and observational schematic diagram

場源體在觀測面上某一點處產生的重力位U(x,y,z0)(簡記為U)可表示為式(1)[27]。

(1)

(2)

式(2)中：Kx、Ky分別表示x、y方向上的波數。將式(2)改變積分順序，根據傅里葉平移性質(3)和Bracewell[28]推導的變換式(4),我們可以將式(2)經變換得到式(5)。

f(x-ξ)?F(k)e-ikξ

(3)

(4)

(5)

將地層的上、下面h1(ξ,η)和h2(ξ,η)代入式(5)并對ζ求積分可得式(6)。

(6)

對重力位求z向偏導，得到上下界面起伏地層模型在水平觀測面上產生的重力異常的頻率域表達式為式(7)。

F[Δg]=-2πG

(7)

式(7)即為快速收斂的起伏地層正演公式，與傳統的Parker公式(8)相比具有更高的數值穩定性。

(8)

采用傳統的傅里葉變換計算時，正演結果會出現不同程度的畸變。為了避免這一問題，均采用4節點的二維Gauss-FFT算法[22]。

1.2 曲面觀測

在前人的研究工作中，正演選取的觀測面均為水平觀測面。受到地形等外界條件的影響，實際情況下的觀測面并非水平面，如野外重力數據采集時，觀測面是隨地形變化的起伏面。因此，快速且精確地正演出場源體在起伏觀測面上的重力異常，是一項具有實際意義的研究工作。

如圖2 (a)所示，場源體是上、下起伏面h1(ξ,η)和h2(ξ,η)之間的地層，觀測面為z=z(x,y)的起伏面。為了計算起伏觀測面上的重力異常，采用泰勒級數展開法，依據泰勒級數展開理論，起伏面上的重力場可以用式(9)展開。

圖2 示意圖Fig.2 Sketch map(a)曲面觀測示意圖；(b)泰勒級數展開法示意圖

g(x,y,z(x,y))=g(x,y,z0)+

(9)

其中：g(x,y,z0)表示平面z=z0上的位場；Δz(x,y)=z(x,y)-z0(圖2(b))。

由式(9)可知，利用泰勒級數法的基本思想，可以根據平面z=z0上的重力位場g(x,y,z0)及其對應的各階導數值，來逼近曲面z=z(x,y)上的重力位場值g(x,y,z(x,y))。依據傅里葉變換的性質，各階導數可以按照式(10)轉換到頻域計算[23]。

式中:FT[]表示傅里葉變換；IFT[]表示傅里葉反變換。

根據起伏地層在水平觀測面上產生的重力異常正演式，可以直接求出起伏面z=z(x,y)對應的泰勒展開面z=z0上的重力異常f(x,y,z0)。根據泰勒展開式(9)和傅里葉變換中導數的轉換性質(10)，可得到起伏觀測面z=z(x,y)上的重力異常為式(10)。

(10)

起伏觀測面上重力異常的精度受泰勒展開面z0的取值和泰勒級數展開項數的影響，一般來說z0取起伏面的平均值最佳，泰勒級數展開項數越多正演結果越精確，但計算量也隨之增加。

2 模型算例

建立簡單的模型體，從模型實驗數據出發檢驗正演算法的可靠性。由于空間域算法是將場源體剖分成單個的直立棱柱體，而每個棱柱體在觀測點上產生的重力異?？梢酝ㄟ^解析式準確計算，因此我們將空間域正演結果作為標準來衡量正演算法的精度。

2.1 起伏地層的正演計算

如圖3(a)所示，在z軸豎直向上的坐標系中，產生重力異常的場源體上下界面為凹凸起伏面，頂、底界面的平均值分別為l1=-1 km,l2=-3 km；模型在x方向和y方向均為256 km,即0≤ξ≤ 256 km,0≤η≤ 256 km；場源體模型的密度為ρ=800 kg/m3；觀測面為z=1.5 km的水平面，x方向和y方向上的觀測點數Nx=Ny=256，網格間距dx=dy=1 km。

首先采用Nagy[3]提出的空間域方法，計算出場源體在觀測面上產生的重力異常，如圖3(b)所示。空間域正演結果的平均值為-61.80 mGal，標準差為8.09 mGal?？傮w來說，重力異常數值在中間較大往四周逐漸減小，在邊界處取得最小值-19.89 mGal。

采用Parker[14]提出的上下界面起伏地層的正演公式和改進的正演公式(7)分別對圖3(a)所示的起伏地層做正演計算(圖4)。

圖4(a)所示的正演結果和空間域結果相比大體一致但在數值上略有偏差，圖4(b) 直觀的給出了采用傳統Parker公式正演的誤差大小和分布規律，由圖4(a) 、圖4(b)可知：傳統Parker公式的正演誤差主要集中分布在起伏地層的邊緣處，對于地層沒有起伏變化的區域，正演誤差相對較小?？傮w而言，絕對誤差數值范圍從-18 mGal至26 mGal，最大相對誤差約為10%。

圖4(c)為改進的Parker公式正演結果，平均重力異常為-61.81 mGal，標準差為8.08 mGal。圖4(c)的正演結果在數值上與空間域方法的正演結果幾乎一致，數值的變化和分布規律幾乎相同。

圖4 正演結果和誤差Fig.4 Forward results and error maps(a)傳統Parker公式正演結果；(b)傳統Parker公式正演誤差；(c)改進的Parker公式正演結果；(d)改進的Parker公式正演誤差

由圖4(d)可知，改進后的Parker公式最大相對誤差約為0.1%。

該模型算例表明，改進后的Parker公式(7)在起伏地層的正演計算中具有更高的數值精度，而正演數值精度提高的關鍵點在于增加上、下界面平均值這兩個參數后，將界面的絕對深度轉化為相對界面平均值的起伏量，加速了正演結果地收斂。

圖5 起伏觀測面和正演結果Fig.5 Undulating observation surface and forward results(a)觀測面示意圖；(b)空間域方法正演結果；(c)本文提出的頻率域方法正演結果；(d)新頻率域方法正演誤差

2.2 起伏地層的曲面觀測

場源體模型與上面一致,不同之處在于觀測面由原來z=1.5 km的水平面變為如圖5(a)所示的起伏面，觀測點數和網格間距保持不變，x方向和y方向上的觀測點數Nx=Ny=256，網格間距dx=dy=1 km。觀測面是相對z=1.5 km水平面存在正弦函數形凸起和凹陷的起伏面，凸起區域是以(64,128)為中心64 km為半徑的圓形區域，凹陷區域則是以(192,128)為中心64 km為半徑的圓形區域。

圖5(b)是空間域方法正演得到的起伏觀測面上的重力異常。重力異常的平均值為-61.86 mGal,最大和最小的重力異常值分別為-83.34 mGal和-19.89 mGal,標準差為8.27 mGal。

圖5(c)是曲面正演理論，從頻率域正演得到的起伏觀測面上的重力異常，泰勒級數的展開面為z=1.5 km的水平面，展開項數為10項?？傮w來看，由圖5(c)可知，其正演結果與空間域方法計算結果5(b)幾乎一致，在數值方面，平均重力異常為-61.85 mGal，最大和最小的重力異常值分別為-82.83 mGal和-19.96 mGal,標準差為8.24 mGal。

圖5(d)為頻率域方法的正演誤差。從誤差圖可知，頻率域正演具有較高的數值精度，平均相對誤差小于百分之一。

空間域方法的正演時間為2 165 s，頻率域方法的正演時間僅需16 s。測試計算機配置為CPU-Intel(R) Core(TM) i5-4590，主頻為3.30 GHz，內存為8 GB。由此可見該方法在保障正演精度的同時也兼具很高的計算效率。

通過大量的模擬實驗發現，曲面觀測的正演計算受泰勒級數展開平面和級數展開項數的影響。當泰勒展開平面取起伏面的平均值時正演精度最高，離平均值越遠精度越低。泰勒展開項數則會影響正演精度和效率，正演精度隨著泰勒展開項數的增加而提高，計算效率則相反。為了取得較高的正演精度，同時保留頻率域算法的效率優勢，泰勒級數的展開項數一般取6項到10項。

圖6 帶地形的密度界面起伏模型Fig.6 Density interface with undulating terrain model

2.3 起伏地形的密度界面模型正演模擬

上述模型算例均是計算起伏地層在某一觀測面上產生的重力異常，而Parker公式除了用于地形校正以外，也常用于密度界面的研究。圖6為帶地形的密度界面起伏模型，觀測面為地表起伏面z=z(x,y)，數值與圖5(a)一致。密度界面為ζ=h1(ξ,η)，在(98,128)和(158,128)兩點處分別存在一個以30 km為半徑，幅值為3 km的正弦形凸起和凹陷。密度界面的平均值h0=-5 km,密度差ρ2-ρ1=400 kg/m3。

圖7 正演結果圖Fig.7 Forward results(a)空間域正演結果；(b)頻率域正演結果

計算起伏密度界面在水平觀測面上產生的重力異常時，公式 (7) 被簡化為式(12)。

F[(h1(ξ,η)-h0)n·ρ(ξ,η)]

(12)

根據式(12)計算出起伏密度界面在z=1.5 km平面上產生的重力異常，計算出地表觀測面z=z(x,y)上的重力異常，正演結果如圖7所示。

圖7(a) 為空間域方法正演結果，圖7(b)為頻率域方法正演結果，兩個正演結果的等值線圖幾乎一致。為了進一步研究頻率域方法的正演精度和計算效率，將兩種方法的正演結果進行了統計分析，整理的結果如表1所示。

表1 空間域方法和頻率域曲面觀測方法正演對比Tab.1 Forward comparison of spatial and requency domain methods

從表1可知，本文提出的曲面觀測面算法在正演圖6所示的起伏密度界面模型時，具有非常高的數值精度，在計算效率方面相對于空間域正演方法更是提高了一百多倍。

3 結論

首先對起伏地層的重力正演方法做了簡單的綜述，總結了空間域和頻率域兩種方法各自的優點與不足，并對方法的發展情況作了簡單的介紹。隨后針對頻率域Parker公式存在的問題展開討論，從重力位場出發對Parker公式進行了詳細的推導，并加入上、下界面平均值兩個參數對傳統的正演公式加以改進。在此基礎上，提出了一種曲面觀測算法，實現了起伏地層起伏觀測面的快速高精度重力正演計算和帶地形密度界面起伏模型的快速高精度正演。從模型實驗可知：

1)改進后的Parker公式在起伏地層的正演計算中具有更高的數值精度，其核心在于將界面的絕對深度值化為相對起伏值，大幅減小了由于數值不穩定而產生的正演誤差。改進后的Parker公式可用于地形校正，沉積層校正以及密度界面的正反演等方面的研究。

2)本文提出的曲面觀測的正演方法，在保證正演精度的同時計算效率相對傳統空間域方法高出兩個數量級，是一種計算起伏觀測面上重力異常場的有效方法。