利用迭代反演方法解決變磁性磁界面反演問題

張向宇，吳健生

（1.廣州海洋地質(zhì)調(diào)查局 ，廣州 510760；2.同濟大學(xué) 海洋地質(zhì)國家重點實驗室，上海 200092）

南海是環(huán)太平洋西部的最大邊緣海，面積約350×104km2。深入開展南海北部中生代沉積盆地研究，對于開拓我國海洋油氣勘探新領(lǐng)域具有十分重要的意義。近年來，隨著勘探工作的開展，由于地震方法受射線覆蓋強度、反射能量有限及分辨尺度等因素的限制，使地震資料難以清晰地反映盆地凹陷下埋深較大的基礎(chǔ)層的構(gòu)造，從而影響對基底結(jié)構(gòu)的了解和構(gòu)造演化的推斷，因此使得重磁勘探的手段得到了進一步的重視。

南海地區(qū)經(jīng)歷了多次構(gòu)造運動，使得地層巖石磁性在橫向上呈現(xiàn)較大差異，目前常用的磁界面反演方法如歐拉算法、頻譜分析法、維納反褶積法以及特征點法和切線法等，都是針對常磁性磁界面的計算，若直接將其應(yīng)用于南海北部磁性界面的反演中必定會帶進很大誤差，這就給常用反演方法的應(yīng)用提出了挑戰(zhàn)。通過以往的研究，對目前較常用的迭代反演方法進行改進，獲得適用于變磁性磁界面反演的方法，為類似南海地區(qū)這種磁性變化較大區(qū)域的磁性基底反演提供了一種解決方法。其次磁性界面可能與基底界面有關(guān)聯(lián)，這樣通過反演磁性界面可以得到基底的情況，這也為進一步研究地質(zhì)—地球物理解釋工作提供了依據(jù)。

國內(nèi)、外有關(guān)學(xué)者對磁界面迭代反演方法有了較多研究，國外方面，Parker［1］在1972年提出了界面正演方法；Oldenburg D W［2］提出了起伏界面的迭代反演方法，此方法一經(jīng)提出，便得到了眾多地球物理學(xué)者的重視，并得到了廣泛運用；Guspi［3］在Oldenberg提出方法的基礎(chǔ)上，提出了一種非線性迭代反演法；Barbosa［4］對Guspi提出的方法進行了改進，解決了不宜編程實現(xiàn)的缺點。國內(nèi)方面，趙百民等［5］對用于磁界面反演的一般方法如Parker迭代反演方法、歐拉算法等進行了簡要介紹；相鵬等［6］進行了基于Parker算法的磁性雙界面正、反演研究，并應(yīng)用理論模型進行研究驗證，且運用到實測數(shù)據(jù)中；董煥成等［7］、王家林等［8］、管志寧等［9］也在相關(guān)書籍中分別對界面迭代反演、直接反演法、統(tǒng)計反演法、線性反演法等反演方法的常用方法原理進行了介紹和應(yīng)用?；谇叭藢Υ沤缑娣囱莘椒ǖ难芯浚瑸楸疚淖兇判源沤缑娴囱萏峁┝艘欢ǖ囊罁?jù)。

1 磁界面迭代反演方法原理

1.1 常磁性磁界面迭代反演方法

設(shè)界面以下磁化強度恒定，以上無磁性（M ＝0），取過平均深度的水平面為起算平面，如圖1所示。根據(jù)引力位公式及場位關(guān)系式，考慮界面起伏遠小于平均深度，得到垂直磁化垂直磁異常公式（在xoy平面上）為式（1）。

圖1 磁界面反演計算示意圖Fig.1 Map of inversion for magnetic interface

其中 h0為界面平均深度；Δh＝h－h(huán)0為界面相對平均深度的起伏；Mz為磁化強度M 在z方向上的分量；μ0為真空磁導(dǎo)率。

式（1）寫成近似積分形式為式（2）。

式（2）為2 m×2 n個垂直磁化有限延伸直角棱柱體的垂直磁異常之和，即垂直磁化磁界面垂直磁異常的近似公式。當(dāng)已知磁異常值和磁化強度值時，通過解線性方程組就可得到界面起伏值，實際計算時應(yīng)采用迭代法計算。其步驟為：①將實際的Za值作為迭代的值，并給出Δhij的初值Δ；②將Zadie（k，l）帶入式（2）解出迭代中的 Δhij的改正值δΔhij，將δΔhij與相加后再代入式（2）計算的理論值；③計算理論值與實際值的均方誤差，當(dāng)滿足迭代終止條件時即完成了反演計算［7］。

1.2 變磁性磁界面迭代反演方法及改進

式（2）表示可將磁界面近似看作由2m×2n個磁化強度相等的垂直磁化直立長方體組合而成。對于變磁性磁界面，每個長方體的磁化強度值不等，則垂直磁化垂直磁異常公式變?yōu)槭剑?）。

其中 Mij為第（i，j）個直立長方體的磁化強度。

同樣采用上述提到的迭代法進行反演即可得到界面深度。式（3）可理解為某一觀測點（x，y，z）處的磁異常值是組成磁界面的所有長方體在點（x，y，z）處產(chǎn)生的磁異常的和。將式（3）改寫成矩陣形式為式（4）。

其中 H為界面深度向量；M 為磁化強度向量；Z為垂直磁化垂直磁異常向量。

在實際計算時，已知Z和M 矩陣，通過矩陣求逆解出H矩陣的結(jié)果，便實現(xiàn)了對界面深度的反演。從這個反演過程不難發(fā)現(xiàn)，若要達到可以將每一個長方體的磁化強度值和高度值看做不變的近似條件，則劃分的長方體數(shù)N必定要足夠大才可行，而隨著N的增大，矩陣計算的工作量也大大增加，這樣就不適宜實際計算的實現(xiàn)，考慮到這個因素，我們對變磁性磁界面反演進行改進，獲得了一種適宜實際計算的方法。

針對這種情況，需要對變磁性磁界面迭代反演方法進行改進?？紤]到以往研究相類似問題時，運用滑動時窗的方法較為便利，故我們也由此入手。若界面的磁化強度值和界面起伏都變化不大的時候，我們可以做進一步近似。取適宜大小的窗口，將窗口處的界面深度值和磁化強度值看做是恒定的，為窗口內(nèi)所有長方體的深度和磁化強度的平均值，在每一個窗口處通過解矩陣求得界面深度值，當(dāng)窗口滑動遍歷整個計算區(qū)域，就得到了整個界面的反演結(jié)果。這樣通過引入滑動窗口的方式，將大大減小反演計算量。

作者采用的變磁性磁界面迭代反演方法，就是建立在這個思想上對常磁性磁界面迭代反演方法改進得到的。首先輸入研究區(qū)化極后的Za磁異常數(shù)據(jù)和事先通過研究區(qū)地質(zhì)資料得到的平均界面深度和磁化強度值；為了減小邊界影響，需要將異常下延到平均界面深度處，計算每一個窗口處對應(yīng)平均界面埋深的磁異常線性正演矩陣的解，得到每個滑動窗口處對應(yīng)的界面深度，這樣通過滑動窗口遍歷整個研究區(qū)域就可以得到整個區(qū)域的界面深度；通過對得到的界面深度值進行磁異常正演與輸入的實際磁異常進行比較，判斷均方誤差是否滿足迭代結(jié)束條件，若不滿足迭代結(jié)束條件，則求取界面深度修正值，返回繼續(xù)進行迭代計算，直到滿足迭代精度或迭代次數(shù)要求后即得到界面深度反演結(jié)果。整個流程見圖2。

在反演過程中，為了減小邊界影響，在計算窗口處矩陣解前，需要先將異常下延到平均界面深度處，這里采用的是加入正則化因子進行空間域向下延拓的方法。

異常場下延可表示為矩陣形式［8－9］：

其中 Z0為已知場值矩陣；A為延拓系數(shù)矩陣；ZH為下延到H深度處平面上場值。

則式（5）可以寫成式（6）的形式。

圖2 變磁性磁界面深度反演流程圖Fig.2 Flowchart of inversion for varying magnetic interface depth

引入一個正則化因子α，有

其中 I為單位矩陣。

由式（7）選取合理α值便可計算得到下延異常，若α的取值過小，則無法較好地控制下延發(fā)散的問題；若α的取值過大，則將導(dǎo)致反演誤差增大。這樣，通過參數(shù)α的引入，較好地完成了磁異常的下延。

2 模型正演驗證

為了驗證方法的可靠性和適用條件，選取模型進行正演計算。作者將選擇規(guī)則球體模型和模擬實際地形的仿真模型分別進行正演計算。

2.1 單一規(guī)則球冠模型正演計算

取球體頂點高度為－4km，平均深度為－4.68 km的單一球冠，測線間距和測點間距均為5km，測線數(shù)和測點數(shù)均為64，界面磁化強度取如圖3（a）所示漸變式，其平均值為10A／m，這里分別取界面磁化強度為真值（變磁性）和為平均值（常磁性）時計算反演結(jié)果，得到磁異常正演結(jié)果、變磁性反演結(jié)果及由常磁性反演結(jié)果見圖3。

為了更清晰地對比反演結(jié)果，取過球冠頂點處的測線數(shù)據(jù)進行誤差分析，剖面反演結(jié)果對比圖見圖4。統(tǒng)計測線數(shù)據(jù)的反演誤差，得到變磁性反演結(jié)果相對誤差為0.2%～2.4%；而常磁性反演結(jié)果相對誤差為0.6%～4.5%。根據(jù)圖3（d）等值線圖顯示，常磁性反演結(jié)果使球缺走勢隨界面磁化強度漸變方向傾斜，而變磁性反演則沒有出現(xiàn)這個問題，同時變磁性磁界面反演的結(jié)果比常磁性磁界面反演結(jié)果相對誤差更小，由此可說明變磁性反演的優(yōu)勢。

圖3 單一球冠模型反演結(jié)果Fig.3 Inversion of single spherical cap model

圖4 球冠模型剖面數(shù)據(jù)對比圖Fig.4 Profile comparison of spherical cap model profile

2.2 仿真模型正演計算

仿真模型模擬的是海盆和海山共存時的海底地形，測線和測點間距均為2km，地形平均高度為－4.25km，其地形等值線圖如圖5（a）所示，區(qū)域磁化強度呈線性漸變，平均值為50A／m，其正演磁異常見圖5（b），變磁性反演結(jié)果見圖5（c）。

同樣為了方便分析方法的效果，取如圖5（c）黑色線所示中心測線上的數(shù)據(jù)分析反演誤差，因測線經(jīng)過測區(qū)地形起伏變化最大的地方，故該測線較具代表性，剖面數(shù)據(jù)對比圖見圖6。分析數(shù)據(jù)得到相對誤差為0.6%～7.4%，其中誤差最大點出現(xiàn)在地形起伏頂點處，可以說該反演方法對地形有一定的“抹平”作用，同時整體測線相對誤差大小較合理，說明針對接近實際情形的仿真模型，該方法也取得了較好的效果。

2.3 影響因素分析

在對幾組模型進行反演驗證的過程中發(fā)現(xiàn)，觀測面高度的選取、測線測點間距的選取和計算范圍的大小也會對反演結(jié)果產(chǎn)生一定的影響，下面通過單一球冠模型對這幾方面影響因素進行簡要說明。

2.3.1 觀測面高度影響。

對上述單一球冠模型固定其形狀，取磁化強度為常值10A／m.，取觀測面高度分別為3km、5km和8km，計算反演結(jié)果。

將得到的三個反演結(jié)果分別與模型地形數(shù)據(jù)進行對比：①當(dāng)觀測面高度為3km時，反演誤差最大為7.1%；②當(dāng)觀測面高度為5km時，反演誤差最大為2.2%；③當(dāng)觀測面高度為8km時，反演誤差最大為1.2%。由此可以看到，隨著觀測面高度的增大，反演誤差逐漸減小。所以觀測面高度取值不宜過小也不宜過大，要根據(jù)研究區(qū)的地質(zhì)情況進行合理選擇。

2.3.2 測線間距影響

當(dāng)測線間距取值不合理時，也會影響反演結(jié)果。為了考察測線間距取值的不同對反演結(jié)果的影響，仍采取單一球缺模型，球缺頂點高度仍為－4km，底界面深度為－5km，磁化強度為10A／m，取測線間距和測點間距相同，分別取測線間距為1km、3 km、5km、10km進行變磁性磁界面反演。

分析反演結(jié)果的相對誤差：①當(dāng)測線間距為1km時，反演結(jié)果誤差最大為10.2%；②當(dāng)測線間距為3km時，反演誤差最大為5.8%；③當(dāng)測線間距為5km時，反演誤差最大為2.2%；④當(dāng)測線間距為10km時，反演誤差為1.1%。

圖6 模擬地形剖面數(shù)據(jù)對比圖Fig.6 Profile comparison of simulative topography

由此得出，測線間距越大，反演誤差越小，但考慮到計算的數(shù)據(jù)量，測線間距也不宜取的過大，要根據(jù)計算區(qū)域的大小，選擇合理的測線間距進行計算。

2.3.3 邊界影響

在計算時應(yīng)將根據(jù)計算數(shù)據(jù)的情況調(diào)整計算區(qū)域大小，不同的計算區(qū)域，其邊界效應(yīng)對反演結(jié)果的影響是不同的。下面仍然固定球缺模型深度、大小、磁化強度值及測線間距值，變化測點測線數(shù)，分別取計算區(qū)域為220km×220km、270km×270km、320km×320km和420km×420km計算反演結(jié)果。

統(tǒng)計反演結(jié)果相對誤差，當(dāng)計算范圍為220km×220時，反演誤差最大為5.4%；當(dāng)計算范圍為270km×270km時，反演誤差最大為3.1%；當(dāng)計算范圍為320km×320km時，反演誤差最大為2.2%；計算范圍增大到420km×420km時，反演誤差最大為1.2%。由此看到，隨著計算范圍的增大，邊界對結(jié)果的影響減小，則反演誤差也逐漸減小。

根據(jù)上面的模型分析及影響因素分析，可以對本文采用的變磁性磁界面迭代反演方法進行以下總結(jié)：

（1）通過規(guī)則形體模型和模擬實際地形地質(zhì)模型的驗證，得到較理想的反演結(jié)果，說明了該方法的可靠性。

（2）采用該方法反演時，對地形存在一定的“抹平”作用，因此對地形起伏較陡峭的區(qū)域不適宜應(yīng)用該方法。

（3）運用變磁性磁界面迭代反演方法的前提是，界面起伏遠小于界面平均深度，并且界面磁化強度不能有正負相差較大的情況，否則會導(dǎo)致反演結(jié)果無法收斂，從而無法達到迭代反演的目的。

（4）在運用該方法時，存在觀測面高度、測線測點間距、計算范圍等幾個主要因素影響反演結(jié)果。因此在對實測數(shù)據(jù)進行反演計算時，應(yīng)充分注意這幾個因素。

（5）觀測面高度取值越大，反演誤差越小，因此應(yīng)該根據(jù)研究區(qū)實際情況，在合理范圍內(nèi)選取最大可取觀測面高度值進行反演。

（6）隨著測線、測點間距取值增大，反演誤差減小，但為了保證計算的數(shù)據(jù)量，測線、測點間距不宜取值過大，應(yīng)根據(jù)計算范圍的大小選擇合理的點距、線距進行反演。

（7）計算范圍越大，邊界影響越小，則反演誤差越小，但為了避免計算量得增大，影響反演計算的速度，也不宜將計算區(qū)域擴的過大，所以在對實測數(shù)據(jù)計算前，應(yīng)該根據(jù)邊界影響大小適當(dāng)調(diào)整計算范圍。

3 實測數(shù)據(jù)磁性基底反演

通過對模型進行反演結(jié)果對比后，對變磁性磁界面迭代反演方法有了一定的了解，將該方法運用于南海東北部某研究區(qū)實測數(shù)據(jù)中，計算該區(qū)域磁性基底深度。研究區(qū)地形沒有起伏較尖銳的點，并且各測點間磁化強度無變化較劇烈的點，適合應(yīng)用變磁性磁界面迭代反演方法進行磁性基底深度反演。

在計算前，需要預(yù)先得到研究區(qū)由磁性基底引起的磁異常及區(qū)域磁化強度值，這里先通過化極及小波分解獲得由磁性基底引起的Za磁異常（圖7），再由研究區(qū)地質(zhì)地球物理資料及相關(guān)方法計算得到區(qū)域磁化強度值（圖8）。將上述數(shù)據(jù)代入變磁性磁界面迭代反演程序中進行計算，得到研究區(qū)磁性基底界面深度分布圖，反演結(jié)果見圖9。圖7～圖9中棕色線段標(biāo)示區(qū)域斷裂，紫色線段標(biāo)示二級構(gòu)造單元，各構(gòu)造單元名稱見圖例。

由圖7－圖9進行對比總結(jié)：

圖7 研究區(qū)化極磁異常小波三階細節(jié)場Fig.7 Wavelet third－order details for Magnetic anomaly of research area

圖8 研究區(qū)基底磁異常Fig.8 Magnetic anomaly for basement of research area

圖9 研究區(qū)磁性基底深度反演結(jié)果Fig.9 Inversion for magnetic basement depth of research area

（1）在圖7中，研究區(qū)Za異常呈條帶狀，近東北－西南方向延展，并出現(xiàn)正磁異常、負磁異常條帶間隔出現(xiàn)的現(xiàn)象，在東沙隆起板緣區(qū)域磁異常呈現(xiàn)最大值約為100nT，在白云坳陷區(qū)域則出現(xiàn)磁異常負值，為－20nT～－100nT。

（2）在圖8中，研究區(qū)磁化強度形態(tài)也呈條帶狀展布，與磁異常延伸方向相近，在圖7中對應(yīng)正磁異常條帶的位置（東沙隆起板緣），區(qū)域磁化強度近乎呈常值，約為2A／m；從磁化強度的這一高值帶的分布特征和強度來看，推論可能是侵入基底的火成巖的反映，其分布位置和展布方向與陽江—統(tǒng)暗沙斷裂東南段非常一致，認為其發(fā)育與這一斷裂有密切關(guān)系；而對應(yīng)磁異常負值的區(qū)域（白云坳陷）磁化強度則出現(xiàn)負值，最小達到－10A／m，并向研究區(qū)東南方向逐漸減小。

（3）圖9中磁性基底在東沙隆起板緣區(qū)域深度近乎不變，約為6.7km，延展區(qū)域也與磁異常形態(tài)相近，向近東北－西南方向延伸，同時對應(yīng)區(qū)域磁異常和磁化強度負值的區(qū)域（白云坳陷）磁性基底深度值則較大，最大約為10km左右；兩部分中間區(qū)域磁性基底深度則較小，為4.5km左右。

從整體反演的磁性基底深度結(jié)果看到，南海東北部地區(qū)除局部地區(qū)外，總體呈現(xiàn)大面積寬緩的等值線特征，一直延伸到洋陸分界處。從磁性基底分布的宏觀特征來看，東北部陸緣的構(gòu)造性質(zhì)更傾向于非火山型［10－15］。

4 結(jié)論

（1）作者使用滑動時窗，對常規(guī)磁界面迭代反演方法進行改進得到適用于變磁性磁界面的迭代反演方法，這是與以往磁界面反演方法不同之處，而應(yīng)用模型進行驗證取得了較好的效果，說明了該方法的可靠性較好。

（2）在運用本文提到的變磁性磁界面迭代反演方法前，需要預(yù)先得到研究區(qū)化極磁異常Za數(shù)據(jù)、磁化強度數(shù)據(jù)和界面平均埋深值。

（3）在對實測數(shù)據(jù)應(yīng)用本文提及的變磁性磁界面反演方法反演時，要注意觀測面高度、邊界影響等因素的影響。觀測面高度取值要合理，計算范圍要根據(jù)實測數(shù)據(jù)正演磁異常情況進行一定調(diào)整，防止邊界效應(yīng)影響反演結(jié)果。

（4）作者研究變磁性磁界面迭代反演方法具有地形“抹平”作用，對于起伏較大的地區(qū)，在地形最值處反演結(jié)果有一定誤差，在應(yīng)用實測數(shù)據(jù)進行計算時要注意這個特點，對于這一缺陷的改進在今后的研究工作中有待解決。

（5）通過實測數(shù)據(jù)計算，得到研究區(qū)磁異常及磁化強度均呈近東北—西南方向條帶狀特征展布，通過變磁性磁界面反演得到的磁性基底分布特征與磁異常特征相似，平均深度為6km左右，最大深度約10km。這為今后地質(zhì)－地球物理進一步的研究提供了一定的依據(jù)。

［1］PARKER R L.The rapid calculation of potential anomalies［J］.G.J.R，1972，37（3）：447－455.

［2］OLDENBURG D M，MOGILIVILAY P R.Generalized subspace methods for 1arge－scale inverse problems［J］.Geophys J Int，1993，114：12－20.

［3］FERNANDO GRUSPI.Nonilerative nonlinear gravity inversion［J］.Geophysics，1993，58（7）：935－940.

［4］BARBOSA V C F，SILVA B C，MEDEIROS W E.Gravity inversion of basement relidf using approximate equality comstrains on depths［J］.Geophysics，1997，62（6）：1745－1757.

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

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

［7］董煥成.重磁勘探教程［M］.北京：地質(zhì)出版社，1993.

［8］王家林，王一新，萬明浩.石油重磁解釋［M］.北京：石油工業(yè)出版社，1991.

［9］管志寧.地磁場與磁力勘探［M］.北京：地質(zhì)出版社，2005.

［10］姚伯初，萬玲.中國南海海域巖石圈三維結(jié)構(gòu)及演化［M］.北京：地質(zhì)出版社，2006.

［11］劉昭蜀，趙煥庭，范時清，等.南海地質(zhì)［M］.北京：科學(xué)出版社，2002.

［12］郝天珧，徐亞，趙百民，等.南海磁性基底分布特征的地球物理研究［J］.地球物理學(xué)報，2009，52：2763－2774.

［13］陳潔，高德章，溫寧，等.南海磁場特征研究［J］.地球物理學(xué)進展，2010，25（2）：376－388.

［14］江凡.南海東北部中生界分布綜合地球物理研究［D］.上海：同濟大學(xué)，2010.

［15］郝天珧，黃松，徐亞，等.南海東北部及臨區(qū)深部結(jié)構(gòu)的綜合地球物理研究［J］.地球物理學(xué)報，2008，51（6）：1785－1796.