鐘煬，管彥武，石甲強(qiáng)，肖鋒

(吉林大學(xué) 地球探測科學(xué)與技術(shù)學(xué)院，吉林 長春 130026)

0 引言

國內(nèi)外諸多學(xué)者針對IGRF模型的研究，早期集中在分析模型精度及驗(yàn)證模型適用范圍等問題上。Barraclough等[1]利用大量地磁臺站數(shù)據(jù)推導(dǎo)出地磁場長期變化的球諧模型，從模型中導(dǎo)出了磁極、地磁極和偏心磁極的位置及其變化速率，并比較了該模型與IGRF模型的差別，指出在模型中加入加速變化系數(shù)對地磁場模型的擬合度提高具有顯著的作用。安振昌[2-9]首次在國內(nèi)較為詳細(xì)地介紹了國際地磁參考場、發(fā)展歷史及其主要內(nèi)容，并利用第七代IGRF模型計(jì)算1900～2000年非偶極子磁場的全球變化和1995年中國及鄰區(qū)地磁場。林云芳等[10]利用東南亞及其鄰近地區(qū)的地磁臺的實(shí)測資料與IGRF模式進(jìn)行了對比，證明IGRF模型存在偏差，并指出應(yīng)用時應(yīng)該選取時間間隔短的IRGF模型。Matteo等[11]通過分析1991～2010年間電離層磁場強(qiáng)度的衛(wèi)星數(shù)據(jù)，與IGRF模型計(jì)算的結(jié)果進(jìn)行比較，結(jié)果表明即使在地磁場活躍的條件下，在一定的位置和高度范圍內(nèi)，IGRF在電離層仍然具有很高的精度。

為適用于高精度地質(zhì)構(gòu)造研究工作，徐文耀[12]指出在1945～1955年IGRF模型中的高階高斯系數(shù)(n>7)與其他時期的模型相比，存在較大偏差，并通過自然正交分量分析法對1945～1955年IGRF模型中的高次高斯系數(shù)進(jìn)行修正，使這些模型的高斯系數(shù)具有較好的平滑時變特性。Maus等[13]美國地球物理數(shù)據(jù)中心(NGDC)的科學(xué)家們計(jì)算出16～720階球諧系數(shù)，成功建立巖石圈磁場模型NGDC-720。張素琴等[14]利用中國部分地磁臺站1990～2003 年的年均值資料， 研究了地磁臺站年均值與IGRF模型計(jì)算值的一致程度，結(jié)果顯示雖然存在差異，但是各地磁要素年均值與IGRF模型值的差值的標(biāo)準(zhǔn)差均低于IGRF模型的誤差水平，故一致性較好。

在IGRF模型的應(yīng)用方面， Molina(1987)等[15]針對意大利等歐洲中部國家，建立了基于IGRF模型的意大利地磁參考場(ITGRF)和歐洲中部地磁參考場(EGRF)。任國泰[16]利用IGRF模型結(jié)合地面測量資料和地磁臺資料分析了東亞大陸磁場。安振昌[17]比較和討論了中國地區(qū)地磁場區(qū)域模型與全球模型，得出地磁場區(qū)域模型，更詳細(xì)準(zhǔn)確地表示我國地磁場的分布。Smart等[18-19]運(yùn)用磁場模型與運(yùn)動軌跡學(xué)原理，利用IGRF模型導(dǎo)出飛機(jī)飛行時截止剛度(VCR)，進(jìn)而通過線性擬合確定其飛行軌跡。Davis[20]利用Matlab編程實(shí)現(xiàn)了IGRF模型關(guān)于場源外任意位置三分量地磁場的計(jì)算。

此外，馮春[21]、柴松均[22]及楊夢雨[23]等人，分別運(yùn)用MATLAB、C及Java語言實(shí)現(xiàn)了任意給定點(diǎn)位的地磁七要素計(jì)算。但隨著航空全張量磁梯度測量技術(shù)的發(fā)展，對利用IGRF模型計(jì)算地球主磁場的全張量磁梯度有著迫切的需求，故筆者梳理了地磁場七要素及張量的計(jì)算流程，推導(dǎo)出計(jì)算任意給定點(diǎn)位全張量地磁梯度的公式，繪制出某地區(qū)地磁場全張量等值線圖，并驗(yàn)證了算法的正確性。

1 利用IGRF模型計(jì)算全張量地磁梯度

目前，球諧函數(shù)法被廣泛應(yīng)用于計(jì)算地磁場參考模型中。在該方法中，認(rèn)為地球是均勻球體，在球坐標(biāo)系下地磁場源區(qū)外任意一點(diǎn)坐標(biāo)(r,θ,φ)磁位V的拉普拉斯方程形式如式(1)所示：

(1)

其中：r、θ和φ是地心球坐標(biāo)系構(gòu)成要素，r是從地心起算的徑向距離，單位為km，θ是地心余緯度(θ=90°-L，L是地心緯度)，φ是從格林尼治子午線算起的地心經(jīng)度。

利用分離變量法求解式(1)的Dirichlet問題(常微分方程第一類邊界條件)，獲得地磁場源區(qū)外任意一點(diǎn)坐標(biāo)(r，θ，φ)的磁位V表達(dá)式[24]如式(2)所示：

(2)

(3)

(4)

在物理學(xué)中，磁場強(qiáng)度定義為磁位V的負(fù)梯度，在地心球坐標(biāo)系中，對式(4)分別求r、θ和φ的偏導(dǎo)數(shù)，并取負(fù)得到磁場強(qiáng)度的三個分量Ur、Uθ、Uφ：

(5)

(6)

(7)

為計(jì)算全張量磁梯度，保持原坐標(biāo)系不變，在式(5)～(7)中，分別對Ur、Uθ和Uφ求r、θ、φ的偏導(dǎo)數(shù)，得到磁位的二階偏導(dǎo)數(shù)：

(8)

(9)

(10)

(11)

(12)

(13)

Balmino等對地心球坐標(biāo)系下磁位的一階及二階偏導(dǎo)數(shù)(即式(5)～(13))進(jìn)行坐標(biāo)變換，將地心球坐標(biāo)系轉(zhuǎn)換為局部“北東上”直角坐標(biāo)系，其變換方式如式(14)、(15)所示[25]：

(14)

(15)

其中：Bx、By、Bz分別表示局部“北東上”直角坐標(biāo)系下地磁場三分量，Bxx、Bxy、Bxz、Byy、Byz、Bzz表示局部“北東上”直角坐標(biāo)系下地磁場的全張量的6個分量。

然后將局部“北東上”直角坐標(biāo)系xyz轉(zhuǎn)換到符合右手定則且更加常用的局部“北東下”直角坐標(biāo)系x1y1z1，轉(zhuǎn)換關(guān)系為：

(16)

再將z1指向地心的x1y1z1坐標(biāo)系轉(zhuǎn)換到z2沿旋轉(zhuǎn)橢球面法線向下的局部直角坐標(biāo)系x2y2z2，即相當(dāng)于x1y1z1坐標(biāo)系繞y1旋轉(zhuǎn)角度d。d為地心余緯度θ與地理余緯度θ′的差值。圖1中P點(diǎn)表示測點(diǎn)。

圖1 地理坐標(biāo)系與地心坐標(biāo)系Fig.1 Geographic coordinate system and geocentric coordinate system

由圖1可知：

d=θ-θ′，

(17)

(18)

(19)

(20)

b=a·(1-f)，

(21)

其中：a為橢球長半軸，a=6 378.137 km，b為橢球短半軸，f為地球扁率，f=1/298.257 223 563，h為測點(diǎn)的大地高(橢球高度)。

由式(17)可知：

sinθ=sinθ′cosd+cosθ′sind，

(22)

cosθ=cosθ′cosd-sinθ′sind，

(23)

(24)

在實(shí)際應(yīng)用中，通常測量得到測點(diǎn)的大地坐標(biāo)，通過由式(22)～(24)可求得地心余緯度的三角函數(shù)值，并將其應(yīng)用于前述公式。

再由圖1所示坐標(biāo)系x1y1z1與x2y2z2的旋轉(zhuǎn)關(guān)系可知：

(25)

(26)

(27)

利用式(26)及(27)進(jìn)行坐標(biāo)變換，獲得更普遍的z軸沿旋轉(zhuǎn)橢球面法線向下的局部直角坐標(biāo)系(北東下)的地磁場三分量及全張量公式。

2 伴隨勒讓德多項(xiàng)式及其遞推關(guān)系

利用IGRF模型，為計(jì)算地磁場三分量及全張量磁梯度，還需要計(jì)算滿足正交性的施密特規(guī)格化伴隨勒讓德多項(xiàng)式。

數(shù)學(xué)上，勒讓德函數(shù)是指以下勒讓德方程的解：

(28)

勒讓德微分方程(28)的解，可以寫成標(biāo)準(zhǔn)的冪級數(shù)形式，并且當(dāng)方程滿足|v|<1且n為非負(fù)整數(shù)時，方程的解將隨n值的變化而變化并構(gòu)成一組由正交多項(xiàng)式組成的多項(xiàng)式序列，即勒讓德多項(xiàng)式[26]。勒讓德多項(xiàng)式的微分形式，可用羅德里格斯公式表示：

(29)

其中：0～3階勒讓德多項(xiàng)式表達(dá)式如表1所示。

表1 0～3階勒讓德多項(xiàng)式Table 1 Legendre polynomials of order 0～3

在球坐標(biāo)系下求解拉普拉斯方程，可得到如下常微分方程：

(30)

常微分方程式(32)的解序列稱為伴隨勒讓德多項(xiàng)式，又稱為締合勒讓德多項(xiàng)式、連帶勒讓德多項(xiàng)式或關(guān)聯(lián)勒讓德多項(xiàng)式。伴隨勒讓德多項(xiàng)式也可以通過對勒讓德多項(xiàng)式求m次導(dǎo)數(shù)，經(jīng)式(31)換算得到：

(31)

如式(4)所示，磁位的表達(dá)式是球諧函數(shù)的k階展開，可表示為與展開階數(shù)k有關(guān)的多項(xiàng)式之和。為使這些多項(xiàng)式在求和計(jì)算中的相對重要性更加接近，并方便數(shù)值計(jì)算，需要將伴隨勒讓德多項(xiàng)式與一個隨階次變化的因子相乘，即對伴隨勒讓德多項(xiàng)式規(guī)格化。伴隨勒讓德多項(xiàng)式的規(guī)格化主要有高斯規(guī)格化和施密特規(guī)格化兩種方式[27]。采用高斯規(guī)格化伴隨勒讓德多項(xiàng)式的表達(dá)式如下：

(32)

其中：Pn,m(cosθ)為高斯規(guī)格化伴隨勒讓德多項(xiàng)式，Pn,m(cosθ)為伴隨勒讓德多項(xiàng)式。

利用施密特方程計(jì)算的規(guī)格化勒讓德多項(xiàng)式表達(dá)式為：

(33)

兩種規(guī)格化方法的關(guān)系見式(34)所示：

(34)

(35)

0～3階規(guī)格化伴隨勒讓德多項(xiàng)式如表2所示。

表2 0～3階伴隨勒讓德多項(xiàng)式Table 2 Order 0-3 adjoint Legendre polynomials

為了方便計(jì)算機(jī)進(jìn)行編程計(jì)算，需根據(jù)施密特半規(guī)則化勒讓德多項(xiàng)式(33)求取伴隨勒讓德多項(xiàng)式的遞推關(guān)系。采用向前列遞推法，施密特?cái)M規(guī)格化伴隨勒讓德多項(xiàng)式的遞推關(guān)系為：

(36)

地磁場三分量及全張量磁梯度表達(dá)式(26)、(27)中，還含有施密特?cái)M規(guī)格化伴隨勒讓德多項(xiàng)式的微分形式。其遞推關(guān)系如式(37)所示。

在全張量磁梯度表達(dá)式中的Bxx和Bzz分量，除含有勒讓德多項(xiàng)式的一階微分形式外，還含有二階微分形式，其遞推關(guān)系如式(38)所示。

3 實(shí)際地磁場模擬

實(shí)際應(yīng)用中，需要選取地球主磁場梯度較小的區(qū)域和高度作為航磁或FTMG測量學(xué)習(xí)飛行的工區(qū)。同時，為驗(yàn)證上述算法在計(jì)算地磁場七要素及全張量地磁梯度的可行性，選取4°×4°的計(jì)算工區(qū)(經(jīng)度變化范圍為：103.305 6°～107.305 6°，緯度變化范圍為：27.305 6°～31.305 6°)，網(wǎng)格大小為0.1°×0.1°。選取大地高1 km，選定時間為2019年4月7日進(jìn)行計(jì)算。

(37)

(38)

圖2為所選區(qū)域的總地磁場等值線圖，圖中等值線光滑，場值變化范圍為48 693～51 200 nT。

圖2 總地磁場B等值線Fig.2 Contour map of the total geomagnetic field

圖3a為地磁偏角D等值線圖，其變化范圍-2.871 7°～-1.652 5°，圖3b為地磁傾角I等值線圖，其變化范圍為42.748°～49.133 6°。

圖4從左至右依次為地磁場在x、y、z軸的分量Bx、By、Bz。 由于目前還沒有其他直接獲得任意給定點(diǎn)位地球主磁場全張量磁梯度數(shù)據(jù)的方法，無法逐點(diǎn)對比圖5中所選測區(qū)內(nèi)全張量磁梯度場值。但根據(jù)地磁場三分量等值線圖(圖4)，其沿3個坐標(biāo)軸方向的導(dǎo)數(shù)即為全張量地磁梯度值，可以大致判斷6個張量分量符合地磁場梯度量級，且等值線形態(tài)合理。

圖3 地磁場偏角D(a)及傾角I(b)等值線Fig.3 Contour map of geomagnetic field deviation(a) and dip angle(b)

為進(jìn)一步驗(yàn)證圖5中全張量地磁梯度計(jì)算值的正確性，將Bxx、Byy、Bzz相加如圖6所示。

圖4 地磁場三分量(a)Bx(b)By(c)Bz等值線Fig.4 Three-component contour map of the geomagnetic field

圖5 全張量地磁梯度等值線Fig.5 Full-tensor geomagnetic gradient contour map

如圖6所示，全張量磁梯度的3個分量Bxx、Byy、Bzz之和的變化范圍為：-0.001 1～0.001 1 nT/km，可以認(rèn)為利用該計(jì)算方法所獲得全張量地磁梯度滿足Laplace方程。

為進(jìn)一步驗(yàn)證算法的正確性，選取計(jì)算工區(qū)內(nèi)4個隨機(jī)點(diǎn)位：①四川省成都市(30.67N，104.07E)，②自貢市(29.35N，104.78E)，③瀘州市(28.87N，105.43E)，④德陽市(31.13N，104.38E))，利用美國國家海洋和大氣管理局(NOAA)官方網(wǎng)站(www.ngdc.noaa.gov)內(nèi)全球地磁場數(shù)據(jù)，對比兩種方式計(jì)算獲得地磁場七要素誤差見表3。三分量矢量場在保留小數(shù)點(diǎn)一位的情況下誤差很小，可忽略不計(jì)，從而驗(yàn)證了算法的正確性。

表3 測區(qū)內(nèi)點(diǎn)位磁場值Table 3 Table of point position magnetic field values in the test area

圖6 主對角線上3個張量分量之和等值線Fig.6 Contour map of the sum of the three tensor components on the main diagonal

4 結(jié)論

1)筆者梳理了利用IGRF模型計(jì)算地磁場的原理與步驟，利用4個點(diǎn)位的地磁七要素計(jì)算結(jié)果與美國國家海洋和大氣管理局計(jì)算數(shù)據(jù)進(jìn)行對比，結(jié)果準(zhǔn)確可靠。

2)依據(jù)IGRF模型，推導(dǎo)出全張量地磁梯度的球諧展開形式，給出了施密特規(guī)格化伴隨勒讓德多項(xiàng)式二階導(dǎo)數(shù)的遞推公式。實(shí)驗(yàn)工區(qū)全張量地磁梯度的計(jì)算結(jié)果形態(tài)與量級合理，且滿足Laplace方程，為航空全張量磁梯度測量中選擇學(xué)習(xí)飛行工區(qū)和飛行高度提供了重要的依據(jù)。