劉卓, 曾昭發(fā), 李靜，何榮欽，霍祉君，林景宜

1.吉林大學(xué) 地球探測科學(xué)與技術(shù)學(xué)院，長春 130026;2.國土資源部 應(yīng)用地球物理重點實驗室，長春 130026

0 引言

居里面是地殼中巖石鐵磁性礦物在溫度升高到居里點變?yōu)轫槾判晕镔|(zhì)時的深度界面，為磁法勘探的磁性層底界面，表征地下溫度場的分布，對地?zé)崽镌u價開發(fā)、油氣資源的預(yù)測、地震火山災(zāi)害防治及原生熱液礦產(chǎn)勘查等具有重要指導(dǎo)意義[1]。居里面深度計算方法包括譜分析方法和界面反演方法等。譜分析方法[2--3]是通過計算磁異常數(shù)據(jù)對數(shù)功率譜來確定磁源深度，但對數(shù)功率譜的復(fù)雜性和混亂性難以有效地劃分磁源平均深度，具有較大的誤差。Parker在位場界面正演計算中引入快速傅里葉變換(FFT)[4]。Oldenburg根據(jù)Parker公式，提出頻率域深度界面迭代反演方法[5]，改變了過去采用棱柱或長方體模型反演界面深度不連續(xù)性，同時計算速度快[6]。Parker--Oldenburg反演算法存在3個問題：①該方法屬于單磁性界面反演方法，磁性體通常是由雙界面或多界面組成的,利用單界面計算與實際情況不符。王萬銀[7]和相鵬[8]采用雙磁性界面模型來提高反演結(jié)果精度。②在反演過程中假設(shè)磁化率參數(shù)為常數(shù)，而實際地質(zhì)介質(zhì)磁化率參數(shù)在橫向和縱向上均具有較復(fù)雜的變化。馮娟等[9]采用了垂向變物性參數(shù)模型來計算基底界面，提高了反演精度。③在Parker--Oldenburg公式中包含向下延拓因子中的不收斂項，影響計算的穩(wěn)定性。盡管可以通過添加濾波器來確保迭代的收斂，但將會使信號的高頻信息損失，降低反演結(jié)果的分辨率和精度。肖鵬飛[10]、張沖等[11]采用多項式迭代方法來改進傳統(tǒng)Parker--Oldenburg算法，提高計算穩(wěn)定性。

針對上述問題，前人研究中分別針對不同問題提出了對應(yīng)的解決方案，沒有考慮綜合影響。筆者綜合了上述3個問題提出了一種改進Parker--Oldenburg界面的反演方法。通過引入垂向變磁化率因子和雙界面模型來提高模型的準(zhǔn)確性，同時采用多項式迭代算法來提高算法穩(wěn)定性，提高反演精度。通過理論模型測試和實測數(shù)據(jù)應(yīng)用驗證了該算法的可靠性以及結(jié)果的準(zhǔn)確性，計算的居里面深度更加具有理論和實際應(yīng)用價值。

1 基于雙磁性界面和垂向磁性變化的Parker--Oldenburg反演方法

1.1 基于雙磁性界面和垂向磁性變化的磁異常計算

垂向磁化率和雙磁性界面模型條件下頻率域正演磁異常公式為：(關(guān)于公式(1)推導(dǎo)過程詳見附錄)

(1)

對公式(1)進行重組，得到如下公式：

(2)

頻率域位場的向上延拓公式為：

F[U0]=e-ωh0F[U]

(3)

可以看出公式(2)和公式(3)在形式上具有相似性。

1.2 磁異常--深度變換的迭代公式

為迭代計算得出界面反演公式，首先假定地下某一深度處的磁異常初始值設(shè)為Δz(x,y,z)(1)=Δz(x,y,0)。利用公式(3)得出地表的異常值計算初值如下：

Δz(x,y,0)(1)=F-1[e-ωh0F[Δz(x,y,z)(1)]]

(4)

結(jié)合位場迭代計算方法[12]，可以利用Δz(x,y,z)(1)、Δz(x,y,0)、Δz(x,y,0)(1)獲得Δz(x,y,z)(2)，即：

Δz(x,y,z)(2)=Δz(x,y,z)(1)+s(Δz(x,y,0)-Δz(x,y,0)(1))

(5)

以此類推，經(jīng)過n次迭代計算后的結(jié)果：

Δz(x,y,z)(n+1)=Δz(x,y,z)(n)+s(Δz(x,y,0)-Δz(x,y,0)(n))

(6)

當(dāng)|Δz(x,y,0)-Δz(x,y,0)(n)|≤ε時，認為Δz(x,y,z)(n+1)≈Δz(x,y,z)(n)，對公式(6)推導(dǎo)可得：

(7)

提取深度結(jié)果得到：

(8)

公式(8)即為雙界面模型約束下的下界面反演公式。

2 數(shù)值模型模擬計算

圖1和圖2為設(shè)計的模型。從圖1可見模型由位于地下30 km和10 km處的兩個磁性界面組成，每個磁性界面均由兩個大小不同、方向相反的半球組成，30 km深度處底界面的較大凹的半球直徑為40 km，方向向下，最大埋深為40 km；較小隆起半球的直徑為30 km，方向向上，最小埋深22.5 km；10 km深度處頂界面的較大凹半球直徑也為40 km，最大埋深為20 km；較小隆起半球的直徑為30 km，最小埋深2.5 km。目標(biāo)體磁化率設(shè)為M0=1 A/m。為了對比分析，本文另設(shè)模型體的變磁化率設(shè)為M=M0e-0.178z。分別采用本文提出算法和傳統(tǒng)Parker--Oldenburg方法對圖1所示模型開展界面深度反演。提取模型俯視圖的界面對角線AB(圖2)剖面反演結(jié)果進行分析。

圖1 雙磁性界面模型體分布圖Fig.1 Dual magnetic interface model

圖2 模型體俯視圖對角線AB位置圖Fig.2 Diagonal AB top view of interface model

圖3為傳統(tǒng)Parker--Oldenburg方法界面深度反演結(jié)果，圖4為本文方法的磁異常界面深度反演結(jié)果。圖中黑實線為磁異常模型值，點劃線和黑虛線分別表示變磁化率和常磁化率情況下計算深度值。表1和表2為傳統(tǒng)算法與本文提出算法的計算參數(shù)對比。

2.1 變磁化率對界面深度反演的影響

對比圖3和圖4中點劃線和虛線的反演結(jié)果，可以看出應(yīng)用同一算法變磁化率模型下的反演結(jié)果要優(yōu)于常磁化率條件下的反演結(jié)果。

對比表1的數(shù)據(jù)可以看出，在應(yīng)用傳統(tǒng)Parker--Oldenburg方法條件下，變磁化率模型反演結(jié)果迭代20次剩余RMS誤差為0.019 2；常磁化率模型反演結(jié)果迭代30次剩余RMS誤差為0.019 3，明顯提高收斂速度。

圖3 傳統(tǒng)Parker--Oldenburg算法反演模型下界面深度結(jié)果Fig.3 Bottom interface depth inverted results obtained by conventional Parker--Oldenburg algorithm

圖4 新算法反演模型下界面深度結(jié)果Fig.4 Bottom interface depth inverted results obtained by new algorithm

數(shù)值模擬最大埋深/km最大深度位置/km最大埋深/km最小深度位置/km(87,87)點差值/km(36,36)點差值/km收斂準(zhǔn)則/km剩余RMS/km迭代次數(shù)P--O算法反演40.00(44,44)22.50(90,90)-----常磁化率模型39.47(44,44)25.82(93,93)3.481.300.020.01930030變磁化率模型40.09(44,44)24.83(91,91)2.481.500.020.01920020新算法反演40.00(44,44)22.50(90,90)-----常磁化率模型40.08(43,43)25.46(90,90)2.901.460.020.0052003變磁化率模型40.05(44,44)22.48(90,90)0.182.120.020.0000655

在應(yīng)用本文提出的算法變磁化率模型反演結(jié)果迭代5次剩余RMS誤差為6.5×10-5，常磁化率模型反演結(jié)果迭代3次剩余RMS誤差為0.005 2；變磁化率模型的RMS值要遠小于常磁化率模型的RMS值，表明變磁化率因子可以顯著地提高新法的反演精度和收斂的速度。

上述對比分析驗證了將變磁化率因素引入到界面反演公式中的確可以提高反演結(jié)果的精度和可靠性。

2.2 迭代算法對界面反演結(jié)果的影響

對比表2的數(shù)據(jù)可以看出，在應(yīng)用變磁化率模型條件下，新算法的反演結(jié)果迭代5次剩余RMS誤差為6.5×10-5；Parker--Oldenburg方法的反演結(jié)果迭代20次剩余RMS誤差為0.019 2；在應(yīng)用常磁化率模型條件下，新算法的反演結(jié)果迭代3次剩余RMS誤差為0.005 2；Parker--Oldenburg方法的反演結(jié)果迭代30次剩余RMS誤差為0.019 3。

上述對比分析驗證了新算法相較于Parker--Oldenburg方法具有更高的精度和效率。

表2 常規(guī)Parker--Oldenburg方法與本文所述方法反演參數(shù)對比

3 松遼盆地航磁數(shù)據(jù)居里面計算

松遼盆地是中國重要的含油氣盆地，也是中國的一個重要熱盆，地?zé)豳Y源豐富。根據(jù)磁異常計算區(qū)域的居里面深度，為分析區(qū)域油氣和地?zé)豳Y源，特別是深部干熱巖資源具有重要的意義。為了對松遼盆地地?zé)豳Y源進行有效的評估，將本文所述方法應(yīng)用于松遼盆地居里面深度計算。

3.1 異常分離

圖5為研究區(qū)化極磁異常。航磁異常的幅值范圍在-200～300 nT之間，異常值較弱。磁異常方向呈現(xiàn)NE向和SN向展布規(guī)律，局部磁異常呈EW向磁異常帶被NE向異常干擾、錯斷。磁異常方向變化反映了板塊構(gòu)造作用的應(yīng)力方向及方式的改變，反映了研究區(qū)后期主要受太平洋板塊俯沖作用的影響[13]?；瘶O磁異常特征分布的復(fù)雜性反映了松遼盆地漫長的地質(zhì)歷史中經(jīng)歷了多期次復(fù)雜構(gòu)造運動[14]，磁異常形態(tài)差異表明不同剛性程度各微板塊在外部應(yīng)力作用下產(chǎn)生了多樣的地質(zhì)構(gòu)造形態(tài)，這種構(gòu)造形態(tài)具有繼承性和疊加性[15]。

圖6為前人研究資料獲得的松遼盆地基底深度等值線圖[16]。為去除松遼盆地上覆地層的磁性影響，其磁性主要來源于泉頭組和登婁庫組的地層磁性，故選取泉頭組和登婁庫組地層的平均磁化率75×10-5SI[17]，并結(jié)合公式(1)計算得到松遼盆地基底及上覆地層的正演磁異常(圖7)，從總磁異常去除松遼盆地基底及上覆地層的正演磁異?？傻玫饺鐖D8所示的剩余磁異常，剩余異?？梢钥醋魇且跃永锩婧突诪殡p界面的模型體所產(chǎn)生的磁異常，為了引入隨深度變化的磁化率因子來約束居里面的反演結(jié)果，筆者參考前人研究資料中有關(guān)松遼盆地的井溫--深度數(shù)據(jù)[13]，通過數(shù)據(jù)擬合得到了松遼盆地地溫隨深度的變化。地溫隨深度變化關(guān)系為T=19.334e0.451 6z[13]，其中T為井溫，z為深度。大多數(shù)磁性物質(zhì)的磁性隨溫度呈反比例變化，故以泉頭組和登婁庫組地層的平均磁化率75×10-5SI為起算點，得到基底與居里面的之間的磁化率隨深度的變化關(guān)系為κ=3.88×10-5e-0.451 6z。

3.2 松遼盆地居里面計算結(jié)果及解釋

采用圖8的剩余磁異常利用變磁化率模型計算居里面深度如圖9所示。反演得到的居里面深度范圍為18～25 km，松遼盆地居里面反演結(jié)果顯示盆地中北部居里面深度最淺，向外依次變深。

圖10是松遼盆地的地溫梯度分布規(guī)律等值線圖，松遼盆地中北部地溫梯度的變化規(guī)律與居里面正好相反，二者的對應(yīng)關(guān)系滿足理論上的居里面深度與地溫梯度的反比例關(guān)系，地溫梯度高值區(qū)域與居里面深度淺的區(qū)域有很好的一致性；此外利用居里點溫度和反演得到的居里面深度與計算得到的地溫梯度和實測的地溫梯度在具體數(shù)值精度上也有良好的對應(yīng)關(guān)系。由此推測居里面上隆是松遼盆地地溫梯度高的深部原因。

圖5 松遼盆地化極磁異常Fig.5 Polarization magnetic anomaly of Songliao Basin

圖6 松遼盆地基底深度Fig.6 Basement depth of Songliao Basin

圖7 松遼盆地基底及上覆地層正演磁異常Fig.7 Forward magnetic anomalies of basement and overlying strata in Songliao Basin

圖8 松遼盆地去除基底及上覆地層原始磁異常后的剩余磁異常Fig.8 Residual magnetic anomalies obtained by removing forward anomalies of basement and overlying strata from original magnetic anomalies

圖9 松遼盆地居里面深度反演結(jié)果Fig.9 Inverted Curie depth of Songliao Basin

圖10 松遼盆地地溫梯度圖[14]Fig.10 Geothermal gradient map of Songliao Basin

居里面反演結(jié)果與通過測井資料獲得的地溫梯度資料無論是在分布形態(tài)還是在數(shù)值精度上均有良好的對應(yīng)關(guān)系，通過對居里面反演結(jié)果與實測地溫梯度之間的對比研究，驗證了本文獲得的居里面反演結(jié)果具有較高的反演精度。

在去除了基底磁性的影響后，剩余磁異常在EW向上的構(gòu)造特征相比于原始化極磁異常EW向上的構(gòu)造特征更加明顯，表明松遼盆地深部EW向構(gòu)造被NS向構(gòu)造所干擾、錯斷[15,18--19]。

He[20]等根據(jù)深部地震Vp/Vs比值和接收函數(shù)成像結(jié)果，發(fā)現(xiàn)在松遼盆地中部地下410 km和660 km深度不連續(xù)帶存在凹陷，推斷深部存在地幔熱柱上涌現(xiàn)象。地幔柱的位置與本文計算居里面隆起(圖9)吻合，認為松遼盆地中北部區(qū)域居里面隆起與地幔熱柱上涌相關(guān)[21--22]。

松遼盆地中北部基底存在深大斷裂構(gòu)造，為深部熱源提供良好的導(dǎo)熱通道；松遼盆地深部基底分布有大量花崗巖，存在背斜構(gòu)造，為熱量的傳導(dǎo)與儲存提供良好的地質(zhì)條件。居里面隆起的區(qū)域是地?zé)岚ǜ蔁釒r重要的有利區(qū)，推斷圖9中居里面深度較淺的松遼盆地中北部區(qū)域地?zé)豳Y源豐富，結(jié)合松遼盆地地溫梯度圖(圖10)，可以將地?zé)峥碧竭h景區(qū)縮小至地溫梯度最高值附近，獲得圖9所示的1、2、3指示的區(qū)域。結(jié)合前人的研究資料[23--24]，遠景區(qū)2區(qū)域存在大量花崗巖分布，同時處于背斜構(gòu)造，利于熱量的傳導(dǎo)和儲存；遠景區(qū)3的花崗巖分布較少，處于向斜構(gòu)造；遠景區(qū)1，沒有明顯的褶皺構(gòu)造，花崗巖分布最少。通過分析1、2、3區(qū)域的花崗巖分布特征和褶皺構(gòu)造，驗證了2區(qū)域的地?zé)衢_發(fā)遠景最優(yōu)，3區(qū)域其次，1區(qū)域最次。因此在進一步勘探中建議優(yōu)先在2區(qū)域中尋找。

4 結(jié)論

(1) 通過模型對比分析，驗證了將變磁化率引入到界面反演公式中可以提高反演結(jié)果精度。

(2) 新迭代算法的應(yīng)用在保證反演結(jié)果收斂的同時，避免了濾波器的使用，反演過程中保留原始數(shù)據(jù)所有頻率信息，提高了反演結(jié)果的準(zhǔn)確性和穩(wěn)定性。

(3) 通過理論模型對比分析，本文提出相比于常規(guī)Parker--Oldenburg方法的反演計算，具有計算效率更快，精度更高，泛化能力更強等優(yōu)勢。

(4) 利用本文提出的界面反演算法獲得了松遼盆地居里面深度。松遼盆地的居里面表現(xiàn)為中北部相對較淺，而南部地區(qū)和北部地區(qū)居里面深度較深。

(5) 通過計算結(jié)果的分析，并綜合區(qū)域地質(zhì)和鉆井資料，為松遼盆地后期深部構(gòu)造活動對前期深部構(gòu)造形態(tài)的繼承和改造方面的地質(zhì)認識提供基礎(chǔ)。

(6)勾畫了3個地?zé)豳Y源勘探遠景區(qū)，推測出松遼盆地中北部居里面較淺的地區(qū)具有良好的地?zé)衢_發(fā)遠景。

附錄

由頻率域泊松公式可得，

(A-1)

式中:F[]為傅里葉變換;Um為磁位;V為重力位;M為剩余磁化強度;ω為圓波數(shù)；G為萬有引力常數(shù);ρ為剩余密度。當(dāng)磁化率隨深度變化時，可以表示為：

M=M0eaζ

(A-2)

式中:M0是地表地質(zhì)介質(zhì)的剩余磁化強度;a是磁化率隨深度變化因子。

將公式(A-2)代入公式(A-1)中，獲得如下式(A-3)：

(A-3)

磁異常的波譜可以表示為如下形式：

F[Δz]=-μ0ωF[Um]

(A-4)

式中:Δz為磁異常;μ0為真空中磁導(dǎo)率。

重力異常的波譜可以表示為如下形式：

F[Δg]=ωF[V]

(A-5)

綜合公式(A-1) (A-2) (A-3) (A-4) (A-5)，可得

(A-6)

頻率域雙界面模型重力異常正演公式的推導(dǎo)過程，其中頻率域重力異常正演公式可表示成如下形式[7]：

(A-7)

式中:h0為上下界面的平均深度;Δh1為下界面和平均深度的差值;Δh2為上界面和平均深度的差值(已知)。結(jié)合公式(A-7)，公式(A-6)中eaζF[Δg]可化為如下形式：

(A-8)

為頻率域密度隨深度呈指數(shù)變化約束下的雙界面模型重力異常正演公式，并將公式中eζ(a-ω)在ζ=h0處進行Talyor展開[9]，即可得：

(A-9)

將式(A-9)代入式(A-6)可得

(A-10)

公式(A-10)即為垂向磁化率和雙磁性界面模型條件下頻率域正演磁異常公式。