一種磁異常化極方法與流程
2023-12-09 17:46:57 3
本發明涉及地球物理勘測技術領域，具體為磁法勘探領域，特別是涉及一種磁異常化極方法。
背景技術：
在磁法勘探中，磁測數據的處理和解釋(如化極、導數換算、分量換算、正演和反演等)通常需要已知磁場方向信息。磁場是一個矢量場，而利用磁力儀測得的磁場強度是標量，沒有磁場方向。磁場方向是感應磁化方向和剩餘磁化方向的矢量和。利用國際地磁參考場獲得感應磁場方向，即地磁場方向；而剩餘磁化方向無法直接獲取。若磁源地質體剩磁不明顯，那麼磁場方向與感應磁場方向基本一致；否則，較強的剩磁會使磁場方向與感磁方向出現很大偏差，使得地下磁性體產生斜磁化現象。此時，若直接利用地磁場方向作為磁場方向進行磁異常化極，以及在此基礎上進行的一系列磁異常數據處理與轉換、正演和反演，其結果都無法反映真實的地質情況。
目前，除採集巖石標本並測定剩磁參數外，國內外學者提出了一些估計剩磁方向的方法。唐俊德(1986)提出利用磁異常三分量確定磁化方向的圓曲線法；Roest和Pilkington(1993)提出基於磁場總梯度模和磁源重力異常水平梯度互相關的估計方法；Medeiros和Silva(1995)提出利用等效源磁矩反演估計磁化方向；甘西(2001)利用磁場的垂直分量確定磁性體磁化方向；Phillips(2005)給出估計磁化方向的Helbig(1963)積分法的直接算法和間接算法；Dannemiller和Li(2006)提出基於化極磁異常垂直梯度與總梯度模互相關的估計方法；Gerovska等(2009)提出基於化極磁異常與總模量磁異常互相關的估計方法。
但是，磁異常化極的準確性還與研究區的緯度變化有關。若研究區的緯度範圍較大，地磁傾角變化也會較大，斜磁化的程度會隨緯度不同而變化。若採用單一的磁傾角作化極處理顯然是不合理的，其結果將無法反映真實的地質情況，降低磁法勘探的準確度。
技術實現要素：
本發明主要解決的技術問題是提供一種磁異常化極方法，能夠利用地磁場方向和剩磁方向對實測磁異常進行變傾角化極，提高磁法勘探的準確度。
為解決上述技術問題，本發明採用的一個技術方案是：提供一種磁異常化極方法，包括以下步驟：採集研究區的磁場數據，對採集到的磁場數據進行校正，得到磁異常數據；將磁異常數據網格化，形成規則的磁異常數據t(x,y)，其中x,y分別表示磁異常數據的平面網格點坐標，記為矩陣T(1:m,1:n)，大小為m×n，m表示行數，n表示點數；獲得研究區的地磁場方向，並將地磁場方向整理成與矩陣T(1:m,1:n)等大的矩陣，記為I(1:m,1:n)；對矩陣T(1:m,1:n)進行鑲邊處理，將矩陣T(1:m,1:n)轉變成行數和列數均為2的最小整數冪的新矩陣T』(1-m』:m+m』,1-n』:n+n』)；設計頻率域化極因子，並對頻率域化極因子做反傅立葉變換，得到化極結果；計算化極結果的垂直梯度和總梯度，並計算其互相關係數，以最大互相關係數確定剩磁方向；根據地磁場方向和剩磁方向，並根據頻率域化極因子對矩陣T』和I在頻率域內進行變傾角化極處理，計算研究區磁異常的變傾角化極結果。
其中，將磁異常數據網格化的步驟包括：利用kriging插值方法將磁異常數據網格化。
其中，將地磁場方向整理成與矩陣T(1:m,1:n)等大的矩陣，記為I(1:m,1:n)的步驟包括：
獲取地磁傾角範圍為a-b；
設置矩陣
利用矩陣A的偽逆矩陣計算磁場方向；
P＝pinv(A)×V，pinv表示偽逆矩陣P為計算磁場方向時的參數矩陣，其使用方法見下式；
設置網格坐標陣[x,y]＝meshgrid(1:n,1:m)；
計算矩陣I＝P(1)+P(2)×y+P(3)×x。
其中，對矩陣T(1:m,1:n)進行鑲邊處理的步驟包括：
採用以下公式進行鑲邊：
其中，T位於T』的中心位置，即T(i,j)＝T』(i,j)，(1≤i≤m,1≤j≤n)。
其中，頻率域化極因子為：
式中，ωx和ωy是x、y方向的圓波數，是地磁方向的單位向量，L、M、N是地磁場的方向餘弦：cosIcosD,cosIsinD,sinI，I和D分別表示地磁場傾角和地磁場偏角；是剩磁方向的單位向量，l,m,n分別是剩磁方向的方向餘弦：cosicosd,cosisind,sini，i和d分別表示剩磁傾角和剩磁偏角。
其中，對頻率域化極因子做反傅立葉變換的步驟包括：
採用以下的公式做反傅立葉變換：
其中，RTP表示化極結果。
其中，根據頻率域化極因子對矩陣T』、I在頻率域內進行變傾角化極處理的步驟之前包括：
從矩陣I中提取平均地磁傾角為平均地磁方向；
計算出研究區中的實際地磁方向與平均地磁方向的差
其中，利用泰勒公式進行變傾角化極，採用的公式是：
式中，RTPdiff是變傾角化極結果，RTPmean是使用研究區平均地磁傾角和地磁偏角得到的化極結果。
本發明的有益效果是：區別於現有技術的情況，本發明提供一種磁異常化極方法，包括以下步驟：採集研究區的磁場數據，對採集到的磁場數據進行校正，得到磁異常數據；將磁異常數據網格化，形成規則的磁異常數據t(x,y)，其中x,y分別表示磁異常數據的平面網格點坐標，記為矩陣T(1:m,1:n)，大小為m×n，m表示行數，n表示點數；獲得研究區的地磁場方向，並將地磁場方向整理成與矩陣T(1:m,1:n)等大的矩陣，記為I(1:m,1:n)；對矩陣T(1:m,1:n)進行鑲邊處理，將矩陣T(1:m,1:n)轉變成行數和列數均為2的最小整數冪的新矩陣T』(1-m』:m+m』,1-n』:n+n』)；設計頻率域化極因子，並對頻率域化極因子做反傅立葉變換，得到化極結果；計算化極結果的垂直梯度和總梯度，並計算其互相關係數，以最大互相關係數確定剩磁方向；根據地磁場方向和剩磁方向，並根據頻率域化極因子對矩陣T』和I在頻率域內進行變傾角化極處理，計算研究區磁異常的變傾角化極結果。因此，本發明能夠利用地磁場方向和剩磁方向對實測磁異常進行變傾角化極，提高磁法勘探的準確度。
附圖說明
圖1是本發明實施例提供的一種磁異常化極方法的流程圖；
圖2為南充-旺蒼地區的磁異常資料；
圖3為南充-旺蒼地區的磁異常垂直梯度與總梯度的互相關係數；
圖4為南充-旺蒼地區考慮剩磁的變傾角化極磁異常；
圖5為石柱-大竹地區的磁異常資料；
圖6為石柱-大竹地區的磁異常垂直梯度與總梯度的互相關係數；
圖7為石柱-大竹地區考慮剩磁的變傾角化極磁異常。
具體實施方式
請參閱圖1，圖1是本發明實施例提供的一種磁異常化極方法的流程圖。如圖1所示，本實施例的方法包括以下步驟：
步驟S1：採集研究區的磁場數據，對採集到的磁場數據進行校正，得到磁異常數據。
本步驟中，具體是通過磁力儀在野外採集研究區的磁場數據，並將獲得的磁場數據導入計算機中，以進行各種校正，以得到異常數據。
步驟S2：將磁異常數據網格化，形成規則的磁異常數據t(x,y)，其中x,y分別表示磁異常數據的平面網格點坐標，記為矩陣T(1:m,1:n)，大小為m×n，m表示行數，n表示點數。
本步驟中，具體是通過surfer軟體中的kriging插值方法將磁異常數據網格化。
步驟S3：獲得研究區的地磁場方向，並將地磁場方向整理成與矩陣T(1:m,1:n)等大的矩陣，記為I(1:m,1:n)。
本步驟中，具體是利用國際地磁參考場，獲得研究區的地磁場方向。本發明選用的國際地磁參考場(international geomagnetic reference field,IGRF)優選是第11代國際地磁參考場(IGRF-11)，由國際地磁學與高空物理學聯合會(IAGA)於2009年12月發布。
其中，本步驟的將地磁場方向整理成與矩陣T(1:m,1:n)等大的矩陣，記為I(1:m,1:n)的步驟具體包括如下：首先獲取地磁傾角範圍為a-b，然後設置矩陣進而利用矩陣A的偽逆矩陣計算磁場方向；
P＝pinv(A)×V，pinv表示偽逆矩陣P為計算磁場方向時的參數矩陣，進一步設置網格坐標陣[x,y]＝meshgrid(1:n,1:m)，最後計算矩陣I＝P(1)+P(2)×y+P(3)×x。
步驟S4：對矩陣T(1:m,1:n)進行鑲邊處理，將矩陣T(1:m,1:n)轉變成行數和列數均為2的最小整數冪的新矩陣T』(1-m』:m+m』,1-n』:n+n』)。
本步驟中，優選以「餘弦衰減置零擴邊法」給矩陣T(1:m,1:n)鑲邊。具體採用以下公式進行鑲邊：
其中，T位於T』的中心位置，即T(i,j)＝T』(i,j)，(1≤i≤m,1≤j≤n)。
本步驟中，利用「餘弦衰減置零擴邊法」給數據進行鑲邊，很好的壓制了邊界效應，提高了數據處理的精度。
步驟S5：設計頻率域化極因子，並對頻率域化極因子做反傅立葉變換，得到化極結果。
本步驟中，頻率域化極因子具體為：
式中，ωx和ωy是x、y方向的圓波數，是地磁方向的單位向量，L、M、N是地磁場的方向餘弦：cosIcosD,cosIsinD,sinI，I和D分別表示地磁場傾角和地磁場偏角；是剩磁方向的單位向量，l,m,n分別是剩磁方向的方向餘弦：cosicosd,cosisind,sini，i和d分別表示剩磁傾角和剩磁偏角。
其中，具體是採用以下的公式做反傅立葉變換：
其中，RTP表示化極結果。
本步驟中，通過頻率域化極因子進行運算，頻率域內的運算時間短，效率高，在計算機上更易編程實現，降低了方法的適用門檻。
步驟S6：計算化極結果的垂直梯度和總梯度，並計算其互相關係數，以最大互相關係數確定剩磁方向。
本步驟中，利用了磁異常的總梯度是垂直梯度的包絡這一物理性質，二者在垂直磁化時具有對稱性，當二者的互相關係數最大時，磁異常最接近垂直磁化狀態；也就是說，此時的選用的剩磁方向能夠使得化極結果最接近垂直磁化，因此該方向就是磁源體的真實剩磁方向。根據步驟S5的頻率域化極因子中進行變傾角化極，求取化極結果的垂直梯度和總梯度，並計算其互相關係數，以最大互相關係數來剩磁方向。
步驟S7：根據地磁場方向和剩磁方向，並根據頻率域化極因子對矩陣T』和I在頻率域內進行變傾角化極處理，計算研究區磁異常的變傾角化極結果。
其中，在本步驟之前，首先從矩陣I中提取平均地磁傾角為平均地磁方向，然後計算出研究區中的實際地磁方向與平均地磁方向的差
本步驟具體利用泰勒公式進行變傾角化極，採用的公式是：
式中，RTPdiff是變傾角化極結果，RTPmean是使用研究區平均地磁傾角和地磁偏角得到的化極結果。
由於本步驟利用泰勒公式，實現了變傾角化極，避免了用單一地磁傾角進行化極時的失真問題，從數學上解決了緯度跨度較大時研究區磁異常化極的問題，因此提高了磁法勘探的準確度。
本實施例的方法適用於大、中、小比例尺的磁異常的處理。
本發明還分別以南充-旺蒼地區和石柱-大竹地區的磁異常資料為例，結合前文所述的方法，給予進一步說明。
首先，以南充-旺蒼地區的磁異常資料為例，按照前文所述的方法，對磁異常資料進行處理。首先請參照圖2-4。
首先利用磁力儀在野外採集資料，獲得南充-旺蒼地區的磁場數據資料，然後將獲得的磁場資料導入計算機，對數據進行各種校正，得到南充-旺蒼地區(南北210km×東西169km)的磁異常資料。然後利用surfer軟體中的kriging插值方法將上述磁異常資料按照1km間距進行網格化，形成規則的磁異常資料，記為矩陣T(1:210,1:169)，詳見圖2。
進一步的，利用國際地磁參考場(IGRF)，獲得該地區的地磁場方向。本地區磁異常數據採集自1971年，9月1日，緯度範圍約為30.53°N～32.26°N，根據IGRF模型，地磁傾角範圍約為44～46.5，地磁偏角約為2°。地磁傾角按照以下算法，整理成與磁異常矩陣T等大的矩陣I。
①設置矩陣
②利用矩陣A的偽逆矩陣計算磁場方向；
③P＝pinv(A)×V，pinv表示偽逆矩陣，P為計算磁場方向時的參數矩陣；
④設置網格坐標陣[x,y]＝meshgrid(1:169,1:210)；
⑤I＝P(1)+P(2)×y+P(3)×x。
進一步的，從矩陣I中提取平均地磁傾角計算出研究區中的實際地磁方向與平均地磁方向的差
進一步的，以「餘弦衰減置零擴邊法」給矩陣T(1:210,1:169)鑲邊，將矩陣T(1:210,1:169)轉變成行數和列數均為2的最小整數冪的新矩陣T』(-22:233,-43:212)。具體鑲邊公式如前文所述，在此不再贅述。
設計的頻率域化極因子，其中頻率域化極因子的公式如前文所述，在此不再贅述。進一步對對頻率域化極因子做反傅立葉變換，得到化極結果。具體反傅立葉變換的公式如前文所述，在此不再贅述。
進一步的，在頻率域內計算化極結果的垂直梯度和總梯度，並求出其互相關係數C。C是剩磁方向的函數。C最大時對應的剩磁方向就是本地區磁源體的剩磁方向(i,d)。圖3是互相關係數的結果，橫坐標是偏角，縱坐標是傾角，只顯示了最大互相關係數附近的圖像。該地區的剩磁方向為剩磁傾角3.3°，剩磁偏角-0.5°。
進一步的，如圖4，利用該地區剩磁方向和地磁場方向，通過矩陣T』、I在頻率域內進行變傾角化極處理。計算該地區磁異常的變傾角化極結果。該地區地磁場傾角，即地磁場方向為45°左右，而剩磁傾角，即剩磁方向僅為3.3°，二者相差約14倍；對比圖2與圖4可以發現，磁異常幅值由-240～+420nT變為-300～+800nT，二者相差較大；正磁異常中心位置南充附近，在考慮剩磁影響後，磁異常峰值位於南部-儀隴之間，北偏達77.5km，且是唯一的等值線閉合處。原始磁異常(圖2)與化極磁異常(圖4)的巨大差異說明剩餘磁化在該地區是很強的。
以下以石柱-大竹地區的磁異常資料為例，按照前文所述的方法，對磁異常資料進行處理。首先請參照圖5-7。
首先利用磁力儀在野外採集資料，獲得石柱-大竹地區的磁場數據資料，然後將獲得的磁場資料導入計算機，對數據進行各種校正，得到石柱-大竹地區(南北168km×東西150km)的磁異常資料。然後利用surfer軟體中的kriging插值方法將上述磁異常資料按照1km間距進行網格化，形成規則的磁異常資料，記為矩陣T(1:168,1:150)，詳見圖5。
進一步的，利用國際地磁參考場(IGRF)，獲得該地區的地磁場方向。本地區磁異常數據採集自1971年，9月1日，緯度範圍約為29.69°N～30.95°N，根據IGRF模型，地磁傾角範圍約為43.5～45，地磁偏角約為2°。地磁傾角按照以下算法，整理成與磁異常矩陣T等大的矩陣I。
①設置矩陣
②利用矩陣A的偽逆矩陣計算磁場方向；
③P＝pinv(A)×V，pinv表示偽逆矩陣，P為計算磁場方向時的參數矩陣；
④設置網格坐標陣[x,y]＝meshgrid(1:150,1:168)；
⑤I＝P(1)+P(2)×y+P(3)×x。
進一步的，從矩陣I中提取平均地磁傾角計算出研究區中的實際地磁方向與平均地磁方向的差
進一步的，以「餘弦衰減置零擴邊法」給矩陣T(1:168,1:150)鑲邊，將矩陣T(1:168,1:150)轉變成行數和列數均為2的最小整數冪的新矩陣T』(-43:212,-52:203)。具體鑲邊公式如前文所述，在此不再贅述。
設計的頻率域化極因子，其中頻率域化極因子的公式如前文所述，在此不再贅述。進一步對對頻率域化極因子做反傅立葉變換，得到化極結果。具體反傅立葉變換的公式如前文所述，在此不再贅述。
進一步的，在頻率域內計算化極結果的垂直梯度和總梯度，並求出其互相關係數C。C是剩磁方向的函數。C最大時對應的剩磁方向就是本地區磁源體的剩磁方向(i,d)。圖6是互相關係數的結果，橫坐標是偏角，縱坐標是傾角，只顯示了最大互相關係數附近的圖像。該地區的剩磁方向為剩磁傾角21°，剩磁偏角-41°。
進一步的，如圖7，利用該地區剩磁方向和地磁場方向，通過矩陣T』、I在頻率域內進行變傾角化極處理。計算該地區磁異常的變傾角化極結果。該地區地磁場傾角為44°左右，而剩磁傾角為21°，二者相差約2倍；對比圖5與圖7可以發現，磁異常幅值由-180～+200nT變為-320～+280nT，二者相差較小；正磁異常中心位置豐都-石柱之間，在考慮剩磁影響後，磁異常峰值位於墊江，北偏約25km。原始磁異常(圖5)與化極磁異常(圖7)差異不大，說明剩餘磁化在該地區較弱。
綜上所述，本發明能夠利用地磁場方向和剩磁方向對實測磁異常進行變傾角化極，提高磁法勘探的準確度。
以上所述僅為本發明的實施例，並非因此限制本發明的專利範圍，凡是利用本發明說明書及附圖內容所作的等效結構或等效流程變換，或直接或間接運用在其他相關的技術領域，均同理包括在本發明的專利保護範圍內。