基于空間域迭代的Parker-Oldenburg密度界面反演算法研究

［關鍵詞］Parker-Oldenburg密度界面反演；空間域；頻率域；Matlab；迭代

密度分界面是很多油氣藏、礦藏等的賦存區域[1-3]，通過重力密度界面反演計算密度分界面的深度有利于研究密度分界面的形狀起伏以及分布[4-6]，對于區域地質構造的研究和石油天然氣的勘探具有重要意義。

根據重力密度界面反演方法的不同，可以歸結為兩大類，一類是空間域的密度界面反演方法，一類是頻率域的密度界面反演方法[1] [4-6]。在空間域反演方法中，有直接反演方法和間接迭代方法。直接反演方法是依據場源的物理和幾何參數先建立一個外部位場函數關系，在這個外部位場函數關系的基礎上，實現位場的反演，得到地下場源的物理和幾何參數，稱之為線性反演法；間接迭代方法則主要通過空間域迭代的途徑，不斷修改反演模型的初值，以逐漸逼近實測異常值，而最終達到與實測重力異常最佳擬合的方法[7-10]。

由于空間域的重力密度界面反演運用的是迭代思想，即通過設計正演模型計算出模型的重力場，再與實際測量的重力場進行比較，通過修改模型參數一步一步逼近實際測量值，在滿足精度要求后則近似認為模型參數就是所求地質體的參數。因此，空間域反演法的精度并不是很高，且由于很大一部分反演所耗時間都用在正演計算模型體的重力異常上，并且這些模型體一般都是單一模型。但實際情況中都需要采用組合模型，以滿足實際的需要，對每個異常點都要進行一次計算并進行迭代，耗時巨大，反演速度也就很慢。因此，頻率域反演方法應運而生，并逐漸成為了目前最流行的密度界面反演方法。但是兩種方法都存在各自的缺陷，頻率域密度界面反演方法需要預設界面平均深度，為保證反演結果的收斂還需要添加一個低通濾波器，而空間域反演方法則存在計算效率不高、精度較低等缺點。基于此，本文提出一種基于空間域迭代的Parker-Oldenburg 密度界面反演算法[7] [8]，反演過程中不僅不需要設置界面平均深度和低通濾波，而且進行頻率域正演具有更快的速度，在密度界面反演過程中擁有更高的效率和精度。

1. 研究內容與研究思路

研究內容分為三個部分，一是基于Granser的以FFT算法為基礎的正演方法[9]，該方法與Oldenburg反演方法相似[8]，主要解決計算每個矩形棱柱（或多邊形的每邊）的重力異常相加得到總重力場時計算時間很長的問題，用以從重力異常中確定密度界面的幾何參數，盡管基于FFT的算法減少了計算時間，但它需要給定密度界面的平均深度和低通濾波器來實現反演過程的收斂；二是基于Cordell和Henderson的空間域迭代反演方法[9] [10-12]，在三維反演中，雖然該方法使用堆疊矩形棱鏡模型計算速度較慢，但該方法是基于布格平板公式，假設每個棱鏡在所有水平維度上都是無限的平板，從而對每個棱鏡的底部深度做出初步猜測；三是基于空間域迭代的Parker-Oldenburg密度界面反演，為了綜合兩者的優勢，解決它們各自的局限性，考慮將Granser方法和Cordell和Henderson方法結合在一個模型中[7-10]。該算法在Matlab環境下實現，通過對本文理論方法的程序編寫（圖1）。

將Parker-Oldenburg密度界面反演方法程序化之后，導入需要進行反演的實測重力異常網格數據文件，輸入初始模型，再輸入在地表觀測到的相對基底巖層的密度差和密度隨深度指數衰減的衰減常數，最后選擇迭代停止條件，完成以上設置之后開始進行反演。反演過程對初始模型進行頻率域正演，然后計算觀測異常與計算異常之間的均方根誤差RMS，再計算更新模型，判斷均方根誤差是否滿足迭代停止條件，不滿足則需要更新模型并重新正演計算均方根差，若滿足迭代停止條件則迭代停止，輸出反演結果并得出反演計算時間。

2. 密度界面反演基本原理

密度界面反演是根據地面實測的重力異常Δg 來確定地下密度分界面的起伏。由于重力異常的復雜性和反演解釋中的非單一性，要使密度界面反演取得良好的效果，必須具備一定的條件，即用來反演計算的重力異常必須是單純由密度界面的起伏引起；界面上下物質層的密度分布比較均勻、穩定且已知它們的密度差；有一個或幾個點的密度界面深度為已知。

首先，建立一個測點與界面異常源位置關系的直角坐標系，坐標原點O即為測點位置，x-y平面在水平面上，z軸沿重力方向垂直向下。界面以上地層密度為ρ₁，界面以下為基底密度為ρ₂（圖2）。

密度界面反演過程中需要由已知模型根據地下場源在地面產生重力異常的計算公式，計算出各測點的重力異常值[1-10]。密度分界面用深度函數h（x， y ）來表示，假設基底界面h（x，y ）的平均深度為z₀，則界面深度可使用界面相對于平均深度的起伏Δh（x，y ）函數來表示，并假設密度分界面以上地層的密度為ρ₁，密度分界面以下基底密度為ρ₂（圖3），界面上下的密度差為Δρ （x，y ），如果界面以上地層密度是變化的，即隨著深度的變化同一測點下方的密度差Δρ （x，y ）也是一個變化函數。為了方便計算，一般假設相對界面的起伏總在界面的一側，且起伏不大。

2.1 頻率域密度界面反演

Parker-Oldenburg公式

最早的密度界面反演方法大多是在空間域中進行的，自Parker在重力密度界面正演計算中使用快速傅里葉變換（FFT），在這之后由他的學生Oldenburg根據Parker公式，提出了在頻率域中的密度界面迭代反演方法Oldenburg。由于頻率域計算的速度遠遠快于空間域的計算，此后，頻率域密度界面反演便成為應用最為廣泛的密度界面反演方法[7-10] [13]。

采用基于Granser的頻率域正演公式計算出沉積盆地的正演重力異常（圖6d），異常以134列×106行網格數據（網格間距為1 km）繪制盆地的界面深度結構圖，在盆地中央是較大的負異常，盆地邊緣異常接近零（圖6c）。

反演部分主要是利用模型的正演重力異常，再根據模型參數，設置反演算法所需的實際網格間距，初始密度差，密度衰減常數以及迭代停止條件，迭代停止條件有兩種不同的方式，一種是當下一次迭代的觀測數據和計算數據之間的均方根誤差（RMS）相對于上一次迭代的RMS增加時，迭代停止；第二種是當迭代的均方根誤差（RMS）低于一個預設的閾值時迭代停止。

此處我們選擇第一種方式，之后開始反演。如果由于邊緣效應的影響導致反演結果不理想，可以選擇忽略幾個邊緣數據。反演結束之后得出反演時間為2.59s，迭代進行了49次達到了設置的最大迭代次數，而不是達到了迭代停止條件，若49次迭代之后的RMS值未能達到反演要求，之后可以設置更大的最大迭代次數。

由于實測數據有測網以外重力的影響，而在實際的模型正反演計算中使用的是有限模型，在反演計算中測網邊緣的測點的重力異常只能通過有限模型的邊緣和中間模型提供，因此在邊緣會出現與實際情況存在較大差異的值以盡量符合實測數據。為了去除邊緣效應的影響，可以選擇忽略邊緣測點的數據，使實際反演區域小于實測區域。如果邊緣效應對結果影響不大也可不選擇忽略邊緣數據。在此模型中邊緣效應影響較小，因此選擇忽略邊緣為0，得出反演結果之后可以導出反演結果的網格數據，并對比模型正演異常與反演界面深度正演異常，來比較反演結果是否可靠。

用于反演的重力異常，是觀測異常通過反演之后得出基底界面深度，然后利用該反演出的基底界面深度再重新計算其在各觀測點上的異常（圖7a），可以看出該異常（圖7a）與觀測異常（圖6d）幾乎完全相同（圖8a）。對重力異常進行反演，基底界面的最大深度為7.5 km左右，最小深度幾百米，反演深度（圖7b）與模型實際深度（圖6c）誤差大面積接近0 km，整體誤差在-0.3 km到0.5 km之間（圖8b）。

本次試驗選擇的迭代停止條件為第一種，即當下一次迭代的均方根差（RMS）大于上一次迭代的均方根差時迭代停止，在迭代開始前需要設置一個最大迭代次數，當最大迭代次數設置為50時，迭代結果一直沒能滿足第一種迭代停止條件，達到了最大迭代次數，在完成第49次迭代后停止，從這49次迭代計算的RMS值來看，均方根差RMS的收斂很快，當達到20次左右時迭代計算的RMS值就已經很小了，最后一次迭代的RMS值已經達到了0.0020475（圖8c）。

綜上，觀測異常與計算異常差值的變化范圍不大，在-0.013～0.012 mGal之間，且靠近沉積盆地邊緣和沉積盆地中心處觀測異常與計算異常的偏差較大一些。在邊緣有三處異常差明顯偏大的區域，與反演界面深度圖對比可以看出，這三處為邊緣較低處。這種在邊緣較低處產生較大正誤差的原因可能受邊緣效應的影響。為了證明這個觀點繼續設計了其他幾組背斜模型和復雜模型。從算法反演結果中能得出反演計算的時間，選擇第一種迭代停止條件反演進行到49次的反演時間為2.59s，選擇第二種迭代停止條件反演進行到迭代的第25次停止，反演計算時間為1.37s，由于在迭代進行的過程中，如果一直沒有達到停止條件而持續迭代有可能出現過度擬合的情況，即會將一些微小起伏的異常放大，反而使結果不可靠。因此需要根據模型擬合情況選擇適當的迭代停止條件，以使反演結果最佳，比如在此次模型中選擇第二種迭代停止條件，也就是迭代25次之后的結果就已經相當可靠了，就不需要選擇第一種迭代停止條件讓它繼續迭代下去了。由此得出，該算法在向斜盆地模型中擁有良好的反演效果。

4. 結論

（1）基于空間域迭代的Parker-Oldenburg 密度界面反演方法具有反演速度快，反演精度較高的優點，結合了頻率域正演方法和空間域迭代反演的優點，不需要給定界面的平均深度，也不需要引入低通濾波來控制迭代的收斂，擺脫了單獨的空間域反演精度不高、反演速度慢和單獨頻率域反演需要計算界面平均深度和加入低通濾波器的不足。

（2）通過盆地模型試驗，該反演算法可以完成地質構造界面的反演，并且對于密度隨深度變化的地層同樣具有較好的效果。

（3）在模型試驗中，采用指數變密度函數時，由于邊緣效應的影響，可能需要忽略更多的邊緣測點數據，以達到理想的反演效果。