Parker法視磁化強度填圖方法研究及應用

秦 熠，楊宇山，劉天佑，朱 丹

(1.中國地質(zhì)大學 地球物理與空間信息學院，湖北 武漢 430074；2.四川省煤田地質(zhì)局 地質(zhì)工程勘察設(shè)計研究院，四川 成都 610000)

1 引 言

視磁化強度填圖方法是一種利用磁異常反演地下地質(zhì)體磁性強弱的方法，利用航磁資料進行視磁化強度填圖可以有效確定磁性體邊界，劃分磁性地層，特別對于存在大規(guī)模隱伏巖體地區(qū)具有良好的應用效果，作為研究區(qū)域構(gòu)造的輔助手段，已被廣泛地應用于地質(zhì)構(gòu)造解釋與研究之中[1]。

1973年P(guān)arker[2]將傅里葉變化引入位場正演計算，提出了連續(xù)單界面模型快速正演方法，即Parker法。隨后，Oldenburg[3](1974)在Parker法基礎(chǔ)上推導出了頻率域密度界面迭代反演方法，由于該方法解決了過去空間域離散模型界面反演方法計算速度慢的問題，因此被廣泛地應用于起伏界面的計算。王萬銀等[4](1993)在單界面重力場正反演基礎(chǔ)之上研究了更加符合實際情況的雙界面重力場正反演問題，并提出了保證反演快速收斂的措施。

視磁化強度填圖興起于20世紀80年代，Pilkington[5](1989)根據(jù)Parker界面正反演方法推導出了磁化強度反演公式，并應用于Athabasca盆地。早期視磁化強度填圖，尤其在化極處理之后存在負值的問題，難以解釋。管志寧等[6](1990)研究了視磁化強度填圖出現(xiàn)負值的原因，并提出一種采用立體角滑動窗口迭代反演方法的視磁化強度填圖，應用于秦巴地區(qū)。這種方法適用于地形起伏不大，觀測平面水平的航磁資料處理。張貴賓等[7](1996)從理論與模型算例兩方面對Parker公式為基礎(chǔ)的頻率域位場反演方法做了深入研究，提出新的快速正則化反演方法，并將其推廣為單界面無限延深視密度和視磁化強度填圖方法。Pilkington[8](1998)提出了等效地形改正方法，實現(xiàn)了起伏地形起伏觀測面視磁化強度填圖，但是并未解決磁化強度存在負值的問題。眭素文[9](2004)結(jié)合前人無負值填圖與等效地形改正兩種方法，實現(xiàn)了小區(qū)域中高山區(qū)起伏測量面航磁資料無負值視磁化強度填圖。此后，大量學者將視磁化強度填圖方法應用于巖體與構(gòu)造研究中，朱曉穎[10](2007)利用視磁化強度填圖方法研究了內(nèi)蒙古北山地區(qū)斷裂構(gòu)造，圈定侵入巖體，深化了研究區(qū)區(qū)域地質(zhì)、成礦規(guī)律和控礦特征的認識。劉天佑等[11](2012)運用密度填圖反演推斷了西藏朗縣秀溝工區(qū)蛇紋石化橄欖巖的范圍以及鉻鐵礦與礦化體的重磁遠景。王賽昕等[12](2017)針對攀枝花地區(qū)進行視磁化強度填圖，利用磁化強度與全鐵質(zhì)量分數(shù)的經(jīng)驗公式獲得了區(qū)內(nèi)全鐵質(zhì)量分數(shù)分布圖。

盡管前人對起伏界面模型視磁化強度填圖方法進行了大量研究，但是受到方法以及資料的限制，在實際應用時仍然主要采用水平界面模型進行反演，對于區(qū)域跨度廣地形起伏大的地區(qū)會帶來較大誤差。針對于這個問題，本文提出一種功率譜分析約束下的雙界面起伏模型Parker法視磁化強度填圖新思路，步驟如下：①根據(jù)化極磁異常形態(tài)，將等軸狀或長軸狀孤立正磁異常單獨分區(qū)；②分別截取每一個小分區(qū)磁異常以及研究區(qū)全區(qū)異常，采用功率譜分析方法反演得到每一個小分區(qū)的上界面與下界面埋深以及全區(qū)平均上下界面深度；③根據(jù)反演得到的小分區(qū)上下界面深度構(gòu)建上界面與下界面起伏形態(tài)，其中負磁異常等未分區(qū)域深度取全區(qū)界面平均深度；④以功率譜分析構(gòu)建的上下起伏界面為約束，采用雙界面Parker視磁化強度填圖法反演得到更加符合實際的磁性分布。

2 起伏界面視磁化強度填圖原理

Parker法頻率域磁性界面正反演方法是一種已知磁化強度參數(shù)，反演界面深度的方法。由Parker法磁性界面正反演公式可知：

(1)

由(1)式可知，若已知均勻密度分布值M，則由(1)式可以計算界面起伏h，也可以把(1)式推廣為已知界面起伏h，求磁化強度分布，這就是Parker法視磁化強度填圖的方法原理。假設(shè)地下磁性層磁性界面深度已知，需要反演求解界面之間磁化強度的水平分布，其中磁化強度M是隨水平位置變化的函數(shù)。類似于Parker法磁性界面深度反演方法中將界面起伏看作平均深度與相對于平均深度起伏兩部分組成一樣，把磁化強度也看作是平均磁化強度與相對平均磁化強度的變化量組成，即：

M=M0+ΔM

(2)

那么磁異常也可以看作是由平均磁化強度產(chǎn)生的磁異常與相對平均磁化強度的變化量產(chǎn)生的磁異常兩部分組成：

(3)

其中：

(4)

(5)

(6)

3 理論模型研究

考慮實際磁性體延深不同的情況，本文建立厚層與薄層巖體模型，分析不同上界面與下界面取值對視磁化強度反演的影響；然后，采用功率譜分析法反演理論模型的上下界面深度，以得到的界面深度約束反演磁化強度，比較其與理論模型的誤差，驗證功率譜分析法結(jié)合視磁化強度填圖反演的有效性；最后，建立理論模型模擬地下不同磁化強度巖體分布，比較起伏界面模型視磁化強度填圖與水平界面模型視磁化強度填圖應用效果。

3.1 界面深度對填圖的影響分析

通過單個直立長方體模型建立厚層巖體與薄層巖體理論模型，分析上下界面取值對反演視磁化強度的影響，理論模型正演參數(shù)選取見表1，正演磁異常如圖1所示。視磁化強度填圖計算時取磁性界面水平(圖2)，分別繪制厚、薄層巖體深度變化率與相對誤差曲線圖，橫坐標為相對界面深度變化率，其中上界面(圖3中藍線)深度變化率為Δh(選取上界面與理論上界面深度差)與理論上界面深度H的比值，下界面(圖3中紅線)深度變化率為Δl(選取下界面與理論下界面深度差)與理論模型延深長度L的比值，縱坐標為磁性體分布范圍內(nèi)反演的磁化強度均值與理論磁化強度的相對誤差。對于上界面深度影響，固定理論下界面深度不變，只改變上界面深度，然后進行視磁化強度填圖計算。當選取上界面深度位于理論上界面之上時，其深度變化率為正。對于下界面深度影響，則固定理論上界面深度不變，只改變下界面深度，再進行視磁化強度填圖計算。當選取下界面深度位于理論下界面之上時，其深度變化率為負。

綜合分析界面深度變化與相對誤差曲線圖，可以得到以下幾點認識：①不論上界面還是下界面，相同界面深度變化下薄層巖體的磁化強度相對誤差始終大于厚層巖體，說明上下界面深度取值對薄層巖體影響更大，對厚層巖體影響相對較小；②不論薄層巖體還是厚層巖體，在30 %深度變化范圍內(nèi)，相同界面深度變化下上界面相對誤差要大于下界面相對誤差，表明上界面深度取值對填圖結(jié)果影響要大于下界面深度取值；③當上界面位于理論上界面之上或下界面位于理論下界面之下時，曲線斜率逐漸變小，相對誤差變化隨深度變化率的增大逐漸變小，即深度變化影響逐漸減弱。當上界面位于理論上界面之下或下界面位于理論下界面之上時，曲線斜率逐漸變大，深度變化影響逐漸增強。

理論模型表明，采用實際資料進行視磁化強度填圖時，磁性上界面深度選取至關(guān)重要，其對反演計算的磁化強度值影響較大。此外，薄層巖體由于其延深較小，其下磁性界面深度選取不準也會造成反演計算較大誤差。

表1 厚層與薄層巖體模型參數(shù)

圖1 理論模型正演磁異常Fig.1 Theoretical model positive magnetic anomaly

圖2 界面深度選取示意圖Fig.2 Interface depth selection diagram

圖3 界面深度對填圖結(jié)果影響Fig.3 Interface depth affects the mapping results 注：當選取界面較理論界面深時，上界面深度變化率為負，下界面深度變化率為正

3.2 功率譜分析法確定界面深度效果對比

利用航空磁測資料進行視磁化強度填圖時，若能較為準確確定磁性體上下界面深度就能取得良好的反演效果。目前常用的反演磁性體深度方法[13]有歐拉齊次方程法和功率譜分析法、矩譜法等。歐拉齊次方程法不需要預先知道反演場源的形狀和參數(shù)，能夠自動快速估算場源埋深，但是很難確定下界面深度，而功率譜分析法利用場源在對數(shù)功率譜曲線上表現(xiàn)出不同梯度的線性下降特征反演場源深度，不僅能夠反演磁性體上界面深度，也可以反演磁性體底界面深度，結(jié)合功率譜分析對磁性界面的約束與視磁化強度填圖法可以得到更加符合實際情況的反演結(jié)果。此外，對于孤立形體，為了準確地確定孤立形體的下界面深度，可以采用矩譜法求質(zhì)心，再通過上界面深度與質(zhì)心深度獲得下界面深度。

采用上述厚層與薄層巖體模型，使用功率譜分析法反演磁性體的上下界面深度，低頻段反映深部界面信息，中高頻段反映淺部界面信息。根據(jù)厚層巖體功率譜曲線(圖4a)得到模型上界面深度為37.1 m，絕對誤差7.1 m，下界面深度140.8 m，絕對誤差10.8 m，然后采用反演得到的界面深度計算厚層巖體邊界范圍內(nèi)平均磁化強度為8.49 A/m，與理論磁化強度的相對誤差為6 %。由薄層巖體功率譜曲線(圖4b)得到模型上界面深度為35.1 m，絕對誤差5.1 m，下界面深度61.3 m，絕對誤差11.3 m，反演得到薄層巖體平均磁化強度為7.16 A/m，與理論磁化強度的相對誤差為11 %(見表2)。

比較反演計算值與理論值，功率譜分析法反演的上界面深度要比下界面深度更準確，但是反演得到的界面深度都要比實際界面深度更深一些。對于厚層巖體來說，采用此方法反演的上下界面誤差相較其延深來說影響較小，因此反演得到的磁化強度值也更加接近真實值。但是對于薄層巖體來說，根據(jù)該方法確定的上下界面誤差相對其延深來說影響較大，造成其視磁化強度填圖相對誤差要大于厚層巖體?？偟膩碚f，運用功率譜分析法能夠較為有效地確定磁性體上下界面深度，以此為約束得到的視磁化強度填圖相對誤差不大于15 %，理論模型應用效果較好。

圖4 功率譜分析Fig.4 Power spectrum analysis

3.3 起伏界面模型與水平界面模型填圖效果對比

表3 理論模型參數(shù)

圖5 不同界面模式填圖效果對比Fig.5 Comparison of different interface mode mapping effects

以不同埋深直立長方體模型模擬實際磁性巖體分布(圖5a)，對比水平界面模型與起伏界面模型視磁化強度反演不同磁性巖體的效果。其中1號模型為厚層弱磁性巖體，磁化強度為2 A/m，2號模型為薄層強磁性巖體，磁化強度為5 A/m，3號模型為厚層中等磁性巖體，磁化強度為3.5 A/m，各磁性體上界面與下界面深度不同，模型參數(shù)取值見表3。從正演磁異常(圖5b)來看，高磁異常范圍與模型分布范圍相對應，其中1號模型磁異常略大于2號模型，異常幅值在500～550 nT之間，3號模型磁異常最強，異常幅值在650～700 nT之間，正演磁異常大小與理論模型磁化強度并無直接對應關(guān)系。

采用水平界面模型計算時，使用功率譜分析法由全區(qū)磁異常得到上界面平均深度30 m，下界面平均深度70 m，選取該參數(shù)作為水平界面模型約束進行視磁化強度填圖，反演結(jié)果見圖5(e)。當采用起伏界面模型計算時，按照區(qū)內(nèi)高磁異常個數(shù)將全區(qū)分為3個異常區(qū)，功率譜分析法分別反演每一個異常區(qū)域內(nèi)頂?shù)捉缑媛裆?，分區(qū)之外的區(qū)域頂?shù)捉缑嫔疃葎t由全區(qū)磁異常功率譜分析求得，然后將反演得到的頂界面與底界面深度擬合連接成面得到磁性層上下界面起伏形態(tài)(圖5c，圖5d)，以得到的起伏界面模型深度為約束對正演磁異常做視磁化強度填圖結(jié)果(圖5f)。

由水平界面模型填圖結(jié)果可知,1號與2 號巖體磁化強度相似，磁化強度值在2.8～3.8 A/m之間，屬于中等磁性巖體，3號巖體磁化強度值略高于1、2號巖體，磁化強度值在3.2～4.2 A/m之間，屬于中高磁性巖體。顯然該結(jié)果與理論模型存在較大偏差，未能正確區(qū)分巖體磁性的相對強弱，計算的磁化強度值也并無對應關(guān)系。對比起伏界面模型反演結(jié)果可知1號巖體磁性相對較弱，磁化強度分布在2～2.8 A/m，2號巖體磁性最強，磁化強度在4～5 A /m，3號巖體磁性強于1號巖體，弱于2號巖體，磁化強度在3～4 A/m。反演計算的磁化強度值與理論模型磁化強度值相對應，能夠較為明確地劃分為弱、強、中等磁性三類不同磁化強度巖體。從理論模型反演結(jié)果來看，起伏界面模型相較水平界面模型不論反演計算的磁化強度值還是不同巖體磁性相對強弱關(guān)系都與理論模型更相似，說明結(jié)合功率譜分析法反演起伏界面深度約束下的視磁化強度填圖能夠更為準確地圈定巖體分布范圍，劃分不同磁化強度巖體。

4 小秦嶺地區(qū)的應用

4.1 研究區(qū)地質(zhì)概況

小秦嶺地區(qū)是秦嶺成礦帶地質(zhì)研究、資源勘查的主要工作區(qū)之一，區(qū)內(nèi)已查明十多處大型鉬礦、金礦等金屬礦床，是我國金、銀、鉬、鉛等多金屬礦的重要產(chǎn)地[14]。研究區(qū)位于華北陸塊南緣東秦嶺鉬礦帶以西，出露地層主要為太古宇太華群變質(zhì)巖系、中元古界熊耳群火山巖系和中元古界官道口群、新元古界欒川群、下古生界陶灣群等一套碳酸鹽巖組合以及中新生代陸相盆地沉積。受板塊間深大斷裂和秦嶺造山帶長期活動影響，區(qū)內(nèi)構(gòu)造復雜，斷裂與褶皺發(fā)育，區(qū)域構(gòu)造線主要以東西向展布，斷裂構(gòu)造可以分為近東西向、北東向和北西向三組。該區(qū)經(jīng)歷元古代-古生代-中生代多期巖漿活動，巖漿巖十分發(fā)育(圖6)，其中主要以中生代和古生代酸性巖為主，地表有較大面積出露。研究表明，區(qū)內(nèi)巖漿巖與鉬礦、金礦等多金屬礦床形成有關(guān)，已查明金屬礦床大多分布在巖體內(nèi)部或其邊部附近[15]。

4.2 地球物理特征

由于不同時期巖漿侵入作用與火山活動強弱差異較大，造成了地層與巖石不同的磁性特征。前人物性統(tǒng)計資料表明(表4)，秦嶺地區(qū)大部分沉積巖地層磁性較弱，只有太古代及元古代地層作為古老的結(jié)晶基底具有較高的磁性，但極不均勻，磁性變化大[16]。區(qū)內(nèi)巖漿巖磁性特征復雜，不同時期存在明顯差異(表5)。早古生代時期形成的巖漿巖磁性最高，中生代次之，古生代最低，其中中生代二長花崗巖和閃長花崗巖較其他時期形成的花崗巖磁性更強，隨著巖漿巖從酸性到基性的變化，巖石磁性逐漸增強。

1.中生代酸性巖;2.古生代酸性巖;3.中性巖;4.基性巖;5.金礦點;6.鉬礦點;7.鐵銅鉛礦點;8.斷裂 圖6 小秦嶺巖漿巖分布圖(修改自《秦嶺成礦帶地質(zhì)礦產(chǎn)圖》李宗會2013)Fig.6 Map of Xiaoqinling magmatic rocks

表4 秦嶺地區(qū)地層磁性

表5 秦嶺地區(qū)巖漿巖磁性

本文采用秦嶺東段地區(qū)1∶20萬航空磁測資料進行處理解釋，平均飛行高度500 m。航磁異常(圖7a)表現(xiàn)為東西向和北東向連續(xù)條帶狀異常，顯示出小秦嶺地區(qū)受華北陸塊與揚子陸塊擠壓形成的東西向的線性構(gòu)造特征。經(jīng)化極處理消除斜磁化影響(圖7b)，磁異常極值北移，高磁異常分布范圍與地表大面積出露的燕山期花崗巖對應關(guān)系較好；印支期花崗巖與沉積性地層磁性較弱，對應負磁異常或低值正磁異常；具有一定規(guī)模的中基性巖能夠產(chǎn)生較強的磁異常，而小規(guī)模條帶狀中基性巖產(chǎn)生磁異常較弱。

1.中生代酸性巖;2.古生代酸性巖;3.中性巖;4.基性巖;5.金礦點;6.鉬礦點;7.鐵銅鉛礦點 圖7 小秦嶺地區(qū)磁異常Fig.7 Magnetic anomalies in the Xiaoqinling area

4.3 小秦嶺地區(qū)視磁化強度填圖應用

由于區(qū)內(nèi)鉬、金、鉛等金屬礦床的形成與巖體關(guān)系密切，尋找隱伏巖體對于尋找金屬礦具有重要指導意義。研究區(qū)內(nèi)巖體分布，首先要提取淺部磁性體異常特征，由于小波多尺度分解對高頻信號分辨率高的特點[17]，本文采用小波多尺度分解來分離磁異常，提取的剩余異常見圖8(a)。

1.中生代酸性巖；2.古生代酸性巖；3.中性巖；4.基性巖；5.金礦點；6.鉬礦點；7.鐵銅鉛礦點； 8.磁異常小分區(qū)；9.推測隱伏中生代花崗巖(C1、C2、C3為出露巖體；Y1、Y2、…Y11為隱伏巖體) 圖8 小秦嶺視磁化強度填圖Fig.8 Magnetization map in Xiaoqinling area

由理論模型可知，對孤立磁異常分區(qū)，然后采用功率譜分析法確定磁性體頂?shù)捉缑孢M行起伏界面視磁化強度填圖比傳統(tǒng)水平界面填圖效果更優(yōu)。本文將全區(qū)剩余磁異常劃分為24個區(qū)域，分別對每一個孤立磁異常做功率譜分析，計算可以得到其對應的上下界面埋深。而分區(qū)之外的區(qū)域上下界面深度則由全區(qū)磁異常功率譜分析得到，反演平均上界面深度為1 462 m，下界面深度為8 642 m。此外，由于該地區(qū)巖漿巖體出露面積較大，為了提高起伏上界面選取精度，可以結(jié)合GPS地形高程數(shù)據(jù)來確定巖體出露地區(qū)起伏上界面的深度，GPS地形相對觀測面轉(zhuǎn)換的界面深度由航空磁測平均飛行高度來確定(圖8b)。即研究區(qū)起伏上界面深度由兩種方法共同獲取，在巖體未出露的分區(qū)，上界面深度由功率譜分析法得到；在巖體出露的分區(qū)，上界面深度直接由GPS地形高程數(shù)據(jù)獲取，得到的起伏上下界面如圖8c、圖8d。

以反演計算的起伏上下界面為約束，對剩余磁異常采用雙界面模型進行填圖，反演的得到磁化強度分布見圖8(e)。不同類型巖體之間存在明顯差異，古生代酸性巖和沉積地層大多分布在中低磁化強度區(qū)域，中性巖與基性巖本身規(guī)模較小，主要分布在中等磁化強度區(qū)域，中生代酸性巖則分布在中高磁化強度區(qū)域。蔡川-三要-官坡中高磁化強度條帶區(qū)與地表大面積出露的東西向中生代花崗巖對應關(guān)系較好，反映了中生代花崗巖磁性相對較高的特征。華陰-潼關(guān)-靈寶高磁化強度帶也與地表斷續(xù)出露的中生代花崗巖對應，推測地層之下可能存在隱伏的中生代花崗巖帶。五里川北部北西向長軸狀高磁化強度區(qū)地表出露較大面積的中性閃長巖以及小面積中生代花崗巖，由中性巖填圖磁性特征推測，中性巖之下可能隱伏有較大規(guī)模中生代花崗巖體。此外，由中生代酸性巖高磁化強度特征，共圈定區(qū)內(nèi)11個隱伏中生代花崗巖體(圖8f)。對比分析研究區(qū)內(nèi)目前已查明多金屬礦床分布位置，鉬、金礦床與中生代花崗巖關(guān)系密切，集中分布在北西區(qū)域出露的中生代花崗巖體(C-1)邊部，其他區(qū)域鉬、金礦床相對較少，大多分散分布于圈定的隱伏中生代花崗巖體邊部(Y-1、Y-2、Y-3、Y-4、Y-10)。隱伏中生代花崗巖體與圍巖接觸部位既為成礦有利區(qū)域，成礦可能性較高，可以作為后期找礦工作開展的重點方向。

5 結(jié) 論

本文以Parker法界面正反演公式為基礎(chǔ)，推導了起伏界面模型視磁化強度填圖迭代公式；針對地形變化劇烈地區(qū)提出一種起伏雙界面模型快速視磁化強度填圖的新思路，即結(jié)合功率譜分析法確定磁性起伏界面，然后以計算的起伏界面來約束視磁化強度填圖；理論模型驗證了該方法的有效性，相較傳統(tǒng)水平界面模型提高了反演精度。將該視磁化強度填圖方法應用于小秦嶺多金屬成礦帶，圈定了區(qū)內(nèi)與鉬、金成礦密切相關(guān)的隱伏中生代花崗巖體11個，其中Y-1、Y-2、Y-3、Y-4、Y-10隱伏巖體邊部均存在已查明的鉬礦或金礦床與之對應，剩余未發(fā)現(xiàn)礦床點的隱伏巖體Y-5～Y-9、Y-11邊部位置成礦可能性較高，這些巖體邊部附近即為成礦有利位置，可以為今后金屬礦產(chǎn)勘查工作提供參考。需要指出的是，本文僅采用航空磁測資料以及地質(zhì)資料等作為推斷依據(jù)，采用資料較為單一，得到的結(jié)論有待更進一步驗證。此外，建議嘗試將功率譜法求下界面深度改用矩譜法，先求質(zhì)心再求下界面深度。