基于雙曲線密度模型的頻率域界面反演

商宇航 邰振華 秦 濤

(黑龍江科技大學礦業工程學院，黑龍江哈爾濱 150022)

1 引言

在深部構造和油氣勘探等研究領域，密度界面反演是位場數據處理解釋的重要方法之一[1,2]。該方法以界面兩側地質體的密度差為反演參數，采用線性或者非線性的模式獲取地下界面深度[3]。Athy[4]提出盆地基底之上巖石與基底的密度差隨深度增加而減小(這與沉積盆地內蓋層的巖石密度隨深度增加而變大是相關的)； Cordell[5]也認為盆地內巖石與基底的密度差隨深度增加而減小，且淺部的速度降低往往快于深部，于是提出指數密度模型，并采用遞推算法反演了加利福尼亞San Jacinto地塹的基底深度； Murthy等[6]、Rao[7]和陳勝早[8]分別提出了線性密度和二項式密度模型，在盆地重力異常擬合反演中取得了較好的效果。目前，這幾種變密度模型及反演模式已得到推廣和應用[9,10]。

同樣在San Jacinto地塹的研究中，Litinsky[11]發現雙曲線密度模型能更好地擬合密度—深度數據； Rao等[12]在雙曲線密度模型研究的基礎上，給出了空間域密度界面重力異常的擬合公式，在San Jacinto地塹和印度的Godavari地塹的重力剖面擬合中得到了更好的效果。但由于計算過程需要在空間域實現，計算較復雜，精度也偏低。隨著快速傅里葉變換(FFT)的提出，密度界面正、反演被引入頻率域計算，運算速度得到大幅提高。Parker[13]率先推導了頻率域單一物性界面的位場異常正演公式； Oldenburg[14]在此基礎上，提出了常密度界面的迭代反演方法(PO法)，使頻率域密度界面的正、反演融為一體，成為常規算法[15,16]； 柯小平等[10]提出了基于指數密度模型的頻率域界面反演方法，在青藏高原莫霍面深度反演中取得了良好效果。

本文在前人工作基礎上，在頻率域推導了基于雙曲線密度模型的界面正、反演公式，并利用模型試驗和實際重力場數據資料驗證該方法的有效性和實用性。

2 基本原理

2.1 基于雙曲線密度模型的界面正演

如圖1所示，在直角坐標系xyz中，假設地下密度界面A上方不同深度巖石與A下方巖石的密度差可用如下雙曲線密度模型表示[11]

(1)

式中： Δρ(z)為深度z處巖石與密度界面下方巖石的密度差； Δρ0為地表巖石與密度界面下方巖石密度差；β為衰減系數，實際應用中，β可由不同深度密度差計算。

圖1 三維雙曲線密度模型界面示意圖

圖1中，令z0為界面A的平均深度， Δh為界面A相對于z0的起伏深度，并令界面A因起伏引起的剩余密度體的體積元dv=dξdηdζ，其位置坐標為Q(ξ，η，ζ)。利用重力異常體積元積分法，可獲得密度界面A起伏引起的剩余密度體在地面上任一點P(x，y，0)產生的重力異常

Δg(x,y,0)

(2)

式中：G代表萬有引力常數；V代表體積分區域。對上式做快速傅里葉變換(FFT)，可得

(3)

式中u、v分別代表x、y方向的圓頻率。將式(1)代入式(3)，整理可得

(4)

(5)

令b=-β-z0，代入式(5)，可得

(6)

對式(6)中的ζ積分，可得頻率域密度界面的正演公式

F(Δg)=2πGΔρ0β2e-kz0{F(M1)+kF(M2)+

(7)

2.2 密度界面反演

將式(7)兩邊同時除以2πGΔρ0β2e-kz0，然后取出等號右邊第一項，得到界面起伏的反演公式

(8)

(9)

界面的迭代反演過程如下：

(1)確定雙曲線密度模型Δρ(z)及平均深度z0，計算重力異常的頻譜F(Δg)；

(2)利用式(9)和FFT逆變換計算界面起伏的一階近似值Δh(0)；

(3)將Δh(0)代入式(7)，經FFT逆變換，獲得近似重力異常Δg(0)及校正譜δg(0)=Δg-Δg(0)；

(5)校正界面起伏Δh(1)=Δh(0)+δh(0)；

(7)最終得到的界面深度為h=z0+Δh(p)。

計算中，可以采用正則化因子壓制向下延拓對高頻成分的放大作用[17]。本文采用的正則化因子為

(10)

式中：k0為需要壓制的最小頻率；α∈[2,7]。k0和α需要依據實際異常特征進行選擇。λx=mΔx是計算區沿x方向的長度， Δx是沿x方向的點距。

3 模型試驗

圖2(左)為三維界面模型，由132×124個計算點構成，x和y方向的點距均為1km。假設模型密度差僅隨深度變化，在觀測面(z=0)的巖石與界面位置處下部巖石密度差為-0.4g/cm3，地下4km(z=4km)處的密度差為-0.2g/cm3，由式(1)可計算出雙曲線密度模型的參數β=9.6569，界面模型因深度起伏引起的剩余密度在地面上產生的重力異常如圖2(右)所示。

為了檢驗本文方法的有效性和計算精度，利用本文方法和PO法對雙曲線密度模型產生的重力異常(圖2右)進行反演。在反演中，為了保證對比的公平性，兩種方法均采用了相同的正則化因子(α=2、k0=0.8)和平均深度(z0=4km)。PO法經6次迭代(圖3左)、雙曲線密度模型界面反演方法經4次迭代(圖3右)后均可取得精度較高的計算結果。對比分析二者與真實深度(圖2左)的擬合誤差(圖4)，可見PO法的反演結果在凸起和凹陷位置均存在較大的誤差(圖4左)，而雙曲線密度模型的反演結果(圖4右)與真實值更加接近，說明針對變密度地質體界面反演問題，本文提出的頻率域雙曲線密度模型界面反演方法有更高的計算精度。

圖2 模型界面深度(左)及理論重力異常(右)

圖3 PO法(左)與本文方法(右)的界面深度反演結果對比

圖4 PO法(左)和本文方法(右)反演深度誤差

4 實例應用

黑龍江省虎林盆地位于那丹哈達嶺燕山褶皺帶與吉黑華力西褶皺帶的交匯處，是在較為復雜的構造環境中形成的中、新生代凹陷帶，內部次級構造發育，也是大慶油田外圍東部油氣勘探的重點地區之一[18,19]。盆地基底的巖石組成、地層巖石密度及鉆井資料表明，研究區基底為遠古界黑龍江群、麻山群的深變質巖系及古生界淺變質的褶皺基底[20]，它們之間沒用明顯的橫向密度差異，地表巖石與基底巖石密度差為-0.3 g/cm3，盆地基底之上巖石與基底的密度差隨深度的變化趨勢可以用雙曲線密度模型模擬，利用數值模擬可以獲得模型參數為Δρ0=-0.3g/cm3，β=6.5574。

通過對盆地異常特征及巖石密度和鉆井資料分析，對比多個高度的向上延拓結果，認為利用虎林盆地布格重力異常(圖5)減去向上延拓5km的區域場獲得的剩余重力異常(圖6)反演盆地沉積蓋層的分布特征是合理的。

為了進一步驗證剩余場分離的合理性，對布格重力異常(圖5)進行了垂向一階導數處理(圖7)，其結果可以反映蓋層分布的基本特征。對比圖6與圖7，可見兩者等值線顯示的圈閉形態、位置基本相似，這也進一步佐證了剩余場分離的合理性。在此基礎上，經過大量試驗，選取平均深度為1.1km和正則化因子α=2.5、k0=1，作為PO法和本文方法的計算參數，反演了盆地的基底深度。相比于PO法反演結果(圖8)，本文方法的計算結果(圖9)在等值線形態上與垂向一階導數(圖7)的相似度更高(如吉祥以西)，說明本文方法對盆地基底形態的反演效果更好。電法勘探剖面可以較好地反映沉積盆地的基底特征[21,22]，為重力反演結果提供比對依據，故收集了盆地內一條大地電磁測深剖面(圖10，編號為DB1線，位置見圖8)，該斷面淺層顯示的低電阻率是沉積蓋層的反映。將DB1線上PO法與本文方法的計算結果做逐段對比，三處位置可以明顯體現出本文方法的優勢：電磁測深剖面和本文方法結果中A點均對應凸起頂部，而PO法結果中A點對應凸起和凹陷的轉折處；電磁測深顯示B點下方的基底深度為2.05km，本文方法為2.1km，PO法為1.8km； 電磁測深結果顯示C點對應凹陷的底部，本文方法中C點為局部最大深度值點，而PO法中C點處于凹陷的翼部(該凹陷底部位于C點以南3km處)。此外，本文方法反演的盆地基底深度(圖8)清晰地反映了盆地基底的凸起和凹陷分布特征，大部分可與文獻[19，20]的結論相互印證。盆地南部基底較深的凹陷為穆棱河凹陷和東林子凹陷，最深處分別為1.75km和2.05km； 北部為虎林河凹陷，最深處達到2.86km。其中，虎林河凹陷南側等深線密集，北側相對平緩，反映了南側斷裂的陡直特征。盆地南部和北部基底(圖8)表現出的走向、深度等細節信息說明了虎林盆地南部和北部基底結構存在差異，這與文獻[23]的結論也一致。上述對比分析結果可以證明本文方法的有效性和優勢。

圖5 虎林盆地布格重力異常圖

圖7 虎林盆地重力垂向一階導數圖

圖8 虎林盆地PO法計算的盆地基底深度圖

圖9 本文方法計算的盆地基底深度圖

5 結論

針對盆地內巖石密度隨深度變化的基本特征以及頻率域界面反演的優勢，本文給出了基于雙曲線密度模型的界面正演公式，以及界面深度的反演公式和迭代計算過程。模型試驗表明，相比于PO方法，本文方法的反演結果與真實深度更加接近。在虎林盆地實際資料處理中，本文方法的計算結果清楚地刻畫了基底的深度特征，基本上可與DB1線大地電磁測深結果及前人的研究結果相互印證，證明了本文反演方法的合理性和實用性。

圖10 虎林盆地DB1線MT電阻率斷面解釋結果(位置見圖8、圖9中點線)