受限觀測場地三維電阻率層析成像

中圖分類號：P319.3；P631.3 文獻標志碼：A

Abstract： Due to the presence of obstacles in the observation site or the unique characteristics of the target object，it is not always possble to lay out survey lines for electrical resistivity tomography in a straight manner.In some cases，even the electrode array cannot be placed directly above the target body. As a result，there are limitations on the application of resistivity tomography method under certain site conditions. This paper focuses on studying 3D electrical resistivity tomography technology for restricted Supported by the National Natural Science Foundation of China （4l774149）

observation sites.The study begins by designing three electrode aray observation systems specifically tailored for restricted observation sites：“” shaped observation system，“U” shaped observation system，and“L” shaped observation system. Subsequently，we investigate the 3D inversion imaging method for determining resistivity using any type of observation system. Finally，we analyze the detection effectiveness of diferent observation systems through model testing and apply“”shaped observation system to detect the base structure of an ancient tower in Zhejiang Province. Based on both geoelectric model detection results and the actual site conditions，allthree special observation systems for 3D electrical resistivity tomography can effectively capture 3D electrical characteristics below the restricted observation site， thereby providing valuable references for analyzing anomalies within the target area.Among them， the“□”shaped observation system demonstrates the best imaging performance，followed by the“U” shaped observation system，while the“L” shaped observation system shows relatively weak performance. In practical applications，it is essential to consider site conditions and exploration objectives when selecting the electrode array geometry and the number of electrodes，in order to optimize the imaging results.

Key words： restricted observation sites； unconventional observation system； electrical resistivity tomography

0 引言

電阻率層析成像法作為直流電阻率法的重要分支，通過在地表或鉆孔中布設電極陣列，可快速獲取地下目標體的地電信息，結合反演成像技術，能夠重建目標體的大致形態[1-2]。該方法廣泛應用于工程勘察、環境監測、地質災害探測及考古等領域，并取得了良好的應用效果[3-10]。電阻率層析成像方法在實際應用中通常采用直線或網格狀觀測系統，一般將電極陣列敷設在目標體上方，數據采集方便快捷，反演成像結果直觀，勘探效果良好[11-13]。但在調查建筑物底部或存在障礙物等受限觀測場景下，難以沿直線或網格狀布設測線，甚至也不能將電極陣列布設于目標體上方，導致方法應用受限，因此需要進一步拓展該方法技術。

為使電阻率層析成像法能夠適應復雜的觀測環境，研究人員從兩種思路展開研究：一種思路是開發新的電極，利用板狀電極代替棒狀電極或在棒狀電極添加輔助材料，以解決電極無法安插在類似水泥地、河流等特殊地表上的困境[14-16]；另一種思路是改變觀測系統的幾何形態，并重新設計相應的數據采集方式，從而達到既可以規避障礙物又不影響探測區域內測點覆蓋的效果。Loke等[12]設計了一種針對網格狀觀測系統的交叉對角線測量方法，僅沿0^°，45^° 和 90^° 發射和接收，可以使測量數據量減少一半，同時兼顧了良好的數據質量和分辨率。Brunner等[17]設計了一種陣列式同心圓環形觀測系統，并利用 SIRT （ simultaneous iterative reconstructiontechnique）法對火山口地質結構進行了電阻率成像。Tsourlos等[18]設計了一種放射狀陣列觀測系統，用以研究封土下方古墓的墓室結構，并與常規網格狀觀測系統進行了對比研究。Nyquist等[19]設計了網格狀布設電極陣列的放射狀單極-偶極觀測系統，使采集數據質量得到明顯改善。Argote-Espino等[20]在WennerSchlumberger（WS）陣列的基礎上設計了一種矩形環狀觀測系統，用以研究古墓遺址的地質結構，較好地解決了調查中涉及的問題。Tejero-Andrade等[21]設計了適應場地情況的“L\"形和“Corner\"形陣列觀測系統，解決了實際工程中由于道路障礙而難以布設電極的問題。綜上所述，目前電阻率反演成像技術在特殊觀測場地的應用越來越廣泛，因此針對復雜觀測場地，研究適應性強的陣列觀測系統和反演成像處理技術是十分必要的。

本文針對電阻率層析成像方法在實際復雜環境中所遇到的困難和障礙，首先設計了“□”形、“U”形、“L\"形等分布形態的電極陣列觀測方法，然后研究了電阻率三維反演方法，最后通過模型試驗和實際場地勘查檢驗方法的可行性，為受限場地開展電阻率勘探提供新途徑。

1 受限場地觀測系統設計

1.1 \"”形觀測系統

“\"形觀測系統是由多個電極圍繞某中心區域組成的一個環形網狀電極陣列，適用于地表存在障礙物或者需要調查某建筑底部區域而無法正常布設常規矩形觀測系統的受限場景。當障礙物四周均滿足布設電極的條件時，可以考慮采用“口\"形方式布設電極陣列。由于“口\"形觀測系統是在異常區域四周布設電極，因此需要事先確定探測目標的大致范圍，根據探測目標的范圍布設電極陣列。為取得較好的探測效果，需要圍繞探測目標布設多層電極陣列，以獲取盡可能多來自目標區的電性信息。

假如圍繞探測目標布設64根電極，并使其構成兩個封閉環，如圖1a所示。數據采集采用二極裝置AM ，為便于定性分析觀測系統在不同方向的分辨能力，可根據供電電極 A 的平面坐標 和測量電極 M 的平面坐標 （x_M，y_M） 換算出數據點的視記錄點坐標，即 x=（x_A+x_M）/2，y=（y_A+y_M）/2 z=AM/2 ，這樣可將所有數據點展布于地表下方。“\"形觀測系統的測點分布如圖1b—d所示。可以看出，“口\"形觀測系統測點在 x 和 y 方向上分布范圍廣且空間覆蓋均勻，在 z 方向上測點數量隨深度增大呈現先增加后減少的趨勢，其中，上中層區域測點密度較高且覆蓋面積大，確保了在橫向和縱向上對目標體對目標體的有效分辨能力。電極陣列環數越多，分辨能力越強。

1.2 “U”形觀測系統

當探測目標區一側無法布設電極時，可將電極按“U\"形方式布設。假如圍繞探測目標布設38根電極，并使電極陣列呈“U\"形分布，數據采集采用二極裝置，“U\"形觀測系統的電極布設和測點分布如圖2所示。從圖2b—d中可以看出，“U\"形觀測系統沿 x 方向測點覆蓋面積較大且分布相對均勻，而沿 y 方向測點存在缺失，對目標體的分辨能力不及“\"形觀測系統。

圖1“□\"形觀測系統和測點分布示意圖

Fig.1Distribution diagram of “” shaped observation system and its measuring points

圖2“U\"形觀測系統和測點分布示意圖

Fig.2Distribution diagram of“U” shaped observation system and its measuring points

1.3 “L\"形觀測系統

當探測目標區兩側無法布設電極時，可將電極按“L\"形方式布設。假如圍繞探測目標布設51根電極，并使電極陣列呈“L\"形分布，采用二極裝置測量，“L\"形觀測系統的電極布設和測點分布如圖3所示。從圖3b—d中可以看出，“L\"形觀測系統沿x 和 y 方向測點均有缺失，對目標體的分辨能力不及“口\"形和“U\"形觀測系統。

2 電阻率三維約束反演

已知電阻率的線性反演方程為態的，為保證反演過程的穩定性，需要對地電模型施加先驗約束（包括光滑和背景約束）。綜合方程（1）和先驗約束，在最小二乘意義下構建電阻率反演的目標函數[22]：

式中： m 為模型參數向量； m_b 為背景模型； c 為光滑度矩陣； 為單位矩陣； λ_c 和 λ_b 為正則化因子。

將方程（2）兩端對 Δm^T 求導：

2λ?cornerC^TC（m+Δm）+2λ?cornerI^TI（m+Δm-m_b）?corner

AΔm=Δd

式中： Δd 為數據殘差矢量； Δm 為模型參數的修正向量； A 為偏導數矩陣。

在電阻率三維反演中，方程（1）通常是欠定且病

令式（3）等于0，經整理可得電阻率反演的線性方程：

a.“L\"形觀測系統;b.測點空間分布；c.測點在 yOz 平面投影；d.測點在 xOz 平面投影。

圖3“L\"形觀測系統和測點分布示意圖

則方程（4）可簡化為

線性反演方程（5）具有滿秩、病態和規模大的特點，共軛梯度法在求解此類方程方面具有明顯優勢，被廣泛地應用于科學研究、工程計算等領域[23]。本文采用最小二乘共軛梯度法求解方程（5）[24]。首先令方程（5）的初始解 Δm₀=0 ，則方程（5）的初始殘差 r₀=b-GΔm₀=b ，梯度向量 和共軛方向向量p 的初始值為 g₀=p₀=G^Tr₀=G^Tb 。其迭代過程歸納如下。

對于 i=0，1，…，i_max ，計算

q_i=Gp_i，

Δm_i+1=Δm_i+a_i+1p_i，

r_i+1=r_i-a_i+1q_i，

g_i+1=G^Tr_i+1，

如果 |g_i+1|²/|g₀|²?ε ，則終止迭代。否則轉下一步：

p_i+1=g_i+1+β_i+1p_io

式中： i_max 為最大迭代次數； ε 為共軛梯度迭代終止誤差； a_i+1 和 β_i+1 分別為解向量 Δm 和共軛方向向量 p 的修正因子。設置 ， ε=10^-6 。為提高反演計算效率，對共軛梯度算法中涉及的For循環進行了OpenMP并行化。

通過上述方法求解方程（5），可得模型參數修正量 Δm ，將其代人下式：

m₁=m₀+Δm

即可得到新的模型參數向量 m₁ 。經多次迭代，直至迭代次數或數據擬合差滿足終止條件為止，此時 m₁ 即為預測的電阻率模型。

3 數據處理過程

受場地觀測條件限制，電極陣列通常是不規則的，甚至電法儀器也可能不是通用的三維數據采集設備。結合實際情況，下面給出電阻率三維反演的數據處理過程。

第一步：野外數據采集時，需要對電極陣列的所有電極進行順次編號 （1，2，...，n） ，并記錄每根電極的空間坐標 （x，y，z） 。電極陣列編排格式見表1。

表1電極陣列的編排格式

Table1 Format of electrode array

第二步：從電法儀器中導出視電阻率數據（因默認電極位置不準確，視電阻率數據換算錯誤），編輯每個視電阻率數據點對應的裝置類型、供電電極 A 和 B 的編號、測量電極 M 和 N 的編號。如果數據采集為二極裝置，裝置類型為2，供電電極 B 和測量電極 N 的編號均為0；如果數據采集為三極裝置，裝置類型為3，供電電極B的編號為0。采用三極裝置測量時，觀測數據的編排格式如表2。

表2視電阻率數據的編排格式

Table2Arrangementformatofapparentresistivitydata

第三步：根據電極陣列的實際空間坐標重新計算準確的視電阻率或阻抗。目前大多陣列式電法儀器按測線為直線、道間距相等的理想情況計算視電阻率。在沒有考慮測線彎折和道間距不相等等情況下，儀器按非實際坐標計算的視電阻率為

式中： K^′ 。為儀器按電極的非實際坐標計算的裝置系數； ΔU 為觀測電位差； I 為發射電流。對式（7）進行變換，可以獲取準確計算的視電阻率 ρ_s 和阻抗 R ：

式中， K 。為按電極的實際坐標計算的裝置系數。

第四步：根據換算的 ρ_s 或 R 反演地下電阻率ρ_c 。反演過程大致可分為輸入實測數據（電極陣列數據和準確計算的視電阻率或阻抗）、設置反演參數與構建反演初始模型、電阻率三維反演計算、輸出反演預測模型。具體反演成像流程如圖4所示。

圖4反演處理流程圖

Fig.4Flow chart for the processof inversion

4模型反演與應用案例

4.1 模型反演

假定在圍巖電阻率為 ρ=100Ω?m 的均勻半空間中存在兩個異常體，其電阻率分別為 ρ₁=10Ω · ，尺寸均為 2m×2m×2m ，兩異常體相距 2m ，頂板距地面均為 2m ，在地表將電極陣列分別布設成“口\"“U\"和\"L\"形態，相鄰電極間距為 1m ，如圖5所示。

基于這三種觀測系統對所建立的異常體模型進行三維反演，各觀測系統所采用的裝置類型、電極數目和模擬數據數目詳見表3。反演過程均迭代5次，其均方根誤差（rootmean squared，RMS）的收斂曲線如圖6a所示。可以看出，隨著迭代次數的增加，三種觀測系統的誤差均穩步下降，而且數據點越多（表3），收斂曲線的RMS越小，數據擬合越好。圖6b一d為模擬的視電阻率與最終反演模型計算的視電阻率數據擬合結果，可以看出，三種觀測系統的數據擬合圖形態相近，并且 ρ_s 與 ρ_c 的交匯點位基本分布于圖形對角線附近，符合地球物理反演數據擬合規律。

圖7為三種觀測系統的反演結果的切片圖。對比三種觀測系統的反演結果，“口\"形觀測系統的勘探效果最佳，反演異常體的平面分布與實際位置吻合較好（圖7a），但在深度上異常體呈現上移的情況（圖7b）。可能存在兩方面的因素：1）反演模型采用了分塊連續介質模型，目標體頂板 2m 處至其上方1m 處物性是線性過渡的；2）電極布設于目標體周圍，導致其頂部地電信息缺失。此外，由圖5a可知，兩異常體在 x 方向靠近電極，而在 y 方向遠離電極，使得“口\"形觀測系統的反演結果在 x 方向的分辨率高于 y 方向的分辨率（ x 方向梯度變化大于 y 方向梯度變化），說明分辨率高低與目標體和電極陣列的距離密切相關，兩者距離越近，探測效果越好。

圖5觀測系統的電極分布和異常體位置圖

Fig.5Electrode distribution of observation systems and location of anomalous bodies

表3不同觀測系統的模擬參數

Table3 Simulationparametersofdifferentobservation systems

對比分析“U\"形和“L\"形觀測系統的反演結果（圖7c—f可以看出，這兩種觀測系統對低阻體的定位結果總體良好，但高阻體的位置和形態出現了較大偏差，只能提供較彌散的高阻特征。這主要源于電阻率勘探方法對低阻體敏感的特性；另一方面，該低阻異常體距離電極陣列較近（圖5a—c），高阻異常體處于電極陣列的觀測盲區。

種觀測系統的誤差收斂曲線;b.“口\"形觀測系統的數據擬合圖;c.“U\"形觀測系統的數據擬合圖；d.“L\"形觀測系統的數據擬合|

圖6三種觀測系統的誤差收斂曲線和數據擬合圖

Fig.6Error convergence curves and data fitting plots for three observation systems

a 和b分別為\"□\"形觀測系統在 z=-2.5m 處的水平切片和 y=0m 處的垂直切片;c和d分別為\"U\"形觀測系統在 z=-2.5m 處的水平切片和 y=0m 處的垂直切片;e和f分別為“L\"形觀測系統在 z=-2.5m 處的水平切片和 y=0m 處的垂直切片。紅色線框表示異常體。

圖7三種觀測系統的反演結果切片圖

Fig.7Sliced images of the inversion results for three observation systems

4.2 應用案例

采用受限場地三維電阻率層析成像技術對浙江省龍游縣某古塔塔基進行探測。該塔始建年代已無從考證，重建于北宋嘉祐三年（公元1058年），現存為北宋塔身，塔高 27.31m ，塔身呈八角形。該古塔是中國佛教建筑史上的重要遺產，也是江南地區的文化遺產之一。它的建筑形制與風格對中國佛教建筑和塔式建筑的發展產生了深遠影響，被譽為“江南古塔之首”。本次探測工作主要查明塔基的基礎規模，為下一步保護提供參考依據。

采用受限場地三維電阻率層析成像技術對浙江某古塔塔基進行探測。電極陣列采用如圖1a所示的“□\"形布設，平面坐標（0，0）點基本為塔基中心位置。受觀測條件的限制，古塔周邊不方便布置無窮遠極，因此觀測裝置采用溫納alpha裝置。電極間距為 1m 。儀器設備默認測線為直線且點距為等間距，導致部分數據無法使用。首先剔除不正常的數據。電纜線環狀布設，部分供電電極對和測量電極對相對位置關系復雜，導致裝置系數非常大，需要手工剔除掉不能用的數據。野外實測數據為552個，剔除掉不能用的數據后，剩余有效數據345個；然后將有效的視電阻率數據換算成電位差，再對電位數據作反演，得到地下三維電阻率信息。圖8為本次探測工作的三維電阻率反演結果，可以看出，不同深度的電阻率水平切片總體呈中間高四周低的趨勢，由淺至深基本反映了塔基及其周圍的電阻率變化特征。從深度 1.6m 的電阻率切片（圖8b）中可以看出，高阻區基本位于塔基中心區域；而隨深度的增加，高阻異常中心有向北和西偏移的趨勢（圖8c、d）。這可能與建造塔基時所使用的材料有關。此外，古塔西側地表出露長條狀大理巖體，可能延伸至底部，進而導致塔心高阻異常出現偏移現象。

a.施工現場；b. z=-1.6m 電阻率切片結果；c. z=-2.5m 電阻率切片結果；d. z=-3.2m 電阻率切片結果。

圖8施工現場和不同深度的電阻率切片結果

Fig.8 Construction site and resistivityslice results at different depths

5結論

1）在觀測場地受限的條件下，可以嘗試采用形如“口”、“U\"和\"L\"等分布形態的電極陣列開展電法勘探工作。為提高受限觀測場地三維電阻率層析成像技術的探測效果，實際勘探中電極陣列布設應盡可能靠近目標體，采用盡可能多的電極數量，并且可以根據場地要求設計電極陣列的分布形態，確保采集到盡可能多的有效數據。

2）不同形態的觀測系統與目標體的耦合關系不同，從算例結果對比來看，“口\"形觀測系統的探測效果優于“U\"形和“L\"形觀測系統。在觀測場地十分苛刻的情況下，“U\"形和“L\"形觀測系統也能提供有意義的參考信息。

3）二極相比于三極和四極裝置，需要剔除的無

效數據較少，推薦使用。

參考文獻（References）：

[1] Shima H. 2 - D and 3 - D Resistivity ImageReconstruction Using Crosshole Data[J]. Geophysics，1992，57（10）：1270-1281.

[2] Paz-Arellano P，Tejero-Andrade A，Argote-Espino D.2D-ERT Survey for the IdentificationofArchaeological and Historical Structures Beneath thePlaza of Santo Domingo，Mexico City，Mexico[J].ArchaeologicalProspection，2017，24（3）：183-194.

[3] 呂惠進，劉少華，劉伯根.高密度電阻率法在地面塌陷調查中的應用[J].地球物理學進展，2005，20（2）：381-386.LüHuijin，Liu Shaohua，Liu Bogen.ApplicationofResistivity Tomography Survey Method in DetectingGIUuI20（2）：381-386.

[4]江玉樂，周清強，黃鑫，等.高密度電阻率法在滑坡探測中的應用[J].成都理工大學學報（自然科學版），2008，35（5）：542-546.Jiang Yule， Zhou Qingqiang，Huang Xin， et al.Application of the High Density Resistivity Method toLandslide Prediction [J]. Journal of ChengduUniversity of Technology （Science amp;.TechnologyEdition），2008，35（5）：542-546.

[5]尹極，王曉東，程邈，等.高密度電阻率法在果香峪溶洞勘查中的應用效果分析[J].CT理論與應用研究，2011，20（3）：301-310.Yin Ji， Wang Xiaodong，Chen Miao，et al. Analysis onApplication Effect of High Density Resistivity Methodto Guoxiangyu Solution Cavity Exploration[J]. CTTheory and Applications，201l，20（3）：301-310.

[6]Mauriello P，Monna D，Patella D. 3D GeoelectricTomography and Archaeological Applications[J].Geophysical Prospecting，1998，46（5）： 543-570.

[7]蘇永軍，王緒本，羅建群.高密度電阻率法在三星堆壕溝考古勘探中應用研究[J].地球物理學進展，2007，22（1）： 268-272.Su Yongjun， Wang Xuben， Luo Jianqun. TheAhaeological Application of High-Density ResistivityMethod to Ditch Exploration on Sanxingdui Site[J].Progress in Geophysics，2007，22（1）：268-272.

[8]戴田宇，謝尚平.高密度電阻率法在葬墓考古的應用[J].地球物理學進展，2015，30（6）：2885－2891.Dai Tianyu， Xie Shangping. Application of ResistivityTechnique in the Relic Archaeology[J]. Progress inGeophysics，2015，30（6）：2885-2891.

[9]余天祥，石戰結，劉杰，等.多方位擬3D 電法技術在古墓葬探測中的應用[J].CT理論與應用研究，2018，27（6）：727-738.Yu Tianxiang，Shi Zhanjie，Liu Jie，et al. Applicationof Multi-Azimuth Pseudo - 3D Electrical ResistivityTomographyinthe Investigation ofAncientMausoleums[J]. CT Theory and Applications，2018，27（6）：727-738.

[10]丁超，解陽波，張家豪，等.三維高密度電法揭示Optimized Arrays for 2-D Resistivity Survey Lineswith a Large Number of Electrodes[J]. Journal ofApplied Geophysics，2015，112：136-146.

[12]Loke M H，Barker R D. Practical Techniques for 3DResistivity Surveysand Data Inversionl [J].Geophysical Prospecting，1996，44（3）： 499-523.

[13]顧觀文，王順吉，李桐林，等.不同供電方向激勵對多極距組合中梯裝置三維激電探測效果的影響[J].吉林大學學報（地球科學版），2024，54（1）：292－309.Gu Guanwen， Wang Shunji， Li Tonglin， et al.Influence ofDifferentPowerSupplyDirectionExcitation on Three-Dimensional Induced PolarizationExploration Effect with Central Gradient Array ofMuti-Group Power Supply Pole Distance Combination[J]. Journal of Jilin University （Earth ScienceEdition），2024，54（1）：292-309.

[14]Athanasiou E N，Tsourlos P I，Vargemezis G N，etal. Non-Destructive DC Resistivity Surveying UsingFlat-Base Electrodes [J]. Near Surface Geophysics，2007，5（4）：263-272.

[15]Tsokas G N，Tsourlos P I，Vargemezis G，et al.Non-Destructive Electrical Resistivity Tomographyfor Indoor Investigation： The Case of KapnikareaChurch in Athens[J]. Archaeological Prospection，2008，15（1）：47-61.

[16]李彬，廣先國，張賡，等.高密度電法在礦山地質災害勘察中的應用[J].金屬礦山，2013，42（6）：87－89.Li Bin，Tuo Xianguo， Zhang Geng， et al. Applicationof High Density Resistivity Method to the Survey ofGeologic Disasters in Mine[J]. Metal Mine，2013，42（6）：87-89.

[17]Brunner I， Friedel S，Jacobs F，et al. Investigation ofa Tertiary Maar Structure Using Three-DimensionalResistivityImaging [J]. Geophysical JournalInternational，1999，136（3）：771－780.

[18]Tsourlos P， Papadopoulos N， Yi M J， et al.Comparison of Measuring Strategies for the 3 - DElectrical Resistivity Imaging of Tumuli[J]. Journalof Applied Geophysics，2014，101：77-85.

[19]Nyquist J E，Roth M J S. Improved 3D Pole-DipoleResistivity Surveys Using Radial Measurement Pairs[J].Geophysical Research Letters，2005，32（21）：44-48.

[20]Argote-Espino D， Tejero-Andrade A，Cifuentes-NavaG， et al. 3D Electrical Prospection intheArchaeological Site of El Pahnu， Hidalgo State，Central Mexico [J]. Journal ofArchaeologicalScience，2013，40（2）：1213-1223.

[21] Tejero-Andrade A，CifuentesG，Chavez RE，etal.L- and CORNER-Arrays for 3D Electric ResistivityTomography： AnAlternative forGeophysicalSurveys inUrbanZones [J].NearSurfaceGeophysics，2015，13（4）：355-368.

[22] 劉海飛，柳建新，麻昌英.直流激電反演成像理論與方法應用[M].長沙：中南大學出版社，2017.Liu Haifei，Liu Jianxin， Ma Changying. InverseTheory and Methodology for Direct Current InducedPolarization Method[M].Changsha：Central SouthUniversityPress，2017.

[23] 劉海飛，柳建新，柳卓.數值計算與程序設計（地球物理類）[M].長沙：中南大學出版社，2021.LiuHaifei，Liu Jianxin，Liu Zhuo.NumericalComputation and ProgramDesign（GeophysicsCategory）[M].Changsha：Central South UniversityPress，2021.

[24] 劉海飛，阮百堯，柳建新.變阻尼共軛梯度算法及其性能分析[J].地球物理學進展，2008，23（1）：89－93.Liu Haifei，Ruan Baiyao，Liu Jianxin. The ConjugateGradient Algorithm of Variational Damping and ItsCapability Analysis[J]. Progress in Geophysics，2008，23（1）：89-93.