高精度磁測梯度改正的 Matlab程序實現(xiàn)

陳林 李暖 陳偉剛

(青海有色地質礦產(chǎn)勘查局地質礦產(chǎn)勘查院，青海西寧 810007)

地球磁場磁化作用于巖(礦)石所產(chǎn)生的磁性，是研究磁異常的理論基礎，通過研究現(xiàn)代樣品的原生剩磁推測地質時期的磁性特征，使得研究地質構造的演化提供了可能(羅效寬等，1991)［1］。當進行大面積高精度磁測工作時，需要進行正常梯度改正。此時要用查全國地磁圖的辦法將不能滿足正常場梯度改正的精度要求，而球諧分析方法則可滿足這一要求。目前，由于球諧分析方法公式復雜，參數(shù)較多，需要花大量的人力物力才能進行正常場的梯度改正，同時，由于市面上出現(xiàn)的同類軟件價格昂貴，且購買一個軟件只能安裝在一臺電腦上。因此，作者利用Matlab軟件編程的簡易性、函數(shù)的多樣性、數(shù)學計算的快速性等特點，編寫了一套高精度磁測梯度改正的簡易程序，解決了這一難題。

1 地球磁場球諧分析的原理

1838年由高斯首先提出球諧分析的方法，該方法是表示全球范圍地磁場的分布及其長期變化的數(shù)學方法。假設地球是均勻磁化球體，球體半徑為R，N為地理北極。這里設地球旋轉軸與地磁軸重合，故N也表示地磁北極。若采用球坐標系，如圖1所示。

圖1 球極坐標系示意圖

坐標原點為球心，球外任一點的地心距為r，余緯度為θ(θ= 90°－ψ，ψ為緯度)，經(jīng)度為λ。則在地磁場源區(qū)之外空間域坐標系(r，θ，ψ)中，磁位u的拉普拉斯方程可以寫成:

對上式采用分離變量法，即令U(r，θ，λ)=R(r)·H(θ)·Φ (λ)，則可解得拉普拉斯方程的一般解，從而可分別獲得其內源場和外源場的磁位球諧表達式。若設外源場磁位為零，則內源場的磁位球諧一般表達式為:

其中，Pn(cosθ)為勒讓德多項式。

由此，便可得到相應三個軸向磁場強度的三分量的球諧表達式。地磁場感應強度的三個分量即北向水平分量X、東向水平分量Y、垂直分量Z(注意這里定義X軸指北為正，Z軸向下為正):

2 Matlab編程及應用

2.1 技術路線

由上述球諧分析方法的原理，需要進行n階(n+1)次的迭代計算，因此，地磁場磁感應強度的三個分量X，Y，Z值在Matlab軟件中將分別進行兩次套合循環(huán)計算。程序運行路線見圖2。

2.2 Matlab軟件編程源代碼

根據(jù)上述編程技術路線，其程序源代碼如下:

圖2 Matlab軟件程序運行路線圖

其中，需要注意的是:

2.3 實際應用

為驗證程序計算結果，實際選取4個點進行計算，計算結果見表1。

表1 實際計算結果表

3 結語

1)本程序經(jīng)過在高精度磁測實際操作中，快速解決了正常場梯度改正的問題，并且有效精度能達到4位小數(shù)，從而快速解決了利用人力所需大量時間的問題。2)Matlab軟件由于含有大量的數(shù)學計算函數(shù)，即便是勒讓德多項式這一復雜度計算，在編程過程中，也只需一個函數(shù)代碼即能解決。3)由于本程序只是為解決人力大量計算的輔助簡易程序，并且Matlab軟件是以C語言及數(shù)理方程為基礎，其編程有相應的數(shù)學及語言格式，人機窗口操作上則為不便，因此，使得程序存在變量較多，操作繁瑣這一局限性。

［1］ 羅孝寬，郭紹雍.應用地球物理教重力：磁法篇［M］.北京：地質出版社，1991.