張 立，劉爭平

西南交通大學地球科學與環(huán)境工程學院地質工程系，成都 610031

1 引 言

自英國學者Rayleigh［1］發(fā)現(xiàn)瑞利面波的存在并揭示瑞利面波在彈性半空間介質中的傳播特性以來，由于瑞利（Rayleigh）面波采集簡便、面波速度VR與橫波速度VS相近、且具有速度頻散［2－3］等運動學特性，而被廣泛研究并應用于地殼結構、石油及工程等各種地下調查領域.如在工程應用中，StokeⅡ和Nazarian等［4］，提出了面波頻譜分析方法SASW，通過分析瑞利面波速度頻散曲線建立近地表S波速剖面.美國Kansas大學地質調查所的Xia等［2］提出了瑞利面波的多道分析方法MASW，反演估算近地表橫波波速.2003年，Xia等［5］又首次提出利用瑞利面波的基階和高階模式進行聯(lián)合反演的方法.目前，瑞利面波方法已成為近地表和大型土木工程基礎探測的重要方法.

瑞利面波理論和方法的研究主要集中在面波激勵場源和觀測參數(shù)方面.根據(jù)場源的激勵方式，瑞利面波可分為主動源面波法和被動源面波法［6］；按觀測參數(shù)的不同，可分為面波速度頻散特征的觀測和質點位移特征的觀測兩類，或者說是對面波運動學特征和動力學特征參數(shù)觀測的兩類.

傳統(tǒng)的瑞利面波的研究和應用主要集中在對主動源瑞利面波的運動學特征——速度頻散特征的研究和應用方面.速度頻散曲線的提取方法主要是利用頻率－波數(shù)（F－K）變換或慢度－頻率（P－F）變換［7］.近年來，研究學者利用空間相關技術（SPAC）［8－9］從地脈動信息中提取面波的速度頻散特征.但無論是主動源還是被動源面波法，提取速度頻散曲線都需要在觀測時布置多個具有一定道間距的檢波器的觀測排列，且勘探深度與排列長度成正比.因此，以觀測速度頻散等運動學特征為主的瑞利面波法，在實際應用中容易受場地限制，不適合人口稠密、空間狹小的城市及交通繁忙的公路等場區(qū).

1989年，Nakamura［10］首次提出地脈動信息的單點譜比法—H／V譜比或HVSR法.野外實際觀測表明單個三分量檢波器采集的地脈動信息的水平分量和垂直分量的頻譜比與場地VS結構之間具有較好的對應關系.Arai和Tokimatsu［11］及張立等［12］給出了利用地脈動的HVSR數(shù)據(jù)反演場地VS剖面的實例.2011年Haney等［13］對實際觀測到的被動源面波的橢圓極化進行了分析研究.HVSR法的工程應用展示了基于單點地脈動信息的動力學特征評價場地VS結構的實際可行性.但地脈動信息所包含波的主要成分具有爭議性［11－23］，如有研究學者認為其主要成分是基階瑞利面波［14－18］，有的認為是高階瑞利面波［11，21－22］，還有的認為是勒夫面波［19－20］，甚至有研究人員認為在某些條件下其主要成分是體波［24－25］，此外，再加上地脈動信息場源所具有的隨機性和復雜性，從而使得對地脈動隨機波場的正演數(shù)值模擬和理論研究不完備［26］等因素，導致HVSR法缺乏系統(tǒng)的應用理論基礎的支持.

本文基于瑞利面波理論解，對水平層狀介質中的基階瑞利面波的橢圓極化特征進行正演數(shù)值模擬，明確了瑞利面波的水平分量和垂直分量的頻譜比（HVSR）物理實質上表征的是瑞利面波的橢圓極化的頻散特征，該頻散特征與場地的泊松比結構密切相關.以期為發(fā)展基于單點瑞利面波的動力學特征的HVSR方法奠定理論基礎.

2 水平層狀介質中Rayleigh面波位移矢量

本節(jié)在朱介壽等［28］研究成果的基礎上，推導了水平層狀介質地表上Rayleigh面波水平位移ux和垂直位移uz計算公式，具體計算過程如下.

平面地震波在n層水平均勻介質組成的層狀空間中傳播，層狀空間如圖1［28］所示，設x軸沿測線方向，z軸垂直向下，各層及界面編號如圖1所示.

對于第m層，其參數(shù)分別為

圖1 水平n層介質模型［28］Fig.1 The model of n－layer horizontal media［28］

ux，uz為質點沿x，z方向的位移；σm，τm為分別法向應力和切向應力；波數(shù)k＝ω／c＝2π／λ，ω為角頻率.

由層狀介質中平面波傳播的波動方程：

式中，φ、ψ為縱、橫波位移位函數(shù)，可表示為

Δ′m、Δ″m、ω′m、ω″m為常系數(shù).當rαm為實數(shù)時，Δ′m表示與z軸正方向夾角沿arctanrαm方向傳播的平面波（下行波），即通過層狀介質傳播的體波；當rαm為虛數(shù)時，Δ′m表示沿x軸正方向傳播的面波，其振幅隨著深度增大而呈指數(shù)規(guī)律衰減.當rαm為實數(shù)時，Δ″m表示與z軸負方向夾角沿arctanrαm方向傳播的平面波（上行波），即通過層狀介質傳播的體波；當rαm為虛數(shù)時，亦表示沿x軸方向傳播的面波.對于ω′m、ω″m，只要將rαm改為rβm，也可按上述類似方式定義.令：

將（3）式中含有jkrz因子的指數(shù)函數(shù)改為三角函數(shù)，只考慮實數(shù)部分.并根據(jù)層狀介質內任一界面上umx，umz，τm，σm四個量連續(xù)條件，可得n層介質位移及應力的4n－2階的線性方程組：

其中，F(xiàn)R（ω，c）系數(shù)矩陣為面波的頻散方程.與位移有關的待求變量項：

即方程（5）為齊次方程組.為求解該齊次方程組，將系數(shù)矩陣FR（ω，c）、未知數(shù)向量X和右端矢量R進行如下處理：設第一層介質中波幅值A1已知，不失一般性，A1＝1.將未知矢量X和右端矢量R同除以A1，然后將其代入（5）式中，并將展開后的第一列，即與A1相關項，移到方程右端，經(jīng)過整理，方程（5）變?yōu)榉驱R次方程組（為方便起見，新方程以F′R（ω，c）、X′、R′表示）：

此時，方程（8）中系數(shù)矩陣F′R（ω，c）為4n－3階方陣.

圖2 質點位移及橢圓極化特性正演數(shù)值模擬研究程序流程框圖Fig.2 The numerical simulation flow diagram for the elliptic polarization characters of particle displacement

3 基階面波橢圓極化頻散特征數(shù)值模擬分析

根據(jù)以上求解，可求得當z＝0時，地表處不同頻率fj條件下，基階面波質點位移的水平分量和垂直分量umx｜z＝0（fj）＝ux（fj）、umz｜z＝0（fj）＝uz（fj）.

基于以上推導的理論公式和分析，本文對以下常見的三種工程地震勘探地層模型［29－32］：（1）均勻半無限空間地層模型，（2）正向結構地層模型（VS隨著埋深增大而增大）和（3）含有軟弱夾層結構的地層模型的橢圓極化頻散特性進行數(shù)值模擬研究.

3.1 均勻半無限空間地層模型

該組地層模型的泊松比結構特征如表1所示.由于一般介質的泊松比在0～0.5之間變化，如流體的泊松比0.5，軟沉積物的泊松比0.45，大多數(shù)巖石的泊松比為0.25，極堅硬的剛性巖石泊松比為0.05，即介質的彈性及強度越大，泊松比越小.因此，設計模型泊松比υ在［0.25，0.48］之間變化，橫波波速VS隨著泊松比的變化而變化，縱波波速VP，密度ρ等其它參數(shù)為常數(shù).圖3所示為該組模型中VP＝600m·s－1，υ＝0.25，VS＝346m·s－1的均勻半無限空間模型.

表1 均勻半無限空間地層模型參數(shù)Table 1 The model parameters for the half infinite homogeneous space

圖3 均勻半無限空間地層模型VP＝600m·s－1，VS＝346m·s－1，υ＝0.25Fig.3 The half infinite homogeneous space stratum model

對該組模型的橢圓極化頻散特征數(shù)值模擬結果見圖4.圖中表明，在泊松比一定的均勻半無限空間地層模型中，瑞利面波的橢圓極化頻散曲線是水平直線，即泊松比一定時，橢圓極化頻散特｜ux／uz｜｜z＝0（f）是一與f無關的常數(shù)，并且對于泊松比不同的均勻半無限空間地層模型，｜ux／uz｜｜z＝0（f）值不同，表明在均勻介質中，瑞利面波不存在橢圓極化頻散現(xiàn)象.

此外，對于泊松比不同的均勻半無限空間地層模型，隨著泊松比的增大，橢圓極化頻散｜ux／uz｜z＝0（f）的值逐漸減小，即隨著泊松比的增大，介質越松軟，橢圓越趨于垂直極化現(xiàn)象.這不難用泊松比是介質橫向應變與縱向應變之比的物理定義來解釋.因此，｜ux／uz｜｜z＝0（f），或其的倒數(shù)實質上反映的是介質泊松比度量，具有定性評價介質彈性性質和強度等物理力學性質的應用前景.

表2 軟弱夾層厚度不同、埋深相同的地層模型參數(shù)Table 2 The parameters of the models with a interbeded soft layer into a background model of elastic parameter increasing with depth increasing．

3.2 軟弱夾層厚度不同、埋深相同的正向結構地層模型

該組地層模型以單層厚度1m、各層VS隨深度遞增的20層正向結構地層模型為背景，均含有埋深為5m、夾層厚度在0.1～2m之間變化的軟弱夾層.所有模型中，第一層參數(shù)均為VS＝126m·s－1，VP＝422m·s－1，ρ＝1.61g·cm－3，h＝1m；軟弱夾層的VS＝102m·s－1，VP＝342m·s－1，ρ＝1.6g·cm－3，h＝0.1～2m.隨著深度的增大，各層VP與VS近似符合一定的線性變化規(guī)律［29－33］，密度ρ1～ρ20在1.6～2.01g·cm－3之間呈線性變化［29－33］（表2）.圖5（a、b）分別顯示了‘正向結構地層模型’和‘模型2＿5’.模型橢圓極化頻散曲線｜ux／uz｜z＝0（f）和速度頻散曲線的數(shù)值模擬結果見圖6所示.

從圖6a可看出，該模型中瑞利面波橢圓極化頻散曲線不再是直線，而是隨著f的變化而變化.這表明，對于非均勻介質，瑞利面波存在橢圓極化頻散特性.并且隨著模型軟弱夾層厚度在一定范圍內變化，與模型相應的橢圓極化頻散曲線也隨之變化，模型的橢圓極化頻散曲線在具有明顯的正負單峰值.如圖6a中正向結構地層模型的｜ux／uz｜｜z＝0（f）曲線在f0＝6Hz附近達到正峰值，在f0＝12Hz附近有負峰值，模型2＿5的曲線在f0＝4Hz附近出現(xiàn)正峰，在主f0＝8Hz附近出現(xiàn)負峰.正、負單峰值對應的頻率隨著模型軟弱夾層厚度的增大均向低頻端移動.此外，橢圓極化頻散曲線｜ux／uz｜｜z＝0（f）異常幅值隨著軟弱夾層的厚度H變化，具有最大軟弱夾層厚度模型2＿5的｜ux／uz｜｜z＝0（f）曲線與正向結構地層模型曲線的差異最明顯，而隨著H→0，橢圓極化頻散曲線｜ux／uz｜｜z＝0（f）也趨近于正向結構地層模型的曲線.

圖6b為該模型速度頻散曲線數(shù)值模擬結果.為對比，圖6（c、d）和表3分別是相對于正向結構地層模型的｜ux／uz｜｜z＝0（f）曲線和速度頻散曲線的相對變化圖和變化量表.在該地層模型中，｜ux／uz｜z＝0（f）曲線的相對變化量較大，最大值在2%～20%之間，平均變化量在0.1%～1.1%之間；相應的速度頻散曲線的相對變化較小，最大值在0～0.6%之間，平均變化量在0～0.1%之間.如模型1＿5，它的｜ux／uz｜（f）曲線的相對變化量最大近20%，平均變化量1.04%，而速度頻散曲線的相對變化量最大為0.38%，平均變化量為0.07%.這表明橢圓極化頻散曲線對模型參數(shù)變化的靈敏度比速度頻散曲線更高.

圖5 正向結構地層模型和軟弱夾層地層模型2＿5（軟弱夾層厚2m、埋深5m）（a）正向結構地層模型VP1＝422m·s－1，VS1＝126m·s－1；（b）軟弱夾層地層模型2＿5.Fig.5 The model with a interbeded soft layer into a background of elastic parameter increasing with depth increasing.（a）The background model of elastic parameter increasing with depth increasing with VP1＝422m·s－1，VS1＝126m·s－1；（b）Model 2＿5with a interbeded soft layer.

表3 含軟弱夾層的理論模型相對正向結構模型的｜ux／uz｜｜z＝0（f）曲線和速度頻散曲線的相對變化量Table 3 The variation of the elliptic polarization dispersion and the velocity dispersion compared with the forward stratum models

3.3 軟弱夾層厚度相同、埋深不同的正向結構地層模型

該地層模型同樣以圖6a中的正向結構地層模型為背景，且所有模型各層的參數(shù)VP、VS及ρ的大小及變化規(guī)律均與第二組模型相同，且模型的名稱及意義同，但模型中軟弱夾層固定厚度為2m、埋深在3～30m之間變化，模型具體參數(shù)見表4.

圖7顯示了該組地層模型橢圓極化頻散｜ux／uz｜｜z＝0（f）和速度頻散數(shù)值模擬結果.從圖7a可看出，該模型中瑞利面波橢圓極化頻散曲線也隨著f的變化而變化，即｜ux／uz｜｜z＝0（f）也與f有關，該結果同樣證明對于非均勻介質，瑞利面波存在橢圓極化頻散特性，且橢圓極化頻散現(xiàn)象與軟弱夾層的埋深有關.

從圖7a還可發(fā)現(xiàn)，該組模型的橢圓極化頻散曲線同樣以正、負單峰值為異常特征.隨著軟弱夾層埋深的逐漸增大，異常的正、負單峰值先向低頻端而后向高頻端移動，｜ux／uz｜｜z＝0（f）曲線也隨之先偏離而后趨近于正向結構地層模型的橢圓極化頻散曲線.

此外，橢圓極化頻散曲線｜ux／uz｜｜z＝0（f）異常幅值隨著軟弱夾層的埋深變化，具有最?。ā?）或最大（→∞）軟弱夾層埋深的模型如2＿3和2＿18的橢圓極化頻散曲線趨近于正向結構地層模型曲線，而隨著埋深在一定范圍時如模型2＿10，橢圓極化頻散曲線與正向結構地層模型曲線的差異最明顯.

圖7b為該地層模型的速度頻散曲線數(shù)值計算結果.圖7（c，d）和表3分別為該組模型相對于正向結構地層模型橢圓極化頻散曲線和速度頻散曲線的相對變化圖和變化量表.圖7及表3數(shù)據(jù)顯示，在該組地層模型中橢圓極化頻散曲線的相對變化較大，最大值在3%～13%之間，平均值在0.4%～1.2%之間；相應的速度頻散曲線的相對變化較小，其最大值范圍在0～0.6%之間，平均變化在0～0.2%之間.如模型2＿3的橢圓極化頻散曲線相對變化量的最大值接近13%，平均變化量為1.14%，而該模型的速度頻散曲線相對變化量最大不到0.6%，平均變化量僅為0.16%.

上述研究表明，相對于正向結構地層模型而言，隨著軟弱夾層埋深在一定范圍內的變化，相應模型的橢圓極化頻散曲線也隨之變化，即曲線與模型之間為一一對應關系，且橢圓極化頻散曲線對模型變化的反應比速度頻散曲線更敏感.

4 結 論

本文以工程中常見的幾種典型地層模型（均勻半無限空間地層模型、正向結構地層模型和含有軟弱夾層的地層模型）為例，應用HVSR方法，對基階瑞利面波的正演理論基礎——面波質點位移水平分量與垂直分量之比隨頻率f的變化規(guī)律｜ux／uz｜z＝0（f），即基階瑞利面波的橢圓極化頻散特征進行了數(shù)值模擬研究，得到如下結論：

表4 軟弱夾層厚度相同、埋深不同的地層模型參數(shù)Table 4 The parameters of the interbeded soft layer models with fixed thickness and various burial depth

（1）均勻介質中基階瑞利面波不存在橢圓極化頻散現(xiàn)象，而在非均勻介質中基階瑞利面波存在橢圓極化頻散現(xiàn)象.

（2）以｜ux／uz｜｜z＝0（f）表征的橢圓極化頻散特性，或其倒數(shù)實質上是介質泊松比度量，具有定性評價介質彈性性質和強度等物理力學性質的應用前景.

（3）在本文研究模型范圍內，基階瑞利面波橢圓極化頻散曲線對地層VS結構反應的靈敏度比速度頻散曲線更高.

（4）瑞利面波動的橢圓極化頻散曲線計算僅需面波單點的多分量數(shù)據(jù)特性，展示了單點面波調查方法的潛在發(fā)展前景.

此外，需要說明的是：

（1）本文所有計算研究過程中不考慮介質的耗散作用，則質點運動不考慮出現(xiàn)橢圓斜極化［28］現(xiàn)象.

（2）本文所選模型是以工程勘探中常遇到的軟地層條件為背景，其波速較低，對于硬地層模型還需要做進一步研究.

（

）

［1］ Rayleigh L.On waves propagated along the plane surface of an elastic solid.Proc.Lond.Math.Soc.，1885，17（1）：4－11.

［2］ Xia J H，Miller R D，Park C B.Estimation of near surface shear wave velocity by inversion of Rayleigh waves.Geophysics，1999，64（3）：691－700.

［3］ Dorman J，Ewing M.Numerical inversion of seismic surface wave dispersion data and crust－mantle structure in the New York－Pennsylvania area.JournalofGeophysicalResearch，1962，67（13）：5227－5241.

［4］ Stoke II K H，Nazarian S.Effectiveness of ground improvement from spectral analysis of surface waves／／Proc.8th.Euro Conf on Soil Mech and Found Engrg.Cambridge：Cambridge University Press，1983.

［5］ Xia J H，Miller R D，Park C B，et al.Inversion of high frequency surface waves with fundamental and higher modes.JournalofAppliedGeophysics，2003，52（1）：45－57.

［6］ Geophysical Service.Active and Passive Surface wave Techniques，www.geovision.com.

［7］ Capon J.High－resolution frequency wave number spectrum analysis.Proc.IEEE，1969，57（8）：1408－1418.

［8］ Aki K.Space and time spectra of statinary stochastic waves with special reference to microtreors.Bull.Earthq.Res.Inst.，1957，35：415－456.

［9］ Hortencia Flores Estrella.Jorge Aguirre Gonzalez，SPAC：An alternative method to estimation earthquake site effects in mexico city.GeofisicaInternacional，2003，42（2）：227－236.

［10］ Nakamura Y.A method for dynamic characteristics estimation of subsurface using microtremor on the ground surface.Quart.Rept.RTRIJapan，1989，30：25－33.

［11］ Arai H，Tokimatsu K.S－wave velocity profiling by inversion of microtremor H／V Spectrum.Bull.Seism.Soc.Am，2004，94（1）：53－63.

［12］ 張立，劉爭平，熊自英.利用地脈動H／V譜比反演地層Vs結構剖面.工程地球物理學報，2009，6（S1）：4－5.

Zhang L，Liu Z P，Xiong Z Y.S－wave velocity profiling by inversion of microtremor H／V spectra.ChineseJournalof EngineeringGeophysics（in Chinese），2009，6（SI）：4－5.

［13］ Haney M M，Nies A，Masterlark T，et al.Interpretation of Rayleigh－wave ellipticity observed with multicomponent passive seismic interferometry at Hekla Volcano，Iceland.TheLeadingEdge，2011，30（5）：526－531.

［14］ Lachet C，Bard P Y.Numerical and theoretical investigations on the possibilities and limitations of nakamura′s technique.J.Phys.Earth，1994，42（5）：377－397.

［15］ Nogoshi M，Igarashi T.On the Amplitude Characteristics of Microtremor－part2（in Japanese with English abstract）.J.Seism.Soc.Japan，1971，24：26－40.

［16］ Konno K，Ohmachi T.Ground－motion characteristics estimated from spectral ratio between horizontal and vertical components of microtremor.Bull.Soc.Am.，1998，88（1）：228－241

［17］ Tokimatsu K，Shinzawa K，Kuwayama S.Use of shortperiod microtremors forVsprofiling.J.Geotech.Rngrg.ASCE，1992，118（10）：1544－1588.

［18］ Horike M.Inversion of phase velocity of long－period microtremors to the S－wave－velocity structure down to the basement in urbanized area.J.Phys.Earth.，1985，33（2）：59－96.

［19］ Matsushima T，Okada H.An exploration methof using microtremors（2）－An pxperiment to identify Love waves in long－period microtremors－.Proc.The 82nd SEGF Conf，1990：5－8.

［20］ Tokimatsu K，Arai H，Asaka Y.Deep shera－wave structure and earthquake ground motion characteristics in Sumiyoshi area，Kobe city，based on microtremor measurements.J.Struct.Constr.EngngAIJ，1997，491：37－45（in Japanese）

［21］ Yutaka Nakamura.On the H／V spectrum.The 14th World Conference on Earthquake Engineering.2008，Beijing，China.

［22］ Tokimatsu K，Miyadera Y.Characteristics of Rayleigh waves in microtremors and their relation to underground structures.J.Struct，Constr，Engng，AIJ，1992，439：81－87.（in Japanese）

［23］ Arai H，Tokimatsu K.Effects of Rayleigh and Love Waves on microtremors H／V spectra.Proc.12th World Conf.on Earthquake Engineering，2000，paper 2232，CD－ROM.

［24］ 郭明珠，謝禮立，凌賢長.彈性介質面波地脈動單點譜比法研究.巖土工程學報，2004，26（4）：450－453.Guo M Z，Xie L L，Ling X Z.Research on spectral ratio of horizontal to vertical component for elastic model and surface microtremors.ChineseJournalofGeotechnicalEngineering（in Chinese），2004，26（4）450－453.

［25］ 郭明珠，謝禮立，閆維明等.體波地脈動單點譜比法研究.巖土力學，2003，24（1）：109－112.

Guo M Z，Xie L L，Yan W M，et al.Spectral ratio of horizontal to vertical component research for body wave microtremors.RockandSoilMechanics（in Chinese），2003，24（1）：109－112.

［26］ 劉爭平，張立，閆勝.局部非均質體的瑞利面波動力學響應數(shù)值模擬研究.工程地球物理學報，2009，6（S1）：1－3.

Liu Z P，Zhang L，Yan S.The dynamic feature analysis of Rayleigh surface wave on a local inhomogeneity using the numerical method.ChineseJournalofEngineeringGeophysics（in Chinese），2009，6（S1）：1－3.

［27］ Thomson W T.Transmission of elastic waves through a stratified solid medium.JournalofAppliedPhysics，1950，21（2）：89－93.

［28］ 朱介壽.地震學中的計算方法.北京：地震出版社，1988.

Zhu J S.Calculation Method in Seismology（in Chinese）.Beijing：Seismological Press，1988.

［29］ 陶蒞寧，侯廣春.波速測試技術在工程地質勘察中的應用.山東煤炭科技，2008，（2）：20－21.

Tao L N，Hou G C.Application of wave－speed measure technology in engineering geology survey.ShandongCoal ScienceandTechnology（in Chinese），2008，（2）：20－21.

［30］ 任健.多道瞬態(tài)面波法檢測技術在工程勘察中的應用.陜西水利水電技術，1999，（1）：45－51.

Ren J.Application of multi－channel transient surface wave method in engineering survey detection.ShaanxiWater ConservancyandHydropowerTechnology（in Chinese），1999，（1）：45－51.

［31］ 趙明杰，張桂玉，劉賓.瑞利面波在地質勘察中的應用.華北地震科學，2005，23（3）：31－33.

Zhao M J，Zhang G Y，Liu B.Application of Rayleigh wave in geological exploration.NorthChinaEarthquakeSciences（in Chinese），2005，23（3）：31－33.

［32］ 杜正濤，劉麗敏，程道偉.瑞雷面波勘探技術在第四系分層方面的應用.物探與化探，1999，23（4）：277－282.

Du Z T，Liu L M，Cheng D W.The application of Rayleigh surface wave Exploration Technique to the quaternary stratification.Geophysical＆GeochemicalExploration（in Chinese），1999，23（4）：277－282.

［33］ 陳仲候，王興泰，杜世漢.工程與環(huán)境物探教程.北京：地質出版社，1999.

Chen Z H，Wang X T，Du S H.Engineering and Environment Geophysical Course（in Chinese）.Beijing：Geological Publishing House，1999.