黏彈性層狀介質(zhì)中瑞雷波的本征值求解

劉雪峰，凡友華，常冬梅

(1.哈爾濱工業(yè)大學(深圳) 理學院,廣東 深圳518055；2.中國民航大學 航空工程學院, 天津 300300；3. 天津市高速切削與精密加工重點實驗室(天津職業(yè)技術師范大學)， 天津 300222)

黏彈性層狀介質(zhì)中瑞雷波的本征值求解

劉雪峰1,2，凡友華1，常冬梅3

(1.哈爾濱工業(yè)大學(深圳) 理學院,廣東 深圳518055；2.中國民航大學 航空工程學院, 天津 300300；3. 天津市高速切削與精密加工重點實驗室(天津職業(yè)技術師范大學)， 天津 300222)

為解決黏彈性層狀介質(zhì)中瑞雷波本征值求解中的漏根問題，提出了一種新的求解方法.通過分析以往方法中漏根和多根的原因，提出采用一系列初值虛部進行搜根后將結果合并，并去掉虛部大于實部的根.由此給出了本征值求解的流程圖，在實例計算中采用此方法成功求解了黏彈性層狀介質(zhì)中瑞雷波的本征值.結果表明，在利用該方法在求解過程中不存在漏根和多根問題.在黏彈性層狀介質(zhì)中瑞雷波的本征值求解時采用搜根初值加密并剔除虛部大于實部的根的方法可以有效避免漏根和多根問題.關鍵詞: 瑞雷波；黏彈性；層狀介質(zhì)；本征值問題；優(yōu)化

求解瑞雷波頻散的本征值問題是瑞雷波理論和應用研究的重要基礎，1953年，Haskell[1]改進了Thompson[2]的算法，提出傳遞矩陣算法，可以快速、有效地求解彈性層狀介質(zhì)中瑞雷波本征值.此后，Knopoff[3]和Abo-Zena[4]分別以此為基礎提出改進算法，提高計算速度和高頻數(shù)值穩(wěn)定性.凡友華等[5]在此基礎上提出了快速標量傳遞算法，提高了計算速度.1993年，Chen[6]在Kennett[7]的反射/透射系數(shù)法基礎上提出了廣義反射/透射系數(shù)法，從根本上解決了以往算法中高頻數(shù)值不穩(wěn)定的問題.在這之后也有學者[8-9]將該算法修正，分別提高了穩(wěn)定性和速度.

以往關于瑞雷波研究更多的是基于彈性假設.但實際上很多介質(zhì)都表現(xiàn)出一定的黏彈性，因此有必要研究黏彈性層狀介質(zhì)中瑞雷波的理論及應用，其中重要的基礎就是本征值問題求解.過去往往用Muller法[10]或牛頓法[11]進行求解.例如，Cao等[12]在2010年提出了一種將Muller法與牛頓法結合使用的Love波本征值問題求解方法，但是每個頻率下本征值計算要依賴于前一個頻率的本征值.而張凱[13]在2011年采用Muller法進行了求解.除此之外，Lai等[14]于2004年基于復分析中的柯西留數(shù)定理提出了一種求解方法.然而，這些方法都普遍存在著漏根及多根的問題，尤其在高頻下更為嚴重.本文將研究黏彈性層狀介質(zhì)中瑞雷波本征值求解中的漏根和多根問題，并基于Muller法給出一種求解方法以避免漏根和多根問題.

1 本征值求解的基本方法

瑞雷波頻散的本征值問題需要求解一系列的本征值，且要避免漏根問題，因此采用簡單的有限元方法或者求解奇異值方法并不適用.目前已有很多種不同形式的彈性層狀介質(zhì)中瑞雷波本征函數(shù)，其本質(zhì)上是等價的，只是在求解對應的本征值問題時在數(shù)值穩(wěn)定性及速度上有差異.Christensen[15]證明，對于黏彈性情況，瑞雷波的本征函數(shù)的形式是相同的，只是函數(shù)中原本取值為實數(shù)的彈性模量變?yōu)槿椭档酿椥越橘|(zhì)復模量.因此，黏彈性層狀介質(zhì)中瑞雷波的頻散方程是一個高度非線性的復函數(shù)方程，這一類方程的求解比彈性情況下的實方程更復雜.黏彈性層狀介質(zhì)中瑞雷波的本征函數(shù)是一個高度非線性的復函數(shù)，因此它的求解比彈性情況更加復雜.而且由于瑞雷波發(fā)生頻散后的同頻率下往往有若干個衰減系數(shù)和相速度，且彼此間差別有時較小，因此采用簡單的將本征值復根的實部和虛部分別細分并進行逼近不僅計算量無法接受，而且很容易發(fā)生漏根現(xiàn)象.Muller法可以用于這類方程的求解，下面簡單敘述利用Muller法求解該方程的過程.

設黏彈性層狀介質(zhì)中瑞雷波的頻散方程為F(c)=0，其中F(c)即為本征函數(shù)，c為本征值，即瑞雷波的復速度.Muller法對方程求解需要給出3個初值(c0,F(c0) ), (c1,F(c1) ), (c2,F(c2) )，然后利用以下迭代公式計算下一個值c3為

(1)

其中：

在得到c3后，令c0=c1,c1=c2,c2=c3，由此可以再計算下一個迭代值c3.直到函數(shù)值F(c3)的絕對值滿足所求精度或達到最大迭代次數(shù).

總的來說，對于一個給定的頻率，為求解頻散方程，大體計算步驟如下：

3)令c0=Vj+δi,c1=Vj+1+δi,c2=Vj+2+δi，其中δ為人為給定的常數(shù).利用式(1)計算c3，并利用Muller法進行迭代計算.

4)當函數(shù)值F(c3)的絕對值滿足所求精度或達到最大迭代次數(shù)時，結束迭代.通過判斷，若c3的值是可以接受的，則將其保存.

5)令j=j+1，繼續(xù)從步驟3)開始進行.

如此由欲求的最小頻率fmin到最大頻率fmax依次進行迭代計算，而每個頻率下的求解結果不受其他頻率下求得本征值的影響.

2 漏根現(xiàn)象及其解決方法

文獻[13]利用Muller法對黏彈性層狀介質(zhì)中瑞雷波的本征值問題進行了求解.在這里首先進行相同的計算，令m=100,δ等于當前頻率下各層橫波復速度的最大虛部.最大迭代次數(shù)設為20，頻散方程形式采用快速標量傳遞算法的頻散方程.該頻散方程具有較快的計算速度，由于求解黏彈性介質(zhì)中瑞雷波的本征值問題的計算量較大，因此采用該算法可以在一定程度上加快計算速度.在計算的穩(wěn)定性上，該算法與其他算法相比沒有本質(zhì)的差別，不影響本文的研究.以表1中的模型為例(參數(shù)與文獻[14]中的例子相同)，利用上述方法計算頻散曲線如圖1所示.顯然，本征值求解過程中發(fā)生了較為嚴重的漏根現(xiàn)象，尤其是在較高頻率區(qū)域漏根現(xiàn)象更嚴重.瑞雷波本征值求解中的漏根現(xiàn)象會對瑞雷波的反演應用造成巨大的影響，因此有必要對此進行深入的研究.

表1 黏彈性層狀介質(zhì)的參數(shù)

圖1 利用舊方法計算得到的頻散曲線

本文通過研究發(fā)現(xiàn)，增加迭代的最大步數(shù)可以得到更多的基階模(即第1條頻散曲線)頻散點，增大m的取值也可以在一定程度上得到更多的頻散點.除此之外，采用其他形式的頻散方程后漏根的區(qū)域也會有所不同.然而，這些處理措施都不能完全解決本征值求解中的漏根問題.

本文研究被漏掉的根附近的迭代步驟，發(fā)現(xiàn)有時迭代收斂到一個遠離初值的其他根，有時收斂到局部極小值點，由此發(fā)生了漏根.而這兩種現(xiàn)象的根本原因在于初值虛部的系數(shù)δ選取不合適，通常情況下都是由于δ太大造成的.由此，選取較小的δ進行求解.令c0=Vj+aδi,c1=Vj+1+aδi,c2=Vj+2+aδi，其中δ仍然等于當前頻率下各層橫波復速度的最大虛部，而令a=0.5.在新的初值下進行搜根，得到的頻散曲線如圖2所示.可以看到，在圖1中基階模以及其他一些模式的部分被漏掉的頻散點被成功計算出來.然而，其他一些圖1中原本能夠被搜到的頻散點反而被漏掉了.再次分析其原因發(fā)現(xiàn)這是由于初值虛部選取不夠大造成的.而還有一些點在圖1和圖2中都被漏掉，說明在計算這些頻散點對應的根時，初值虛部的選取仍然不合適.

圖2 將初值虛部減小后計算得到的頻散曲線

Fig.2 Dispersion curves obtained with initial values that have smaller imaginary part

圖3 選取10個不同初值虛部后計算得到的頻散曲線

Fig.3 Dispersion curves obtained with 10 initial values that have different imaginary part

3 多根現(xiàn)象及其解決方法

由圖3可以看到，盡管幾乎所有的頻散點都被成功得到，但是卻出現(xiàn)個別明顯不應存在的點.例如7 Hz頻率的第1個點和8 Hz頻率的第2個點.以7 Hz頻率的第1個點為例，其對應的本征值為111.59+152.51i,，盡管虛部的系數(shù)遠大于其他正常根，但是這個根卻可以驗證是近似滿足頻散方程的.也就是說，在搜根中會發(fā)生多根現(xiàn)象，盡管不如漏根現(xiàn)象顯然，但是對瑞雷波的理論及應用仍然有較大的影響,而這一現(xiàn)象在以往文獻中并沒有學者進行研究.由于多出的根是近似滿足頻散方程的，因此在現(xiàn)有的精度條件下無法從數(shù)值上直接剔除，只能根據(jù)經(jīng)驗人為去除.通過大量實驗發(fā)現(xiàn)，多根的虛部通常大于實部，而正常的根通常虛部遠小于實部，因此在目前不了解多根形成機理的情況下，可以將虛部大于實部的根去除，采用大量實例驗證發(fā)現(xiàn)用這種方法一般不會將正常的根錯誤地去除.由此，總結上述搜根方法得到最終的搜根流程如圖4所示.

利用圖4中所示的計算流程，計算得到表1中模型對應的頻散曲線和衰減系數(shù)曲線如圖5所示.在本例中，采用的最大迭代次數(shù)為50，在區(qū)間[0.1, 2.0]中等距地選取20個值作為a的取值，m的值取為1 000.由圖5可見，頻散曲線和衰減系數(shù)曲線上所有的點都成功地計算出來，即使高頻下也沒有發(fā)生漏根現(xiàn)象，而多根現(xiàn)象也沒有出現(xiàn).本文在對其他模型進行計算時用本方法也能夠避免漏根和多根現(xiàn)象，只要在頻散點較密集時增加a的取值個數(shù)并增大m即可.因此，利用圖4中的流程進行計算是可以成功地求解黏彈性層狀介質(zhì)中瑞雷波的本征值問題的.

圖4 求解本征值的計算流程

圖5 由圖4中流程計算得到的頻散曲線和衰減系數(shù)曲線

Fig.5 Dispersion curves and attenuation curves obtained with the flowchart shown in Fig.4

4 結 論

1) 分析了Muller法求解黏彈性層狀介質(zhì)中瑞雷波本征值過程中漏根和多根的原因，發(fā)現(xiàn)漏根是由于搜根時初值的虛部與真實值差別較大而陷入局部極小值造成的，采用一系列的虛部初值分別進行搜根，將搜根結果結合可以有效避免漏根現(xiàn)象.

2)發(fā)現(xiàn)盡管多根的機理并不清楚，但多根的虛部通常大于實部，而正常的根虛部通常遠小于實部，因此采用將虛部大于實部的根剔除的方法可以有效去除多根.

3)由此給出了本征值求解的具體流程，通過實例計算發(fā)現(xiàn)利用該方法可以成功求解黏彈性層狀介質(zhì)中瑞雷波的本征值問題，并避免漏根和多根現(xiàn)象.

[1] HASKELL N A. The dispersion of surface waves on multilayered media[J]. Bulletin of the Seismological Society of America, 1953, 43(1): 17-34. http://www.bssaonline.org/content/43/1/17.short.

[2] THOMSON W T. Transmission of elastic waves through a stratified solid medium[J]. Journal of Applied Physics, 1950, 21(2): 89-93. DOI: 10.1063/1.1699629.

[3] KNOPOFF L. A matrix method for elastic wave problems[J] . Bulletin of the Seismological Society of America, 1964, 54(1): 431-438.

[4] ABO-ZENA A M. Dispersion function computations for unlimited frequency values[J]. Geophysical Journal Royal Astronomical Society, 1979, 58(1), 91-105.DOI: 10.1111/j.1365-246X.1979.tb01011.x.

[5] 凡友華, 劉家琦. 層狀介質(zhì)中瑞雷面波的頻散研究[J]. 哈爾濱工業(yè)大學學報, 2001, 33(5): 577-581. FAN Youhua, LIU Jiaqi. Research on the dispersion of rayleigh waves in multilayered media[J]. Journal of Harbin Institute of Technology, 2001, 33(5): 577-581.

[6] CHEN Xiaofei. A systematic and efficient method of computing normal modes for multi-layered half-space[J]. Geophysical Journal International, 1993, 115(2): 391-409. DOI: 10.1111/j.1365-246X.1993.tb01194.x.

[7] KENNETT B L N. Reflection rays and reverberations[J]. Bulletin of the Seismological Society of America, 1974, 64(6): 1685-1696.

[8] 何耀鋒， 陳蔚天，陳曉非. 利用廣義反射-透射系數(shù)方法求解含低速層水平層狀介質(zhì)模型中面波頻散曲線問題[J]. 地球物理學報, 2006, 49(4): 1074-1081.DOI: 10.3321/j.issn:0001-5733.2006.04.020. HE Yaofeng, CHEN Weitian, CHEN Xiaofei. Normal mode computation by the generalized reflection-transmission coefficient method in planar layered half space[J]. Chinese Journal of Geophysics, 2006, 49(4): 1074-1081.DOI: 10.3321/j.issn:0001-5733.2006.04.020.

[9] PEI Donghong, LOUIE J L, PULLAMMANAPPALLIL S K. Improvements on computation of phase velocities of Rayleigh waves based on the generalized R/T coefficient method[J]. Bulletin of the Seismological Society of America, 2008, 98 (1): 280-287. DOI: 10.1785/0120070057.

[10]MULLER D E. A method for solving algebraic equations using an automatic computer[J]. Mathematical Tables and Other Aids to Computation, 1956, 10(56): 208-215. DOI: 10.2307/2001916.

[11]YAU L, BEN-ISRAEL A. The Newton and Halley methods for complex roots[J]. The American Mathematical Monthly, 1998, 105(9): 806-818. DOI: 10.2307/2589209.

[12]CAO Zhengliang, DONG Hefeng. Attenuation dispersion of Love waves in a viscoelastic multilayered half-space[C]//Proceedings of the SEG Technical Program Expanded Abstracts 2010. Denver, Colorado: Society of Exploration Geophysicists, 2010： 1924-1928. DOI: 10.1190/1.3513469.

[13]張凱. 粘彈性介質(zhì)瑞雷波模擬及其頻散曲線反演[D].武漢：中國地質(zhì)大學，2011. ZHANG Kai. Modeling and dispersion curve inversion of Rayleigh waves in viscoelastic media[D]. Wuhan: China University of Geosciences, 2011.

[14]LAI C G, RIX G J. Solution of the Rayleigh eigenproblem in viscoelastic media[J]. Bulletin of the Seismological Society of America, 2002, 92(6): 2297-2309. DOI:10.1785/0120010165.

[15]CHRISTENSEN R. Theory of viscoelasticity: an introduction[M]. New York: Academic Press, 1971, 245.

(編輯 張 紅)

Eigenproblem of Rayleigh wave in multilayered viscoelastic medium

LIU Xuefeng1,2, FAN Youhua1, CHANG Dongmei3

(1.School of Science, Harbin Institute of Technology (Shenzhen), Shenzhen 518055, Guangdong, China; 2.School of Aeronautical Engineering, Civil Aviation University of China, Tianjin 300300, China; 3. Tianjin Key Laboratory of High Speed Cutting and Precision Machining(Tianjin University of Technology and Education), Tianjin 300222, China)

To avoid the root lost problem in the former methods, a new method of solving the eigenproblem of Rayleigh wave in multilayered viscoelastic medium. The reason of root-lost and pseudo-root phenomenon is analyzed. From this, a series of different initial imaginary parts are employed, and the roots whose imaginary part is larger than real part are omitted. A flowchart of solving the eigenproblem is presented based on this, and the eigenproblem of Rayleigh wave in multilayered viscoelastic medium in an example is successfully solved with the method. The results show that there is no root-lost and pseudo-root phenomenon when solving the eigenproblem with the method. In the eigenproblem of Rayleigh wave in multilayered viscoelastic medium, root-lost and pseudo-root phenomenon can be avoided by employing a series of different initial imaginary parts and omitting the roots whose imaginary part is larger than real part.

Rayleigh wave; viscoelastic; multilayered medium; eigenproblem; optimization

10.11918/j.issn.0367-6234.201508063

2015-08-26

國家自然科學基金 (41204042)；中央高?；究蒲袠I(yè)務(3122015D009)；中國民航大學科研啟動基金(2014QD03S)

劉雪峰(1981—)，男，博士，講師; 凡友華(1975—)，男，教授，博士生導師

凡友華，yhfan@hit.edu.cn

O347.4

A

0367-6234(2017)04-0122-04