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

張向宇，吳健生

（1.廣州海洋地質調查局 ，廣州 510760；2.同濟大學 海洋地質國家重點實驗室，上海 200092）

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

其中 I為單位矩陣。

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

2 模型正演驗證

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

2.1 單一規則球冠模型正演計算

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

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

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

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

2.2 仿真模型正演計算

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

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

2.3 影響因素分析

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

2.3.1 觀測面高度影響。

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

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

2.3.2 測線間距影響

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

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

圖6 模擬地形剖面數據對比圖Fig.6 Profile comparison of simulative topography

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

2.3.3 邊界影響

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

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

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

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

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

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

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

（5）觀測面高度取值越大，反演誤差越小，因此應該根據研究區實際情況，在合理范圍內選取最大可取觀測面高度值進行反演。

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

（7）計算范圍越大，邊界影響越小，則反演誤差越小，但為了避免計算量得增大，影響反演計算的速度，也不宜將計算區域擴的過大，所以在對實測數據計算前，應該根據邊界影響大小適當調整計算范圍。

3 實測數據磁性基底反演

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

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

由圖7－圖9進行對比總結：

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

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

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

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

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

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

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

4 結論

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

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

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

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

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

［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］趙百民，郝天珧.反演磁性地質界面的意義和方法［J］.地球物理學進展，2006，21（2）：353－359.

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

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

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

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

［10］姚伯初，萬玲.中國南海海域巖石圈三維結構及演化［M］.北京：地質出版社，2006.

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

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

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

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

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