商宇航 邰振華 秦 濤

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

1 引言

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

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

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

2 基本原理

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

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

(1)

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

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

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

Δg(x,y,0)

(2)

式中：G代表萬有引力常數(shù)；V代表體積分區(qū)域。對上式做快速傅里葉變換(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)兩邊同時(shí)除以2πGΔρ0β2e-kz0，然后取出等號右邊第一項(xiàng)，得到界面起伏的反演公式

(8)

(9)

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

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

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

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

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

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

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

(10)

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

3 模型試驗(yàn)

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

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

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

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

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

4 實(shí)例應(yīng)用

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

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

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

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

圖7 虎林盆地重力垂向一階導(dǎo)數(shù)圖

圖8 虎林盆地PO法計(jì)算的盆地基底深度圖

圖9 本文方法計(jì)算的盆地基底深度圖

5 結(jié)論

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

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