地磁場(chǎng)中磁目標(biāo)的等效衰減磁矩定位方法

賈文抖， 林春生， 孫玉繪， 翟國(guó)君

(1.海軍工程大學(xué) 兵器工程學(xué)院， 湖北 武漢 430033； 2.海軍工程大學(xué) 導(dǎo)航工程系， 湖北 武漢 430033；3.海軍海洋測(cè)繪研究所， 天津 300061)

0 引言

磁梯度張量具有能夠有效克服地磁場(chǎng)的干擾、突出目標(biāo)磁信號(hào)、更加全面地反映磁場(chǎng)細(xì)節(jié)等突出的優(yōu)勢(shì)，因此近年來(lái)基于磁梯度張量的探測(cè)技術(shù)逐步成為磁探測(cè)領(lǐng)域的研究熱點(diǎn)，德國(guó)、美國(guó)、澳大利亞等國(guó)已相繼研制出磁梯度張量探測(cè)系統(tǒng)，進(jìn)行了大量的野外試驗(yàn)[1]，并逐步應(yīng)用于軍事和民用領(lǐng)域。

基于歐拉齊次方程的磁梯度張量定位方法是早期人們研究的一個(gè)重點(diǎn)，但該方法需要獲取磁目標(biāo)磁場(chǎng)中測(cè)量點(diǎn)上干凈的目標(biāo)三分量磁信號(hào)，這在地磁場(chǎng)環(huán)境中是難以實(shí)現(xiàn)的，實(shí)際應(yīng)用效果也不夠理想。磁梯度張量不變量是利用磁梯度張量數(shù)據(jù)計(jì)算出來(lái)的一些不隨測(cè)量坐標(biāo)系變化而變化的量，它反映了磁目標(biāo)的磁場(chǎng)特征，不易受地磁場(chǎng)的干擾，因此適用于地磁場(chǎng)中磁目標(biāo)的探測(cè)。另外，由于磁梯度張量不變量還具有不隨測(cè)量坐標(biāo)系的變化而改變的性質(zhì)，適合于搭載在移動(dòng)探測(cè)平臺(tái)上。常見(jiàn)的不變量有磁梯度張量的特征值、跡、等效衰減磁矩μ等[2-3]。

江勝華等[4]依據(jù)磁梯度張量縮并理論對(duì)磁梯度張量的模Gt及參數(shù)k的空間分布特征進(jìn)行了分析，結(jié)果顯示Gt和參數(shù)k呈橢球形分布，并擬合得到k值分布的近似計(jì)算公式。陳謹(jǐn)飛等[5]建立了正六面體測(cè)量模型，利用磁梯度張量的模Gt進(jìn)行磁目標(biāo)定位，但由于參數(shù)k的存在使得定位結(jié)果有一定的偏差。若能對(duì)參數(shù)k的影響進(jìn)行補(bǔ)償，或者將其從磁梯度張量的模中完全消除掉，則可以進(jìn)一步提高磁定位的準(zhǔn)確性。Nara等[6]研究發(fā)現(xiàn)影響參數(shù)k分布變化的因素是磁矩與位置矢量之間的夾角，并給出了參數(shù)k隨夾角變化的計(jì)算公式。呂俊偉等[2]研究發(fā)現(xiàn)參數(shù)k可以用磁梯度張量特征值的組合公式進(jìn)行計(jì)算，從磁梯度張量的模Gt中消除參數(shù)k后，得到了一個(gè)僅與距離有關(guān)的不變量。該不變量就是等效衰減磁矩μ，而磁梯度張量的模Gt與等效衰減磁矩μ之間的比值就是參數(shù)k. 此外，Sui等[7]也對(duì)如何消除參數(shù)k帶來(lái)的橢圓誤差進(jìn)行了研究，并進(jìn)行了相應(yīng)的實(shí)驗(yàn)，修正效果較好。實(shí)際上，關(guān)于不變量——等效衰減磁矩μ的研究在國(guó)外文獻(xiàn)中早有出現(xiàn)，在1985年Wilson[8]就提出過(guò)與磁目標(biāo)磁矩方向無(wú)關(guān)的不變量。Beiki等[9]研究發(fā)現(xiàn)當(dāng)直線上的測(cè)點(diǎn)距離磁目標(biāo)最近時(shí)，不變量等效衰減磁矩μ取得極大值這一特性。Clark[10]計(jì)算了磁梯度張量的3個(gè)特征值，發(fā)現(xiàn)可以用磁梯度張量的特征值直接計(jì)算出等效衰減磁矩μ，繼而可以根據(jù)μ的分布規(guī)律確定磁目標(biāo)的位置參數(shù)。Beiki等[11]進(jìn)一步研究發(fā)現(xiàn)等效衰減磁矩μ同樣滿足歐拉齊次方程，并可以利用多點(diǎn)處μ的梯度值反演得到磁目標(biāo)的位置。

在計(jì)算文獻(xiàn)[2，4]中正六面體模型的六面體側(cè)面磁梯度張量數(shù)據(jù)時(shí)，首先利用口字型布置的磁力儀測(cè)得的磁場(chǎng)數(shù)據(jù)近似計(jì)算磁梯度張量中的9個(gè)量，然后計(jì)算出相應(yīng)的磁梯度張量不變量；計(jì)算目標(biāo)源與測(cè)點(diǎn)的距離r時(shí)，根據(jù)上下2個(gè)側(cè)面的磁梯度張量不變量與距離之間的比例關(guān)系，借助麥克勞林公式簡(jiǎn)化近似計(jì)算得到距離r的計(jì)算公式，本文稱(chēng)這種方法為方法1. 為了進(jìn)一步降低上述計(jì)算誤差，得到更精確的磁梯度張量信息和距離r，本文對(duì)方法1進(jìn)行改進(jìn)，在計(jì)算正六面?zhèn)让嬷行狞c(diǎn)的磁梯度張量中的分量時(shí)，采用坐標(biāo)系旋轉(zhuǎn)的方法將口字型布置的磁力儀陣列變成十字型測(cè)量陣列，相比方法1，所用到的磁力儀數(shù)目由4減少到2. 由于正六面體中心點(diǎn)既不位于正六面體中某個(gè)正方形中心點(diǎn)也不位于十字型交叉點(diǎn)，該中心點(diǎn)的磁梯度張量既不能利用口字型磁力儀陣列也不能利用十字型磁力儀陣列差分計(jì)算得到。為此根據(jù)全微分概念，提出一種微分求解方法，利用正六面體上下對(duì)角磁力儀測(cè)量的磁梯度數(shù)據(jù)建立微分方程組，求解計(jì)算磁梯度張量中的9個(gè)元素量，解決了正六面體中心點(diǎn)磁梯度張量難以計(jì)算的問(wèn)題。在計(jì)算目標(biāo)源與測(cè)點(diǎn)的距離r時(shí)，直接利用等效衰減磁矩μ的梯度公式推導(dǎo)得到r，消除了方法1中借助麥克勞林公式簡(jiǎn)化計(jì)算帶來(lái)的誤差。最后根據(jù)μ的圓周分布特征進(jìn)行磁目標(biāo)位置參數(shù)的反演，實(shí)現(xiàn)磁目標(biāo)的定位。

1 基于等效衰減磁矩的定位方法

1.1 磁梯度張量特征值的旋轉(zhuǎn)不變特性

將磁偶極子視為研究對(duì)象，以磁偶極子位置為原點(diǎn)(0, 0, 0)建立直角坐標(biāo)系Oxyz，磁偶極子的磁矩M=(mx,my,mz)，位置矢量為r=(x,y,z)的測(cè)量點(diǎn)Q處磁感應(yīng)強(qiáng)度為

(1)

B的三分量形式為

(2)

磁梯度張量G為磁場(chǎng)三分量在空間x、y、z方向上的導(dǎo)數(shù)，其形式為

(3)

式中：

(4)

下標(biāo)i、j表示x、y、z方向中的任意2個(gè)，δij是克羅內(nèi)克函數(shù)，當(dāng)i=j時(shí)δij=1，當(dāng)i≠j時(shí)δij=0，ri、rj表示距離矢量r=(x,y,z)中的任意兩個(gè)分量。

根據(jù)特征方程det (G-λI)=0(I是與G維數(shù)相同的單位對(duì)角陣，λ是待求的磁梯度張量特征值)，計(jì)算得到G的特征值為

(5)

式中：

假設(shè)坐標(biāo)系Oxyz繞x軸旋轉(zhuǎn)角后得到坐標(biāo)系Oxαyαzα. 在坐標(biāo)系Oxαyαzα中，測(cè)量點(diǎn)Q的坐標(biāo)、磁矩M可分別表示為

(6)

(7)

將(6)式、(7)式代入(5)式，計(jì)算得到坐標(biāo)系Oxαyαzα中測(cè)量點(diǎn)Q處磁梯度張量對(duì)應(yīng)的特征值為

(8)

由(8)式可知，磁偶極子上坐標(biāo)系Oxyz繞x軸旋轉(zhuǎn)前后磁梯度張量G的特征值不變。當(dāng)坐標(biāo)系繞y軸、z軸旋轉(zhuǎn)時(shí)，同樣可得磁梯度張量特征值不隨測(cè)量坐標(biāo)系的旋轉(zhuǎn)而變化。由于位于磁偶極子處的任意直角坐標(biāo)系都可以由坐標(biāo)系Oxyz繞經(jīng)x軸、y軸、z軸旋轉(zhuǎn)得到，根據(jù)以上分析可知，當(dāng)固定了測(cè)量點(diǎn)與磁目標(biāo)點(diǎn)的相對(duì)位置，同一測(cè)量點(diǎn)的磁梯度張量特征值在不同直角坐標(biāo)系中都是相等的。因此，磁梯度張量的特征值是磁梯度張量的旋轉(zhuǎn)不變量。

1.2 等效衰減磁矩的定位原理

(9)

式中：

(10)

根據(jù)特征方程

det (G-λI)=-(μcosβ-λ)[(2μcosβ+
λ)(μcosβ-λ)+μ2sin2β]=0,

求得特征值為

(11)

由(11)式可得

(12)

當(dāng)矢量r與磁矩M同向(β=0)或反向(β=180)時(shí)G退化為對(duì)角矩陣，對(duì)角線上的元素即為特征值，磁梯度張量的特征值與等效衰減磁矩之間的關(guān)系同樣滿足(12)式。由于特征值是磁梯度張量的不變量，因此μ也是磁梯度張量的不變量，且與位置矢量r的模的4次方呈反比。

(13)

(14)

根據(jù)(13)式可得，測(cè)量點(diǎn)Q與磁偶極子之間的距離為

(15)

則位置矢量r為

r=-rμe.

(16)

2 正六面體定位模型

(17)

以計(jì)算x軸正向?qū)?yīng)面中心的μx-為例，磁力儀S1、S2、S3和S4測(cè)量的磁感應(yīng)強(qiáng)度分別為

(18)

由第1節(jié)的分析可知，μ不隨計(jì)算坐標(biāo)系的變化而變化。將磁力儀S1、S2、S3和S4測(cè)得的3個(gè)分量數(shù)據(jù)對(duì)應(yīng)到坐標(biāo)系Ox45°y45°z45°中(見(jiàn)圖3)，坐標(biāo)系Ox45°y45°z45°由坐標(biāo)系Oxyz繞x軸順時(shí)針旋轉(zhuǎn)45°得到。

測(cè)量數(shù)據(jù)由坐標(biāo)系Oxyz對(duì)應(yīng)到坐標(biāo)系Ox45°y45°z45°中的轉(zhuǎn)換公式為

(19)

式中：Bx、By、Bz為坐標(biāo)系Oxyz中磁力儀直接測(cè)得磁場(chǎng)的3個(gè)分量；Bx45°、By45°、Bz45°為對(duì)應(yīng)到坐標(biāo)系Ox45°y45°z45°中的3個(gè)分量。

在坐標(biāo)系Ox45°y45°z45°中，直接利用差分方法計(jì)算磁梯度張量：

G45°=

計(jì)算正六面體中心點(diǎn)P處的等效衰減磁矩μP時(shí)，由于正六面體中心點(diǎn)處的磁梯度張量無(wú)法按照差分方法計(jì)算得到，在此利用全微分的概念，建立微分方程組間接求解磁梯度張量各分量[12]。將圖2中的8個(gè)磁力儀按交叉線分成4組：S1和S8、S2和S7、S3和S6、S4和S5. 這4組磁力儀各自的連線均過(guò)正六面體中心點(diǎn)P. 以磁力儀S1和S8這一組為例，磁力儀S1和S8的測(cè)量數(shù)據(jù)差值為

ΔB81=(ΔB81x,ΔB81y,ΔB81z)=

(B8x-B1x,B8y-B1y,B8z-B1z).

(20)

根據(jù)全微分的概念，結(jié)合S1、S8的空間位置，有

(21)

結(jié)合其他3組磁力儀數(shù)據(jù)，則有方程組：

(22)

求方程組(22)式即可得到Gxx、Gxy和Gxz. 類(lèi)似地，利用ΔBy和ΔBz的全微分表達(dá)式，可計(jì)算得到Gyx、Gyy、Gyz、Gzx、Gzy和Gzz. 得到正六面體中心點(diǎn)P處的磁梯度張量G后，可計(jì)算出等效衰減磁矩μP.

3 仿真分析

仿真條件設(shè)置為磁偶極子磁矩M=(3 000 A·m2，3 300 A·m2，4 200 A·m2)，磁目標(biāo)的真實(shí)位置坐標(biāo)r=(-4 m，-6 m，-11 m)，計(jì)算得到的位置參數(shù)用rP表示，位置參數(shù)的相對(duì)誤差為

下面分析地磁場(chǎng)、磁梯度測(cè)量噪聲、正六面體基線以及磁力儀測(cè)量精度等因素對(duì)位置參數(shù)反演結(jié)果的影響，并與方法1的反演結(jié)果進(jìn)行比較。

3.1 地磁場(chǎng)對(duì)定位結(jié)果的影響

假設(shè)研究環(huán)境中地磁場(chǎng)大小為50 000 nT，地磁場(chǎng)梯度0.02 nT/m[13]，地磁傾角I=42°，磁偏角T=3°(設(shè)北偏東為正)。正六面體模型中的坐標(biāo)系x軸正向朝東，y軸正向朝北，z軸正向朝上。測(cè)量數(shù)據(jù)分兩種：一組數(shù)據(jù)中包含地磁場(chǎng)；一組數(shù)據(jù)中不含地磁場(chǎng)。根據(jù)設(shè)定的磁目標(biāo)參數(shù)，利用仿真的兩組數(shù)據(jù)計(jì)算磁目標(biāo)的位置參數(shù)，計(jì)算結(jié)果見(jiàn)表1.

表1 有無(wú)地磁場(chǎng)時(shí)的定位結(jié)果

從表1中的計(jì)算結(jié)果可以看出，以磁偶極子作為磁目標(biāo)，利用本文方法進(jìn)行磁目標(biāo)位置參數(shù)的反演計(jì)算時(shí)，地磁場(chǎng)對(duì)定位結(jié)果的影響很微弱，在應(yīng)用中地磁場(chǎng)的影響完全可以忽略。

3.2 磁梯度噪聲對(duì)定位結(jié)果的影響

在測(cè)量過(guò)程中，由于周?chē)拇艌?chǎng)環(huán)境比較復(fù)雜，各種設(shè)備及人為干擾難以避免，且各磁力儀之間還存在匹配誤差，使得測(cè)量的目標(biāo)磁場(chǎng)梯度信號(hào)中不可避免地包含一定的磁干擾量，測(cè)得的磁梯度信號(hào)不純凈。本文將這些磁干擾等效為磁梯度噪聲，添加到目標(biāo)磁信號(hào)測(cè)量數(shù)據(jù)中，磁梯度測(cè)量噪聲對(duì)定位結(jié)果的影響見(jiàn)圖4.

圖4中藍(lán)線是磁梯度測(cè)量噪聲水平為0.353 6 nT/m、隨機(jī)進(jìn)行100次計(jì)算得到的定位誤差曲線，黑線是磁梯度測(cè)量噪聲水平為0.707 2 nT/m、隨機(jī)進(jìn)行100次計(jì)算得到的定位誤差曲線。由圖4可見(jiàn)：磁梯度噪聲水平為0.353 6 nT/m時(shí)，定位誤差不超過(guò)1 m，大部分誤差在0.5 m以內(nèi)；磁梯度噪聲為0.707 2 nT/m時(shí)，定位誤差幾乎全在1.5 m以內(nèi)，大部分不超過(guò)1 m. 由于初始設(shè)置磁目標(biāo)的磁矩不大，在測(cè)量點(diǎn)上目標(biāo)信號(hào)的磁梯度大小約為幾十納特斯拉每米，相當(dāng)于測(cè)量信號(hào)中的磁梯度噪聲不超過(guò)目標(biāo)信號(hào)強(qiáng)度的10%，定位相對(duì)誤差也基本維持在10%以內(nèi)。

3.3 基線長(zhǎng)度對(duì)定位結(jié)果的影響

根據(jù)磁梯度張量G45°的表達(dá)式可知，采用差分替代微分的方法計(jì)算磁梯度張量時(shí)，基線長(zhǎng)度會(huì)影響G45°的計(jì)算結(jié)果，并最終影響位置參數(shù)的反演結(jié)果。定位誤差隨基線長(zhǎng)度變化的曲線見(jiàn)圖5.

從圖5可以看出，定位誤差隨正六面體邊長(zhǎng)基線的增加而變大。主要是因?yàn)橛?jì)算六面體側(cè)面的磁梯度張量時(shí)采用的是差分計(jì)算方法，與磁梯度張量的理論值之間存在一定的偏差，差分基線越長(zhǎng)，偏差越大，定位誤差也越大。在本節(jié)設(shè)置的條件下，當(dāng)邊長(zhǎng)基線不超過(guò)1.5 m時(shí)，定位誤差在0.7 m以內(nèi)。

3.4 磁力儀測(cè)量精度對(duì)定位結(jié)果的影響

當(dāng)磁力儀測(cè)量精度不同時(shí)，從表2中的定位結(jié)果可以看出，測(cè)量精度從0.01～1 nT逐漸降低時(shí)定位誤差逐漸加大，當(dāng)磁力儀的測(cè)量精度為1 nT時(shí)，定位誤差在1.4 m左右，相對(duì)誤差約為11%.

表2 不同測(cè)量精度時(shí)的定位誤差

3.5 本文方法與方法1的比較

表3 本文方法與方法1的計(jì)算結(jié)果比較

由表3可知，在磁力儀測(cè)量精度相同情況下，利用本文方法反演計(jì)算得到的位置參數(shù)與真實(shí)的位置參數(shù)更加接近，定位誤差為0.101 m，相對(duì)誤差為0.77%. 利用方法1反演計(jì)算磁目標(biāo)位置參數(shù)時(shí)，定位誤差為0.418 m，相對(duì)誤差為3.18%. 由此可見(jiàn)，本文方法的定位精度高于方法1. 另外，在磁梯度計(jì)算基線相等的情況下，本文正六面體模型的空間體積僅為方法1中模型體積的35.3%，顯著降低了模型所占的空間體積。

4 結(jié)論

本文研究了在地磁場(chǎng)中利用磁目標(biāo)的等效衰減磁矩信息實(shí)現(xiàn)磁目標(biāo)定位的方法。通過(guò)計(jì)算分析了等效衰減磁矩μ不隨測(cè)量坐標(biāo)系的改變而變化的特征。根據(jù)這一特征，利用了坐標(biāo)系旋轉(zhuǎn)的方法將口字型布置的磁力儀陣列轉(zhuǎn)變成十字型，在差分基線相等的情況下使正六面體的邊長(zhǎng)縮減為原來(lái)的0.707倍。利用求解全微分方程組的方法計(jì)算得到正六面體中心點(diǎn)的磁梯度張量，解決了正六面體中心點(diǎn)磁梯度張量難以計(jì)算的問(wèn)題。對(duì)磁目標(biāo)位置參數(shù)進(jìn)行了反演，結(jié)果表明在地磁場(chǎng)中利用等效衰減磁矩的定位方法可實(shí)現(xiàn)對(duì)磁目標(biāo)的有效定位。與之前的計(jì)算方法相比，在計(jì)算磁梯度的差分基線相等的情況下，本文方法正六面體模型的體積減小了64.7%，且定位效果更好。

)

[1] 于振濤，呂俊偉，畢波，等．四面體磁梯度張量系統(tǒng)的載體磁干擾補(bǔ)償方法[J]．物理學(xué)報(bào)，2014，63(11)：110702-1-110702-6．

YU Zhen-tao, LYU Jun-wei, BI Bo, et al. A vehicle magnetic noise compensation method for the tetrahedron magnetic gradiometer[J]. Acta Physica Sinica，2014, 63(11): 110702-1-110702-6. (in Chinese)

[2] 呂俊偉，遲鋮，于振濤，等．磁梯度張量不變量的橢圓誤差消除方法研究[J]．物理學(xué)報(bào)，2015，64(19)：190701-1-190701-8．

LYU Jun-wei, CHI Cheng, YU Zhen-tao, et al. Research on the asphericity error elimination of the invariant of magnetic gradient tensor[J]. Acta Physica Sinica，2015, 64(19): 190701-1-190701-8. (in Chinese)

[3] 萬(wàn)成彪，潘孟春，張琦，等. 基于張量特征值和特征向量的磁性目標(biāo)定位[J]. 吉林大學(xué)學(xué)報(bào)：工學(xué)版，2017，47(2)：655-660．

WAN Cheng-biao, PAN Meng-chun, ZHANG Qi, et al. Magnetic object location with eigenvalue and eigenvector of tensor[J]. Journal of Jilin University: Engineering and Technology Edition, 2017, 47(2):655-660. (in Chinese)

[4] 江勝華，申宇，褚玉程．基于磁偶極子的磁場(chǎng)梯度張量縮并的試驗(yàn)驗(yàn)證及相關(guān)參數(shù)確定[J]．中國(guó)慣性技術(shù)學(xué)報(bào)，2015，23(1)：103-106．

JIANG Sheng-hua, SHEN Yu, CHU Yu-cheng. Experimental verification and related parameter's determination for magnetic gradient tensor contraction using magnetic dipole[J]. Journal of Chinese Inertial Technology, 2015, 23(1): 103-106. (in Chinese)

[5] 陳謹(jǐn)飛，張琦，潘孟春，等．基于正六面體結(jié)構(gòu)測(cè)量陣列的磁異常定位技術(shù)研究[J]．傳感技術(shù)學(xué)報(bào)，2012，25(8)：1088-1092．

CHEN Jin-fei, ZHANG Qi, PAN Meng-chun, et al. Research on geomagnetic anomaly localization based on cubic measurement array[J]. Chinese Journal of Sensors and Actuators, 2012, 25(8): 1088-1092. (in Chinese)

[6] Nara T, Ito W. Moore-penrose generalized inverse of the gradient tensor in Euler's equation for locating a magnetic dipole[J]. Journal of Applied Physics, 2014, 115(17): 17E504-1-17E504-3.

[7] Sui Y Y, Li G, Wang S L, et al. Asphericity errors correction of magnetic gradient tensor invariants method for magnetic dipole localization[J]. IEEE Transactions on Magnetics, 2012, 48(12):4701-4706.

[8] Wilson H. Analysis of the magnetic gradient tensor[R]. Canada: Defence Research Establishment Pacific, 1985: 5-13.

[9] Beiki M, Keating P, Clark D. Depth estimation of magnetic and gravity sources using normalized source strength calculated from gradient tensor[C]∥2012 SEG Annual Meeting.Las Vegas, CA，US: The Society of Exploration Geophysicists, 2012: 996.

[10] Clark D A. New methods for interpretation of magnetic vector and gradient tensor data I: eigenvector analysis and the normalised source strength[J]. Exploration Geophysics, 2012, 43(4):267-282.

[11] Beiki M, Clark D A, Austin J R, et al. Estimating source location using normalized magnetic source strength calculated from magnetic gradient tensor data[J]. Geophysics, 2012, 77(6):23-37.

[12] 賈文抖，林春生，孫玉繪，等．基于單個(gè)磁梯度計(jì)的磁目標(biāo)定位方法研究[J]．兵工學(xué)報(bào)，2017，38(8)：1572-1577．

JIA Wen-dou, LIN Chun-sheng, SUN Yu-hui, et al. Research on magnetic target location method based on a single magnetic gradiometer[J]. Acta Armamentarii, 2017, 38(8): 1572-1577. (in Chinese)

[13] 張光，張英堂，李志寧，等．載體平動(dòng)條件下的磁梯度張量定位方法[J]．華中科技大學(xué)學(xué)報(bào)：自然科學(xué)版，2013，41(1)：21-24．

ZHANG Guang, ZHANG Ying-tang, LI Zhi-ning, et al. Localizing method of magnetic field gradient tensor under carriers moving parallelly[J]. Journal of Huazhong University of Science and Technology: Natural Science Edition, 2013, 41(1): 21-24. (in Chinese)