本發(fā)明涉及地球物理勘測技術(shù)領(lǐng)域，具體為磁法勘探領(lǐng)域，特別是涉及一種磁異?；瘶O方法。

背景技術(shù)：

在磁法勘探中，磁測數(shù)據(jù)的處理和解釋(如化極、導(dǎo)數(shù)換算、分量換算、正演和反演等)通常需要已知磁場方向信息。磁場是一個(gè)矢量場，而利用磁力儀測得的磁場強(qiáng)度是標(biāo)量，沒有磁場方向。磁場方向是感應(yīng)磁化方向和剩余磁化方向的矢量和。利用國際地磁參考場獲得感應(yīng)磁場方向，即地磁場方向；而剩余磁化方向無法直接獲取。若磁源地質(zhì)體剩磁不明顯，那么磁場方向與感應(yīng)磁場方向基本一致；否則，較強(qiáng)的剩磁會(huì)使磁場方向與感磁方向出現(xiàn)很大偏差，使得地下磁性體產(chǎn)生斜磁化現(xiàn)象。此時(shí)，若直接利用地磁場方向作為磁場方向進(jìn)行磁異常化極，以及在此基礎(chǔ)上進(jìn)行的一系列磁異常數(shù)據(jù)處理與轉(zhuǎn)換、正演和反演，其結(jié)果都無法反映真實(shí)的地質(zhì)情況。

目前，除采集巖石標(biāo)本并測定剩磁參數(shù)外，國內(nèi)外學(xué)者提出了一些估計(jì)剩磁方向的方法。唐俊德(1986)提出利用磁異常三分量確定磁化方向的圓曲線法；Roest和Pilkington(1993)提出基于磁場總梯度模和磁源重力異常水平梯度互相關(guān)的估計(jì)方法；Medeiros和Silva(1995)提出利用等效源磁矩反演估計(jì)磁化方向；甘西(2001)利用磁場的垂直分量確定磁性體磁化方向；Phillips(2005)給出估計(jì)磁化方向的Helbig(1963)積分法的直接算法和間接算法；Dannemiller和Li(2006)提出基于化極磁異常垂直梯度與總梯度?；ハ嚓P(guān)的估計(jì)方法；Gerovska等(2009)提出基于化極磁異常與總模量磁異?；ハ嚓P(guān)的估計(jì)方法。

但是，磁異?；瘶O的準(zhǔn)確性還與研究區(qū)的緯度變化有關(guān)。若研究區(qū)的緯度范圍較大，地磁傾角變化也會(huì)較大，斜磁化的程度會(huì)隨緯度不同而變化。若采用單一的磁傾角作化極處理顯然是不合理的，其結(jié)果將無法反映真實(shí)的地質(zhì)情況，降低磁法勘探的準(zhǔn)確度。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明主要解決的技術(shù)問題是提供一種磁異?；瘶O方法，能夠利用地磁場方向和剩磁方向?qū)?shí)測磁異常進(jìn)行變傾角化極，提高磁法勘探的準(zhǔn)確度。

為解決上述技術(shù)問題，本發(fā)明采用的一個(gè)技術(shù)方案是：提供一種磁異?；瘶O方法，包括以下步驟：采集研究區(qū)的磁場數(shù)據(jù)，對采集到的磁場數(shù)據(jù)進(jìn)行校正，得到磁異常數(shù)據(jù)；將磁異常數(shù)據(jù)網(wǎng)格化，形成規(guī)則的磁異常數(shù)據(jù)t(x,y)，其中x,y分別表示磁異常數(shù)據(jù)的平面網(wǎng)格點(diǎn)坐標(biāo)，記為矩陣T(1:m,1:n)，大小為m×n，m表示行數(shù)，n表示點(diǎn)數(shù)；獲得研究區(qū)的地磁場方向，并將地磁場方向整理成與矩陣T(1:m,1:n)等大的矩陣，記為I(1:m,1:n)；對矩陣T(1:m,1:n)進(jìn)行鑲邊處理，將矩陣T(1:m,1:n)轉(zhuǎn)變成行數(shù)和列數(shù)均為2的最小整數(shù)冪的新矩陣T’(1-m’:m+m’,1-n’:n+n’)；設(shè)計(jì)頻率域化極因子，并對頻率域化極因子做反傅里葉變換，得到化極結(jié)果；計(jì)算化極結(jié)果的垂直梯度和總梯度，并計(jì)算其互相關(guān)系數(shù)，以最大互相關(guān)系數(shù)確定剩磁方向；根據(jù)地磁場方向和剩磁方向，并根據(jù)頻率域化極因子對矩陣T’和I在頻率域內(nèi)進(jìn)行變傾角化極處理，計(jì)算研究區(qū)磁異常的變傾角化極結(jié)果。

其中，將磁異常數(shù)據(jù)網(wǎng)格化的步驟包括：利用kriging插值方法將磁異常數(shù)據(jù)網(wǎng)格化。

其中，將地磁場方向整理成與矩陣T(1:m,1:n)等大的矩陣，記為I(1:m,1:n)的步驟包括：

獲取地磁傾角范圍為a-b；

設(shè)置矩陣

利用矩陣A的偽逆矩陣計(jì)算磁場方向；

P＝pinv(A)×V，pinv表示偽逆矩陣P為計(jì)算磁場方向時(shí)的參數(shù)矩陣，其使用方法見下式；

設(shè)置網(wǎng)格坐標(biāo)陣[x,y]＝meshgrid(1:n,1:m)；

計(jì)算矩陣I＝P(1)+P(2)×y+P(3)×x。

其中，對矩陣T(1:m,1:n)進(jìn)行鑲邊處理的步驟包括：

采用以下公式進(jìn)行鑲邊：

其中，T位于T’的中心位置，即T(i,j)＝T’(i,j)，(1≤i≤m,1≤j≤n)。

其中，頻率域化極因子為：

式中，ω_x和ω_y是x、y方向的圓波數(shù)，是地磁方向的單位向量，L、M、N是地磁場的方向余弦：cosIcosD,cosIsinD,sinI，I和D分別表示地磁場傾角和地磁場偏角；是剩磁方向的單位向量，l,m,n分別是剩磁方向的方向余弦：cosicosd,cosisind,sini，i和d分別表示剩磁傾角和剩磁偏角。

其中，對頻率域化極因子做反傅里葉變換的步驟包括：

采用以下的公式做反傅里葉變換：

其中，RTP表示化極結(jié)果。

其中，根據(jù)頻率域化極因子對矩陣T’、I在頻率域內(nèi)進(jìn)行變傾角化極處理的步驟之前包括：

從矩陣I中提取平均地磁傾角為平均地磁方向；

計(jì)算出研究區(qū)中的實(shí)際地磁方向與平均地磁方向的差

其中，利用泰勒公式進(jìn)行變傾角化極，采用的公式是：

式中，RTP_diff是變傾角化極結(jié)果，RTP_mean是使用研究區(qū)平均地磁傾角和地磁偏角得到的化極結(jié)果。

本發(fā)明的有益效果是：區(qū)別于現(xiàn)有技術(shù)的情況，本發(fā)明提供一種磁異?；瘶O方法，包括以下步驟：采集研究區(qū)的磁場數(shù)據(jù)，對采集到的磁場數(shù)據(jù)進(jìn)行校正，得到磁異常數(shù)據(jù)；將磁異常數(shù)據(jù)網(wǎng)格化，形成規(guī)則的磁異常數(shù)據(jù)t(x,y)，其中x,y分別表示磁異常數(shù)據(jù)的平面網(wǎng)格點(diǎn)坐標(biāo)，記為矩陣T(1:m,1:n)，大小為m×n，m表示行數(shù)，n表示點(diǎn)數(shù)；獲得研究區(qū)的地磁場方向，并將地磁場方向整理成與矩陣T(1:m,1:n)等大的矩陣，記為I(1:m,1:n)；對矩陣T(1:m,1:n)進(jìn)行鑲邊處理，將矩陣T(1:m,1:n)轉(zhuǎn)變成行數(shù)和列數(shù)均為2的最小整數(shù)冪的新矩陣T’(1-m’:m+m’,1-n’:n+n’)；設(shè)計(jì)頻率域化極因子，并對頻率域化極因子做反傅里葉變換，得到化極結(jié)果；計(jì)算化極結(jié)果的垂直梯度和總梯度，并計(jì)算其互相關(guān)系數(shù)，以最大互相關(guān)系數(shù)確定剩磁方向；根據(jù)地磁場方向和剩磁方向，并根據(jù)頻率域化極因子對矩陣T’和I在頻率域內(nèi)進(jìn)行變傾角化極處理，計(jì)算研究區(qū)磁異常的變傾角化極結(jié)果。因此，本發(fā)明能夠利用地磁場方向和剩磁方向?qū)?shí)測磁異常進(jìn)行變傾角化極，提高磁法勘探的準(zhǔn)確度。

附圖說明

圖1是本發(fā)明實(shí)施例提供的一種磁異?；瘶O方法的流程圖；

圖2為南充-旺蒼地區(qū)的磁異常資料；

圖3為南充-旺蒼地區(qū)的磁異常垂直梯度與總梯度的互相關(guān)系數(shù)；

圖4為南充-旺蒼地區(qū)考慮剩磁的變傾角化極磁異常；

圖5為石柱-大竹地區(qū)的磁異常資料；

圖6為石柱-大竹地區(qū)的磁異常垂直梯度與總梯度的互相關(guān)系數(shù)；

圖7為石柱-大竹地區(qū)考慮剩磁的變傾角化極磁異常。

具體實(shí)施方式

請參閱圖1，圖1是本發(fā)明實(shí)施例提供的一種磁異常化極方法的流程圖。如圖1所示，本實(shí)施例的方法包括以下步驟：

步驟S1：采集研究區(qū)的磁場數(shù)據(jù)，對采集到的磁場數(shù)據(jù)進(jìn)行校正，得到磁異常數(shù)據(jù)。

本步驟中，具體是通過磁力儀在野外采集研究區(qū)的磁場數(shù)據(jù)，并將獲得的磁場數(shù)據(jù)導(dǎo)入計(jì)算機(jī)中，以進(jìn)行各種校正，以得到異常數(shù)據(jù)。

步驟S2：將磁異常數(shù)據(jù)網(wǎng)格化，形成規(guī)則的磁異常數(shù)據(jù)t(x,y)，其中x,y分別表示磁異常數(shù)據(jù)的平面網(wǎng)格點(diǎn)坐標(biāo)，記為矩陣T(1:m,1:n)，大小為m×n，m表示行數(shù)，n表示點(diǎn)數(shù)。

本步驟中，具體是通過surfer軟件中的kriging插值方法將磁異常數(shù)據(jù)網(wǎng)格化。

步驟S3：獲得研究區(qū)的地磁場方向，并將地磁場方向整理成與矩陣T(1:m,1:n)等大的矩陣，記為I(1:m,1:n)。

本步驟中，具體是利用國際地磁參考場，獲得研究區(qū)的地磁場方向。本發(fā)明選用的國際地磁參考場(international geomagnetic reference field,IGRF)優(yōu)選是第11代國際地磁參考場(IGRF-11)，由國際地磁學(xué)與高空物理學(xué)聯(lián)合會(huì)(IAGA)于2009年12月發(fā)布。

其中，本步驟的將地磁場方向整理成與矩陣T(1:m,1:n)等大的矩陣，記為I(1:m,1:n)的步驟具體包括如下：首先獲取地磁傾角范圍為a-b，然后設(shè)置矩陣進(jìn)而利用矩陣A的偽逆矩陣計(jì)算磁場方向；

P＝pinv(A)×V，pinv表示偽逆矩陣P為計(jì)算磁場方向時(shí)的參數(shù)矩陣，進(jìn)一步設(shè)置網(wǎng)格坐標(biāo)陣[x,y]＝meshgrid(1:n,1:m)，最后計(jì)算矩陣I＝P(1)+P(2)×y+P(3)×x。

步驟S4：對矩陣T(1:m,1:n)進(jìn)行鑲邊處理，將矩陣T(1:m,1:n)轉(zhuǎn)變成行數(shù)和列數(shù)均為2的最小整數(shù)冪的新矩陣T’(1-m’:m+m’,1-n’:n+n’)。

本步驟中，優(yōu)選以“余弦衰減置零擴(kuò)邊法”給矩陣T(1:m,1:n)鑲邊。具體采用以下公式進(jìn)行鑲邊：

其中，T位于T’的中心位置，即T(i,j)＝T’(i,j)，(1≤i≤m,1≤j≤n)。

本步驟中，利用“余弦衰減置零擴(kuò)邊法”給數(shù)據(jù)進(jìn)行鑲邊，很好的壓制了邊界效應(yīng)，提高了數(shù)據(jù)處理的精度。

步驟S5：設(shè)計(jì)頻率域化極因子，并對頻率域化極因子做反傅里葉變換，得到化極結(jié)果。

本步驟中，頻率域化極因子具體為：

式中，ω_x和ω_y是x、y方向的圓波數(shù)，是地磁方向的單位向量，L、M、N是地磁場的方向余弦：cosIcosD,cosIsinD,sinI，I和D分別表示地磁場傾角和地磁場偏角；是剩磁方向的單位向量，l,m,n分別是剩磁方向的方向余弦：cosicosd,cosisind,sini，i和d分別表示剩磁傾角和剩磁偏角。

其中，具體是采用以下的公式做反傅里葉變換：

其中，RTP表示化極結(jié)果。

本步驟中，通過頻率域化極因子進(jìn)行運(yùn)算，頻率域內(nèi)的運(yùn)算時(shí)間短，效率高，在計(jì)算機(jī)上更易編程實(shí)現(xiàn)，降低了方法的適用門檻。

步驟S6：計(jì)算化極結(jié)果的垂直梯度和總梯度，并計(jì)算其互相關(guān)系數(shù)，以最大互相關(guān)系數(shù)確定剩磁方向。

本步驟中，利用了磁異常的總梯度是垂直梯度的包絡(luò)這一物理性質(zhì)，二者在垂直磁化時(shí)具有對稱性，當(dāng)二者的互相關(guān)系數(shù)最大時(shí)，磁異常最接近垂直磁化狀態(tài)；也就是說，此時(shí)的選用的剩磁方向能夠使得化極結(jié)果最接近垂直磁化，因此該方向就是磁源體的真實(shí)剩磁方向。根據(jù)步驟S5的頻率域化極因子中進(jìn)行變傾角化極，求取化極結(jié)果的垂直梯度和總梯度，并計(jì)算其互相關(guān)系數(shù)，以最大互相關(guān)系數(shù)來剩磁方向。

步驟S7：根據(jù)地磁場方向和剩磁方向，并根據(jù)頻率域化極因子對矩陣T’和I在頻率域內(nèi)進(jìn)行變傾角化極處理，計(jì)算研究區(qū)磁異常的變傾角化極結(jié)果。

其中，在本步驟之前，首先從矩陣I中提取平均地磁傾角為平均地磁方向，然后計(jì)算出研究區(qū)中的實(shí)際地磁方向與平均地磁方向的差

本步驟具體利用泰勒公式進(jìn)行變傾角化極，采用的公式是：

式中，RTP_diff是變傾角化極結(jié)果，RTP_mean是使用研究區(qū)平均地磁傾角和地磁偏角得到的化極結(jié)果。

由于本步驟利用泰勒公式，實(shí)現(xiàn)了變傾角化極，避免了用單一地磁傾角進(jìn)行化極時(shí)的失真問題，從數(shù)學(xué)上解決了緯度跨度較大時(shí)研究區(qū)磁異?；瘶O的問題，因此提高了磁法勘探的準(zhǔn)確度。

本實(shí)施例的方法適用于大、中、小比例尺的磁異常的處理。

本發(fā)明還分別以南充-旺蒼地區(qū)和石柱-大竹地區(qū)的磁異常資料為例，結(jié)合前文所述的方法，給予進(jìn)一步說明。

首先，以南充-旺蒼地區(qū)的磁異常資料為例，按照前文所述的方法，對磁異常資料進(jìn)行處理。首先請參照圖2-4。

首先利用磁力儀在野外采集資料，獲得南充-旺蒼地區(qū)的磁場數(shù)據(jù)資料，然后將獲得的磁場資料導(dǎo)入計(jì)算機(jī)，對數(shù)據(jù)進(jìn)行各種校正，得到南充-旺蒼地區(qū)(南北210km×東西169km)的磁異常資料。然后利用surfer軟件中的kriging插值方法將上述磁異常資料按照1km間距進(jìn)行網(wǎng)格化，形成規(guī)則的磁異常資料，記為矩陣T(1:210,1:169)，詳見圖2。

進(jìn)一步的，利用國際地磁參考場(IGRF)，獲得該地區(qū)的地磁場方向。本地區(qū)磁異常數(shù)據(jù)采集自1971年，9月1日，緯度范圍約為30.53°N～32.26°N，根據(jù)IGRF模型，地磁傾角范圍約為44～46.5，地磁偏角約為2°。地磁傾角按照以下算法，整理成與磁異常矩陣T等大的矩陣I。

①設(shè)置矩陣

②利用矩陣A的偽逆矩陣計(jì)算磁場方向；

③P＝pinv(A)×V，pinv表示偽逆矩陣，P為計(jì)算磁場方向時(shí)的參數(shù)矩陣；

④設(shè)置網(wǎng)格坐標(biāo)陣[x,y]＝meshgrid(1:169,1:210)；

⑤I＝P(1)+P(2)×y+P(3)×x。

進(jìn)一步的，從矩陣I中提取平均地磁傾角計(jì)算出研究區(qū)中的實(shí)際地磁方向與平均地磁方向的差

進(jìn)一步的，以“余弦衰減置零擴(kuò)邊法”給矩陣T(1:210,1:169)鑲邊，將矩陣T(1:210,1:169)轉(zhuǎn)變成行數(shù)和列數(shù)均為2的最小整數(shù)冪的新矩陣T’(-22:233,-43:212)。具體鑲邊公式如前文所述，在此不再贅述。

設(shè)計(jì)的頻率域化極因子，其中頻率域化極因子的公式如前文所述，在此不再贅述。進(jìn)一步對對頻率域化極因子做反傅里葉變換，得到化極結(jié)果。具體反傅里葉變換的公式如前文所述，在此不再贅述。

進(jìn)一步的，在頻率域內(nèi)計(jì)算化極結(jié)果的垂直梯度和總梯度，并求出其互相關(guān)系數(shù)C。C是剩磁方向的函數(shù)。C最大時(shí)對應(yīng)的剩磁方向就是本地區(qū)磁源體的剩磁方向(i,d)。圖3是互相關(guān)系數(shù)的結(jié)果，橫坐標(biāo)是偏角，縱坐標(biāo)是傾角，只顯示了最大互相關(guān)系數(shù)附近的圖像。該地區(qū)的剩磁方向?yàn)槭４艃A角3.3°，剩磁偏角-0.5°。

進(jìn)一步的，如圖4，利用該地區(qū)剩磁方向和地磁場方向，通過矩陣T’、I在頻率域內(nèi)進(jìn)行變傾角化極處理。計(jì)算該地區(qū)磁異常的變傾角化極結(jié)果。該地區(qū)地磁場傾角，即地磁場方向?yàn)?5°左右，而剩磁傾角，即剩磁方向僅為3.3°，二者相差約14倍；對比圖2與圖4可以發(fā)現(xiàn)，磁異常幅值由-240～+420nT變?yōu)?300～+800nT，二者相差較大；正磁異常中心位置南充附近，在考慮剩磁影響后，磁異常峰值位于南部-儀隴之間，北偏達(dá)77.5km，且是唯一的等值線閉合處。原始磁異常(圖2)與化極磁異常(圖4)的巨大差異說明剩余磁化在該地區(qū)是很強(qiáng)的。

以下以石柱-大竹地區(qū)的磁異常資料為例，按照前文所述的方法，對磁異常資料進(jìn)行處理。首先請參照圖5-7。

首先利用磁力儀在野外采集資料，獲得石柱-大竹地區(qū)的磁場數(shù)據(jù)資料，然后將獲得的磁場資料導(dǎo)入計(jì)算機(jī)，對數(shù)據(jù)進(jìn)行各種校正，得到石柱-大竹地區(qū)(南北168km×東西150km)的磁異常資料。然后利用surfer軟件中的kriging插值方法將上述磁異常資料按照1km間距進(jìn)行網(wǎng)格化，形成規(guī)則的磁異常資料，記為矩陣T(1:168,1:150)，詳見圖5。

進(jìn)一步的，利用國際地磁參考場(IGRF)，獲得該地區(qū)的地磁場方向。本地區(qū)磁異常數(shù)據(jù)采集自1971年，9月1日，緯度范圍約為29.69°N～30.95°N，根據(jù)IGRF模型，地磁傾角范圍約為43.5～45，地磁偏角約為2°。地磁傾角按照以下算法，整理成與磁異常矩陣T等大的矩陣I。

①設(shè)置矩陣

②利用矩陣A的偽逆矩陣計(jì)算磁場方向；

③P＝pinv(A)×V，pinv表示偽逆矩陣，P為計(jì)算磁場方向時(shí)的參數(shù)矩陣；

④設(shè)置網(wǎng)格坐標(biāo)陣[x,y]＝meshgrid(1:150,1:168)；

⑤I＝P(1)+P(2)×y+P(3)×x。

進(jìn)一步的，從矩陣I中提取平均地磁傾角計(jì)算出研究區(qū)中的實(shí)際地磁方向與平均地磁方向的差

進(jìn)一步的，以“余弦衰減置零擴(kuò)邊法”給矩陣T(1:168,1:150)鑲邊，將矩陣T(1:168,1:150)轉(zhuǎn)變成行數(shù)和列數(shù)均為2的最小整數(shù)冪的新矩陣T’(-43:212,-52:203)。具體鑲邊公式如前文所述，在此不再贅述。

設(shè)計(jì)的頻率域化極因子，其中頻率域化極因子的公式如前文所述，在此不再贅述。進(jìn)一步對對頻率域化極因子做反傅里葉變換，得到化極結(jié)果。具體反傅里葉變換的公式如前文所述，在此不再贅述。

進(jìn)一步的，在頻率域內(nèi)計(jì)算化極結(jié)果的垂直梯度和總梯度，并求出其互相關(guān)系數(shù)C。C是剩磁方向的函數(shù)。C最大時(shí)對應(yīng)的剩磁方向就是本地區(qū)磁源體的剩磁方向(i,d)。圖6是互相關(guān)系數(shù)的結(jié)果，橫坐標(biāo)是偏角，縱坐標(biāo)是傾角，只顯示了最大互相關(guān)系數(shù)附近的圖像。該地區(qū)的剩磁方向?yàn)槭４艃A角21°，剩磁偏角-41°。

進(jìn)一步的，如圖7，利用該地區(qū)剩磁方向和地磁場方向，通過矩陣T’、I在頻率域內(nèi)進(jìn)行變傾角化極處理。計(jì)算該地區(qū)磁異常的變傾角化極結(jié)果。該地區(qū)地磁場傾角為44°左右，而剩磁傾角為21°，二者相差約2倍；對比圖5與圖7可以發(fā)現(xiàn)，磁異常幅值由-180～+200nT變?yōu)?320～+280nT，二者相差較小；正磁異常中心位置豐都-石柱之間，在考慮剩磁影響后，磁異常峰值位于墊江，北偏約25km。原始磁異常(圖5)與化極磁異常(圖7)差異不大，說明剩余磁化在該地區(qū)較弱。

綜上所述，本發(fā)明能夠利用地磁場方向和剩磁方向?qū)?shí)測磁異常進(jìn)行變傾角化極，提高磁法勘探的準(zhǔn)確度。

以上所述僅為本發(fā)明的實(shí)施例，并非因此限制本發(fā)明的專利范圍，凡是利用本發(fā)明說明書及附圖內(nèi)容所作的等效結(jié)構(gòu)或等效流程變換，或直接或間接運(yùn)用在其他相關(guān)的技術(shù)領(lǐng)域，均同理包括在本發(fā)明的專利保護(hù)范圍內(nèi)。