基于結(jié)構(gòu)模型的磁性基底反演

于會臻，劉 展

（中國石油大學(xué) 地球資源與信息學(xué)院，山東青島 266555）

基于結(jié)構(gòu)模型的磁性基底反演

于會臻，劉 展

（中國石油大學(xué) 地球資源與信息學(xué)院，山東青島 266555）

常規(guī)磁性基底反演方法往往忽略基底磁性變化采用常磁化強度模型，這里分析了其中的局限性，提出了先構(gòu)建相關(guān)搜索重磁三維定量反演技術(shù)的變磁化強度模型，然后通過Parker界面反演算法進(jìn)行基底反演。其中，對反演磁化強度所需磁異常分離采用了小波變換方法，并分析了該方法的關(guān)鍵技術(shù)和措施。通過利用該方法對某工區(qū)進(jìn)行了實際資料處理后，得到較好效果。

磁性基底；變磁化強度；Parker；相關(guān)系數(shù)；多尺度分析

0 前言

由于古中生代及第三紀(jì)沉積巖層和第四紀(jì)疏松沉積為無磁性、或者弱磁性的，所以磁性基底作為磁性地層的上界面時，通常為前寒武紀(jì)的結(jié)晶基底，由巖漿巖或變質(zhì)巖組成。在一般情況下認(rèn)為，居里等溫面為磁性地層的底界面。研究磁性基底界面無論對于研究盆地的形成與演化，還是對于研究礦藏的分布（如油氣形成、運移、賦存等），都有極其重要的理論和實際價值［1］。目前，磁性界面反演方法主要有切線法、人機交互解釋法、Parker迭代反演算法、歐拉算法、模型反演法、頻譜分析法、樣條函數(shù)法［2］。在實際應(yīng)用中，主要采用的是切線法和人機交互解釋法。常規(guī)反演方法通常采用簡單的常磁化強度模型進(jìn)行處理，但由于磁性基底相關(guān)異常是由基底界面起伏和基底磁性變化共同引起的，因此，往往出現(xiàn)反演結(jié)果與實際情況相悖的現(xiàn)象。針對該問題，作者在本文提出了先構(gòu)建變磁化強度模型，再進(jìn)行Parker磁性界面反演算法［6］的新磁性基底反演方法。

1 方法原理

Parker算法是七十年代R.L.Parker提出的一種界面重磁正反演公式，由于它能計算物性橫向變化的連續(xù)界面，速度快，所以很快得到了廣泛的應(yīng)用。其基本原理為：

設(shè)地下界面S上部磁化強度J1為零、下部磁化強度J2不為零，平均深度為H，其在海拔水平面上的觀測磁異常為式（1）。

式中 m表示分步積分的項數(shù)；n是指數(shù)項展開成級數(shù)的項數(shù)；Ml成為磁化強度函數(shù)。

ζn－m－l為界面函數(shù)，F(xiàn)［Mlζn－m－l］為界面函數(shù)的頻譜，則界面的反演迭代公式見式（2）及式（3）。

式中 ω為徑向圓頻率，第（j－1）次得近似值，右端的為第i次的近似值，一般先給定△h初值，代入公式右端，求出△h，經(jīng)反變換得到△h1；然后作為△h的下一次迭代的初始值，如此這樣反復(fù)直到結(jié)果滿足給定精度［2］。

針對算法中高頻因子eωh影響迭代收斂性的問題，前人做出了有效的解決方案［6］，其算法雖然對磁化強度Ml對進(jìn)行了假設(shè)，但用一個通用公式很難描述實際情況的磁化強度的分布。事實上，不僅僅是Parker界面反演算法，目前磁性基底反演中都或多或少地存在缺乏磁化強度約束問題，通常采用均一值視磁化強度模型（即設(shè)磁性基底的地層磁性參數(shù)為一常數(shù)），很少對基底磁性變化進(jìn)行討論［2］。

按均一值視磁化強度模型做界面正演（見圖1）可知，忽略基底磁性變化，會使得對觀測磁異常的處理全部假設(shè)為完全由磁性界面的構(gòu)造引起。圖1中的磁性基底淺會產(chǎn)生磁異常高值，磁性基底深會出現(xiàn)磁異常低值。因此，不合理的磁化強度模型假設(shè)得到的磁性基底反演結(jié)果，可能會在某些區(qū)域呈構(gòu)造反轉(zhuǎn)（見圖1）。為解決該問題，作者將利用相關(guān)搜索重磁三維定量反演方法，構(gòu)建一個接近實際情況的磁化強度模型，以便在Parker基底界面反演方法中增加有效的基底磁性約束。

圖1 構(gòu)造反轉(zhuǎn)示意圖Fig.1 Tectonic reversal

1.1 視磁化強度模型的構(gòu)建

1.1.1 視磁化強度反演方法原理

與界面反演一樣，磁化強度反演也存在嚴(yán)重的多解性問題，從公式（1）中可以看出，它們兩者其實互為約束，即在求解過程中存在著互相制約的矛盾關(guān)系。同時，由兩者產(chǎn)生的磁異常是無法由觀測異常中準(zhǔn)確分離的。因此，需要構(gòu)建的磁化強度模型不可能是一個精確的模型，其目標(biāo)應(yīng)是找到由已知地質(zhì)資料及觀測數(shù)據(jù)所反映磁性基底磁性宏觀分布，采用的方法為概率成像與位場約束擴展結(jié)合的方法。

在作者提出以下模型構(gòu)建方法之前，首先假設(shè)已得到了磁化強度所引起的基底相關(guān)磁性磁異常磁異常的分離方法，將在1.2節(jié)中進(jìn)行詳細(xì)討論。

由本節(jié)開頭所得結(jié)論，磁化強度模型構(gòu)建的目標(biāo)是一種基底磁性的宏觀反映，因此對物性單元的剖分也應(yīng)體現(xiàn)這一點。將磁性基底界面與居里面之間地層剖分成若干個規(guī)則直立六面體形體單元（見下頁圖2），六面體高度的不同反映了地層的起伏，反演所得六面體磁化強度值主要反映了基底磁性橫向的宏觀變化，該處磁化強度定量反演方法采用了文獻(xiàn)［3］提出的，以單相關(guān)系數(shù)為基礎(chǔ)構(gòu)制目標(biāo)函數(shù)的多維黃金分割法。其關(guān)鍵技術(shù)是引入了多元線性回歸分析中的相關(guān)系數(shù)R，具體表達(dá)式見參考文獻(xiàn)［3］。R的物理意義是模型理論場與觀測場之間的誤差方差與觀測場本身方差的逼近程度其目的是求取模型與觀測數(shù)據(jù)的擬合組成的目標(biāo)函數(shù)的極小來獲得物性值的解。R的取值為0～1越接近“1”說明觀測場與核函數(shù)的相關(guān)性越好；反之則越差。使得R最接近“1”的物性模型，就是反演最優(yōu)解。算法在解的逼近問題上采用了n維黃金分割最優(yōu)化方法，反演速度較快。但該算法如大多數(shù)磁化強度反演算法一樣，對已知鉆井信息的利用率較低且反演采取由淺到深的順序易產(chǎn)生上層剖分網(wǎng)格占據(jù)大部份異常值的“上漂”現(xiàn)象。

為得到更合理的反演結(jié)果，作者提出了雙重約束機制。該機制的提出鑒于以下假設(shè)：①由于剖分網(wǎng)格與觀測場之間有一定的映射關(guān)系，找到這種關(guān)系即可確定任一剖分網(wǎng)格對于整個場的貢獻(xiàn)大小按照貢獻(xiàn)值大小進(jìn)行反演，可使存在場源體可能性大的區(qū)域先分配場值從而避免“上漂”；②由于地質(zhì)體賦存空間內(nèi)磁化強度值必然存在著一定的連續(xù)性且能體現(xiàn)在所引起的磁場相關(guān)性特征中，以往簡單的鉆井點約束可利用該特征對附近剖分網(wǎng)格進(jìn)行磁化強度約束空間的擴展。根據(jù)以上假設(shè)，約束分為磁化強度概率成像約束和基于鉆井的位場特征約束。

圖2 磁性基底剖分圖Fig.2 Magnetic basement dissection figure

（1）視磁化強度相關(guān)概率成像約束。利用概率成像方法［4］，通過位場數(shù)據(jù)與場源體核函數(shù)相關(guān)獲得場源體成像范圍進(jìn)行約束信息提取。當(dāng)觀測面為水平面時，位于地下q（xq，yq，zq）點磁體元在XYZ三個方向分量的概率函數(shù)為：

其意義是，通過計算場源體與實測磁異常的核函數(shù)與實測磁異常，在一定窗口范圍內(nèi)的歸一化互相關(guān)得到觀測場，是由任一剖分單元產(chǎn)生的概率，并利用該值搜索場源體磁化強度分布。由于ηq為一互相關(guān)結(jié)果，其絕對值不大于“1”。該值越大，表明q點越可能存在異常源；反之，表明存在場源體的可能性越小。利用磁異常垂直分量，可計算出，取其中絕對值最大者的數(shù)值和符號作為q點場源發(fā)生的概率。ηq值為正說明存在的異常源與掃描函數(shù)中的場源極性相同，為負(fù)值表明極性相反。

當(dāng)?shù)玫礁怕手岛?，選取0～1間的閾值對其進(jìn)行二值化。經(jīng)多次試驗證明，閾值選0.5～0.8之間比較合適。將概率值小于閾值的賦“0”，大于閾值的賦“1”，之后選取一定的搜素半徑對二值化結(jié)果進(jìn)行八方位領(lǐng)域搜索，將在搜索半徑內(nèi)相鄰都為“1”的剖分單元組成一個閉合空間，便可找出場源體磁化強度的成像范圍。利用該成像范圍既可以概率值的大小確定上述多維黃金搜索算法中的搜索順序，又可為后續(xù)磁場特征約束擴展方法做準(zhǔn)備。

（2）基于鉆井的位場特征約束。一般的反演方法僅局限于沿井筒進(jìn)行線約束，即只對井筒穿過的單元進(jìn)行約束。針對該問題，作者借助場源體分布與磁異常有較強的相關(guān)性，將鉆筒約束擴展為體約束（見圖3）。

圖3 鉆井約束擴展示意圖Fig.3 Quantitative constraint expansion

任一場源點（ξ，γ，ζ）在計算點（x，y，z）處的磁異常寫成如下形式（僅考慮垂直磁化的情況）：

其中 △T為或磁力異常；C是不依賴于xyz的常數(shù)；ρ為測點和場源點的距離；ρ＝［（ξ－x）2＋（γy）2＋（ζ－z）2］1/2。

當(dāng)場源與觀測面足夠遠(yuǎn)或場源足夠小時，近似為點源，N＝2。對式（8）分別對x、y、z求導(dǎo)，得：

表1（見后面）是通過模型試驗研究得到的磁異?！鱐的半極值、水平一階導(dǎo)數(shù)△Tx的正負(fù)極值、垂向?qū)?shù)△Tz和△Tzz的零值水平區(qū)間距L與場源體橫向分布L0的之間的關(guān)系。

由下頁表1可知，L/L0的比值越接近“1”，說明對異常體水平邊界刻畫越準(zhǔn)確。因此在同一埋深時，利用磁異常的垂向一次導(dǎo)數(shù)零值范圍圈定場源體橫向賦存空間最準(zhǔn)確。按該結(jié)論便可對鉆井附近的剖分網(wǎng)格磁化強度值約束區(qū)間進(jìn)一步擴展，視位場特征邊界范圍內(nèi)的剖分網(wǎng)格為磁化強度值接近的場源體。

1.1.2 磁性基底磁性磁異常分離與基底反演

針對磁性基底反演可將實測總異常劃分為：

（1）基底上部地層引起的磁異常。該異常趨勢變化主要是由淺部地層的強磁性體（如侵入性火成巖）引起的，由磁異常正演公式知，其離觀測面距離相對較近，因此通常具有頻率高、變化劇烈的特點。而且其中相當(dāng)一部份強磁性體分布與斷裂有較好的相關(guān)性。

（2）基底界面與居里等溫面之間地層引起的磁性基底相關(guān)磁異常。該異常中包含了磁性基底起伏與基底磁性變化共同的影響。不同磁性的變質(zhì)巖作用不同：副變質(zhì)巖磁性較弱，正變質(zhì)巖磁性較強。大部份區(qū)域地層磁性相對較弱，其引起的磁異常可看作是一種構(gòu)造異常，反映的是基底起伏的貢獻(xiàn)，為區(qū)域異常。而各種強磁性變質(zhì)巖的橫向分布變化所引起的磁異常，更多地反映了基底磁性分布情況，為局部異常場。因此，磁性基底相關(guān)異常可表達(dá)為弱磁性區(qū)域異常f0（M0）、中等磁性異常f1（M1）和強磁性異常f2（M2）三部份。其中M0、M1、M2分別代表弱磁性、中等磁性、強磁性區(qū)域。同時從實際勘探情況來看，基底起伏通常比較平緩，f0（M0）頻率相對較低，而f1（M1）和f2（M2）頻率相對較高。要對基底磁性進(jìn)行宏觀描述，其實就是尋找與中、強磁性異常特征。

上述的討論可按如下步驟對觀測異常進(jìn)行分離。首先消除火成巖等淺部地層磁性體的影響，得到磁性基底相關(guān)磁異常，該異常是基底界面反演所對應(yīng)的磁異常。圈定火成巖的方法為利用重力數(shù)據(jù)垂向二次導(dǎo)數(shù)處理、地震資料劃分工區(qū)內(nèi)主要斷裂構(gòu)造，并借助航磁垂向一次導(dǎo)數(shù)精確圈定出火成巖位置。之后，利用多次滑動趨勢分析選取合適的半徑逐步從實測異常中消除火成巖的干擾。

其次，從磁性基底相關(guān)磁異常中分離f0（M0）與f1（M1）、f2（M2），該異常分離得到的是磁化強度反演對應(yīng)的磁異常。作者提出了小波變換的互相關(guān)系數(shù)方法，小波分離的關(guān)鍵在于小波函數(shù)的選取與分解階次的確定。磁性基底相關(guān)磁異常的n階尺度小波逼近反映的是具有低頻特征的區(qū)域磁異常隨著尺度的增大，小波逼近對應(yīng)的頻率信息也逐漸降低。與之對應(yīng)的磁性基底相關(guān)磁異常與小波逼近的差值為剩余磁異常，反映的是淺層高頻信息。由本節(jié)開始的定性分析討論可預(yù)測，當(dāng)小波分解尺度到某一階數(shù)n后，小波逼近對應(yīng)的低頻信息將對應(yīng)僅由弱磁性地質(zhì)體及界面引起，而且，由于小波具有時～頻定位特點，小波逼近系數(shù)（an＋i（i＝1，2…，m））分布不會再隨尺度繼續(xù)增大出現(xiàn)明顯變化。因此，可采用相鄰尺度下的磁性基底相關(guān)磁異常小波逼近系數(shù)的互相關(guān)系數(shù)來判斷分解尺度。隨著分解階數(shù)的增加，an＋i（i＝1，2，…，m）相鄰互相關(guān)系數(shù)越來越接近“1”；尺度an之前的an－1次小波重構(gòu)函數(shù)所反映的是基底中、強磁性區(qū)域所引起的磁異常，因此在兩個尺度的小波系數(shù)分布之間，必然存在較小的相關(guān)性且呈震蕩起伏。在互相關(guān)系數(shù)呈漸進(jìn)增大時，前an－1尺度下的小波重構(gòu)函數(shù)就可以作為基底磁性相關(guān)異常。需要注意的是，在計算互相關(guān)系數(shù)前，需要對小波逼近系數(shù)進(jìn)行歸一化處理。在互相關(guān)系數(shù)突變點不明顯時，可以通過功率譜分析的方法，確定其各階所反映場源的平均深度。利用Bowin的公式，計算相應(yīng)深度的大地水準(zhǔn)面異常與小波分解的細(xì)節(jié)結(jié)果對比，從地球物理意義方面，可增加最高階次分解結(jié)果的合理性［14］。

2 模型試驗

測網(wǎng)為20 000m×20 000m，點距、線距都為500m，參數(shù)如下頁表2所示。模型1～3為三處淺層火成巖磁化強度值，模型4～10為磁性基底部份強磁化強度值。磁性基底平均磁化強度為50× 10－1A/m，試驗中設(shè)實際磁性基底為深度5km的平界面，居里等溫面為深度20km的平界面。對視測磁異常進(jìn)行3km趨勢化分析，剝離掉火成巖引起的磁異常，則磁性基底相關(guān)區(qū)域磁異常如下頁圖4所示。作者對該結(jié)果選用了“Bior3.5”小波進(jìn)行16階的離散小波變換，得到小波系數(shù)的互相關(guān)系數(shù)如下頁圖5所示，可知尺度8的重構(gòu)函數(shù)為基底磁性相關(guān)磁異常。概率成像閾值選為0.5，進(jìn)行重磁三維定量反演得到視磁化強度模型（見下頁圖4（b））。最后利用圖4（b）的磁異常和基底及居里等溫面的初始模型進(jìn)行Parker磁性界面反演算法，得到磁性基底界面如圖6（b）（見后面）。

表1 異常特征值水平范圍L與異常體寬度2a的關(guān)系Tab.1 Abnormal eigenvalue level Land abnormal body width 2arelationship

表2 模型參數(shù)Tab.2 Model parameters

圖4 模型磁異常處理Fig.4 Models of magnetic anomaly processing

結(jié)果分析可知：①下頁圖6（a）為利用均一化磁化強度模型的磁性基底反演結(jié)果，反響基底反演結(jié)果起伏較大，深度分布在5 000m±230m左右，且走勢與觀測磁異常變化相近；②利用本文中的方法構(gòu)建的基底磁化強度模型，與初始模型比較吻合，所得基底界面深度（見下頁的圖6（b））分布在5 000m ±20m。采用本文的方法可較好地降低出現(xiàn)構(gòu)造反轉(zhuǎn)的可能性，增加了反演結(jié)果的合理性。

3 野外資料處理

表3 鉆井深度Tab.3 Drilling depth

通過對中國遼東灣某工區(qū)進(jìn)行磁性基底反演驗證了方法的適用性。該區(qū)磁測資料已進(jìn)行化極和曲化平處理。該區(qū)域七口鉆井集中于工區(qū)中南部，深度達(dá)2km左右（如表3所示），同時結(jié)合該工區(qū)構(gòu)造，經(jīng)分析知磁性基底界面為西低東高，西部深處可達(dá)為5km左右，居里面深度在20km左右，對其中間地層進(jìn)行網(wǎng)格剖分可知，X、Y剖分間隔為500m，Z方向不剖分。

作者首先利用重磁資料對工區(qū)內(nèi)斷裂分布進(jìn)行劃分，同時由化極磁異常的垂向一次導(dǎo)數(shù)進(jìn)一步確定火成巖位置，選用5km半徑趨勢化分析得到去除火成巖后的磁性基底相關(guān)磁異常（見下頁圖7）；接著，選用“Bior3.5”小波對磁性基底相關(guān)磁異常進(jìn)行16階的離散小波變換，則由小波變換尺度9的重構(gòu)函數(shù)得到基底磁性相關(guān)磁異常。由已知的地質(zhì)資料可粗略確定每個剖分單元的視磁化強度值的上限和下限，視磁化強度相關(guān)概率成像約束中閾值選0.7，經(jīng)反演得到視磁化強度分布；最后，借助得到的磁化強度模型，利用Parker算法進(jìn)行磁性基底深度界面，反演結(jié)果如上頁圖8所示。由此可見，工區(qū)內(nèi)磁性基底呈東西成帶、南北分塊的特征，并且南低北高、西低東高，對比已知井資料可知誤差在100m～200m之間，結(jié)合已知地質(zhì)認(rèn)識，反演結(jié)果趨勢及深度比較合理。

4 結(jié)束語

作者在本文對磁性基底反演采用了一種新的思路和嘗試，在研究Parker界面反演算法的基礎(chǔ)上，提出了基于單相關(guān)系數(shù)多維黃金分割搜索方法的磁場概率成像及位場特征值約束機制的磁化強度模型構(gòu)建方法。其中，反演過程中的磁異常分離利用了小波不同尺度變換之間的互相關(guān)關(guān)系，對模型及實際的工區(qū)的處理得到了較好的效果。

由于位場反演本身存在的嚴(yán)重多解性問題，還有許多處理步驟需要進(jìn)一步研究，構(gòu)建更加符合磁異常特征的小波基函數(shù)，將會對異常的提取具有更強的物理意義和更好的實際效果。磁化強度模型構(gòu)建所采用的概率成像掃描函數(shù)，可以采用導(dǎo)數(shù)方程，這樣更有利于對磁性體邊界的確定，鉆井約束擴展方程的進(jìn)一步修改也將會提高約束精度。

［1］ 劉光鼎，郝天珧，劉伊克.重磁研究對認(rèn)識盆地的意義［J］.地球物理學(xué)進(jìn)展，1996，11（2）：1.

［2］ 趙百民，郝天珧.反演磁性地質(zhì)界面的意義與方法［J］.地球物理學(xué)進(jìn)展，2006，21（2）：353.

［3］ 劉展，班麗，魏巍，等.濟陽坳陷花溝地區(qū)火成巖重磁成像解釋方法［J］.中國石油大學(xué)學(xué)報，2007，31（1）：30.

［4］ RAFFAELE ALAIA，DOMENICO PATELLA，PAOLO MAURIELLO.Imaging multipole gravity anomaly sources by 3Dprobalility tomography［J］Journal of Geophysics and Engineering，2009，3（6）299.

［5］ 班麗.相關(guān)約束重磁三維定量反演方法研究［D］.東營：中國石油大學(xué)，2009.

［6］ 相鵬，劉展.雙界面模式Parker算法磁性界面正反演方法［J］.中國石油大學(xué)學(xué)報：自然科學(xué)版，2009，33（1）：37.

［7］ 汪炳柱，徐世浙，劉保華，等.多次插值切割法分場的一個實例［J］.石油地球物理勘探，1997，32（3）：229.

［8］ 段木春，徐世浙，閻漢杰，等.劃分磁異常場的插值切割法在研究火成巖體分布中的應(yīng)用［J］.石油地球物理勘探，1998，33（1）：125.

［9］ MARK E.ANDER，STEPHEN P.Huestisl.Gravity ideal bodies［J］.Geophysics.1987（52）：1265.

［10］孫魯平.基于密度成像的綜合速度建模方法研究［D］.東營：中國石油大學(xué)，2007.

［11］PAOLO MAURIELLO，DOMENICO PATELLA Localization of maximum－depth gravity anomaly sources by a distribution of equivalent point masse［J］.Geophysics，2001，66（5）：1431.

［12］姚長利，郝天珧，管志寧，等.重磁反演約束條件及三維物性反演技術(shù)策略［J］.物探與化探，2002，26（4）253.

［13］段木春，范典高.利用插值切割法研究磁性基底局部起伏特征［J］.石油物探，1999，38（4）：89.

［14］刁博，王家林，程順友.重力異常小波多分辨分析分解階次的確定［J］.中國地質(zhì)大學(xué)學(xué)報，2007，32（4）564.

［15］于德武.用等效磁源法進(jìn)行磁異常轉(zhuǎn)換［J］.物探化探計算技術(shù)，2004，26（2）：133.

［16］吳文鸝，管志寧，高艷芳，等.重磁異常數(shù)據(jù)三維人機聯(lián)作模擬［J］.物探化探計算技術(shù)，2005，27（3）：227.

［17］梅巖輝.均質(zhì)多面體重磁異常正演計算表達(dá)式的一致性［J］.物探化探計算技術(shù)，2007，29（1）：33.

book=107,ebook=107

1001—1749（2012）03—0295—08

P 631.4

A

10.3969/j.issn.1001－1749.2012.03.10

于會臻（1981－），男，漢族，山東桓臺人博士，研究方向為地球探測與信息技術(shù)。

國家重大專項（2008ZX05020－006）

2011－07－31改回日期：2012－03－20