劉雪峰，凡友華，常冬梅

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

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

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

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

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

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

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

1 本征值求解的基本方法

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

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

(1)

其中：

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

總的來(lái)說(shuō)，對(duì)于一個(gè)給定的頻率，為求解頻散方程，大體計(jì)算步驟如下：

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

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

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

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

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

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

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

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

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

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

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

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

圖3 選取10個(gè)不同初值虛部后計(jì)算得到的頻散曲線

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

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

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

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

圖4 求解本征值的計(jì)算流程

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

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

4 結(jié) 論

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

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

3)由此給出了本征值求解的具體流程，通過(guò)實(shí)例計(jì)算發(fā)現(xiàn)利用該方法可以成功求解黏彈性層狀介質(zhì)中瑞雷波的本征值問(wèn)題，并避免漏根和多根現(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è)大學(xué)學(xué)報(bào), 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ì)模型中面波頻散曲線問(wèn)題[J]. 地球物理學(xué)報(bào), 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].武漢：中國(guó)地質(zhì)大學(xué)，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

國(guó)家自然科學(xué)基金 (41204042)；中央高校基本科研業(yè)務(wù)(3122015D009)；中國(guó)民航大學(xué)科研啟動(dòng)基金(2014QD03S)

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

凡友華，yhfan@hit.edu.cn

O347.4

A

0367-6234(2017)04-0122-04