基于Schwab-Knopoff算法的優化改進

王國富，許崇彩，張法全

(桂林電子科技大學信息與通信學院，廣西桂林541004)

0 引言

由于瞬態瑞雷波勘探具有能量強、衰減小、輕便、快捷、高效、分辨率高、不受各地層速度關系影響等優點，廣泛應用在土層分層、地震勘探、淺層煤田勘探和地下煤層勘探[1－4]。研究表明:瑞雷面波在非均勻介質中傳播時產生頻散現象，即其傳播相速度VR是面波w的函數，F(w，VR)為面波的頻散函數，曲線w－F是面波的頻散函數曲線，而瑞雷波頻散曲線與地下介質結構密切相關，所以如何提取準確可靠的頻散曲線是瑞雷波勘探技術的關鍵問題之一[5－7]。

前人對瑞雷面波曲線的提取已作過一些研究。較好的瑞雷波頻散函數理論計算算法主要有Thomoson-Haskell算法、Schwab-Knopoff算法、δ矩陣算法和RT矩陣法。Thomoson-Haskell方法在Thomson研究基礎上，通過相鄰兩界面的傳遞矩陣公式以及自由表面邊界條件和無窮遠處的輻射條件推導出了層狀介質中平面瑞雷波的方程，但是此方法應用于高頻容易出現數值溢出及精度丟失問題[8];Schwab和Knopoff在Knopoff快速計算的基礎上進一步研究和改進，提出“歸一化”、“細分層法”，降低了有效數字的損失，計算溢出得到了改善，但是仍沒有完全消除[9－10]。文獻[11]采用的柱坐標系法降低了數值的不穩定，但并沒有完全消除;δ矩陣算法吸收了Schwab的變換程序思想，具有快速穩定的特點[12]，但是頻散曲線粗糙，尤其在高頻處曲線更加彎曲;RT矩陣法采用了反射和透射子矩陣技術，巧妙地避開了數值的不穩定性，但是這種方法要求所有的計算機均使用復數計算形式，故計算量大[13]。

針對Schwab-Knopoff容易造成數值溢出和丟失的缺陷，筆者在Schwab-Knopoff快速算法的基礎上提出了算法的簡化方法，并編程予以實現。

1 Schwab-Knopoff快速算法原理

Schwab-Knopoff快速算法的主要計算思路是根據上一層地質參數不斷由層參數遞推出新的矩陣元素，即由m層的矩陣元素推算出m+1層的矩陣元素，一直計算到n－1層界面為止。根據Schwab-Knopoff快速算法的推導，n層半無限空間介質的瑞雷波頻散函數為4n－2的行列式，頻散函數F(w，VR)最終形式為[14]:

2 Schwab-Knopoff算法改進

Schwab-Knopoff快速算法雖然應用比較廣泛，但是在高頻率計算頻散曲線時很容易出現高頻數值溢出的問題，對于實際的勘探有一定的限制，無法很好地解釋實際勘探中的問題，使反演受到了限制。仔細分析該算法原理公式時發現，有些計算因子會對數值產生很大的影響，頻率數值的一點變化就能導致頻散函數數值產生很大的變化(采用的層狀地質模型、地質參數如表1所示)。

表1 層狀介質各參數值Table 1 Layered media parameter values

當m為奇數時，公式的展開式為[14]

由以上公式推導可看出，矩陣函數AX為同一數量級的矩陣，矩陣函數BX除了第1行外，其他4行為kεiξj的形式。如果VR小于VSm，那么該層對應的矩陣含有cos(Pm)和cos(Qm)兩項，由于Pm=kγαmdm、Qm=kγβmdm(k是波數;d是厚度;γ為表達式)，那么為頻率的指數函數，當頻率比較高時，cos(Pm)、cos(Qm)為一很大的數量級，很可能超過計算機的最大數限制，出現無窮大，使計算不能順利進行，即出現高頻數值的溢出問題。本文對于層狀介質模型(介質參數見表1)計算了傳遞矩陣中用到的參數ξi(i=7，8，…，15)(公式迭代中運用到的元素為ξ7，ξ8，…，ξ15)、εi(i=1，2，…，15)值的最大和最小值，如表2所示。

表2 ξi、εi最大值和最小值Table 2 Maximum and minimum of ξiand εi

可以看出:當f=2 000 Hz、VR=100 m/s(不是高頻階段)，ξi(i=1，2，…，6)為無窮大，說明已經超出了計算機的數值上限，在搜根時可能會造成有效根的丟失;當f=5 000 Hz、VR=2 200 m/s時(情況與f=8 000 Hz，VR=2 100 m/s時相差不大)，ξi最小數量級為1013，而εi與ξi相差1012，相對于ξi來說，εi與0可以完全忽略不計;當f=5 000Hz，VR=3 000 m/s時，ξimin與εi在同一數量級，ξimax相對大一些，但不會造成很大影響。總之可以看出，當VR＜VS(i)時，在某段頻率內εi與ξi相差比較大，完全可以省略εi，使整體的數值減小，提高有效數字的長度，避免或者減小有效數字的丟失，提高精度，從而得到有效可靠的頻散函數圖形。當高頻時也能起到相同的提高精度的作用。所以可以采用如下的化簡方法:

因為ξm=kcos(Pm)cos(Qm)(m=7，8，…，15)，在高頻和低頻階段情況下，當VR小于VSm時，那么該層對應的矩陣含有cos(Pm)和cos(Qm)，為一大數量級。當cos(Pm)為大數量級時(設大于1040)，sin(Pm)=icos(Pm)，而矩陣的每個元素都可以表示為k1cos(Pm)cos(Qm)+k0，因為此時的cos(Pm)為非常大的數據，且k0相對于k1cos(Pm)cos(Qm)項完全可以忽略不計(由表2可見)，故在此層的傳遞矩陣中的元素可以化為cos(Pm)=1，sin(Pm)=i，k0=0，則完全不影響頻散函數的求根;同樣，當cos(Qm)為大數量級時，sin(Qm)必為大數量級，且有sin(Qm)=icos(Qm)，而矩陣的每個元素都可以表示為k1cos(Pm)cos(Qm)+k0，因為此時的cos(Qm)為非常大的數據，則cos(Pm)必為很大的數據(橫波波速VSm小于縱波波速VPm)且k0相對于k1cos(Pm)cos(Qm)項完全可以忽略不計，故在此層的傳遞矩陣中的元素可以化為cos(Pm)cos(Qm)=1，sin(Qm)sin(Pm)=－1，k0=0，則完全不影響頻散函數的求根，而且能夠提高有效數字的精度。

根據簡化算法，上述公式的元素簡化如下:

可見元素數值大為簡化，能夠減少數據的計算時間，并且εi與ξi能夠在一個數量級內(表3)，有效地減少有效數字的損失和數值的溢出。

表3 簡化后ξi、εi最大值和最小值Table 3 Maximum and minimum of simplified ξiand εi

3 改進算法的結果分析

實驗采用的層狀介質模型如表1所示。實驗模擬時，當計算數值大于1040時調用簡化算法，否則就用原算法。當頻散函數絕對值小于10就畫出該點，頻散函數絕對誤差如表4所示，0～5 000 Hz頻率范圍內的頻散函數分布如圖1、2所示。

從表4可以看出，原Schwab-Knopoff算法的函數絕對值偏大，有的超出范圍;在頻散函數曲線上(圖1)，超過3 500 Hz以后無法準確的分辨出正常的頻散曲線;而修改后的算法和δ矩陣法能夠得到準確而清晰的頻散曲線(圖2、3)，為反演及其地質層解釋打下良好基礎。

表4 頻散函數絕對誤差Table 4 Dispersion function of absolute error

圖1 原Schwab-Knopoff算法頻散圖Fig.1 Dispersion diagram of original Schwab-Knopoff

圖2 改進后的Schwab-Knopoff頻散圖Fig.2 Dispersion diagram of improved Schwab-Knopoff

從原算法頻散圖(圖1)還可以看出:該算法在低頻(＜1 000 Hz)時因為有效數字丟失或者數值溢出，導致頻散函數失真，不能準確提供信息;在中頻(1 000 Hz＜f＜3 500 Hz)由于數值在計算機數值限制范圍內，εi與ξi相差不大，能夠提供準確的頻散圖形;但是高頻(＞3 500 Hz)階段由于εi與ξi相差比較大，出現嚴重的數值精度丟失問題。這是因為計算中存在一些大數量級的量，在理論上能夠相減消去，但是數值上的誤差積累不能消去，即使有原算法的歸一化，也會導致數值溢出不能得到正確的結果，使計算不能順利進行，即出現了數值的溢出。這將導致頻散函數嚴重失真，頻率越大，高階頻散圖失真越厲害，無法提供可靠準確的頻散圖，導致后續工作無法正常進行。修改后的Schwab-Knopoff快速算法的頻散圖如圖2所示:該算法在低頻(＜1 000 Hz)時能夠很好地處理計算機數值溢出和有效數字，得到很好的頻散圖形;在中頻(1 000 Hz＜f＜3 500 Hz)和原算法的處理結果一樣;在高頻(＞3 500 Hz時)階段調用的簡化算法，能夠很好地處理高頻部分數值的溢出，從而確保了數值的精度，實現了正演算法的正確運行，為反演打下了很好的基礎。圖3為δ矩陣法的頻散圖，與改進后的算法頻散圖一致，驗證了算法的正確性，為地球物理勘探提供了準確的頻散曲線和理論研究。

4 結論

圖3 δ矩陣法的頻散圖Fig.3 Dispersion diagram of δ matrix method

在仔細研究原Schwab-Knopoff快速算法的基礎上，提出了該快速算法的簡化方案。實驗結果表明:該簡化算法簡單易行、數據簡單，不僅能夠較快速的計算實驗的頻散圖，還可解決低頻和高頻階段的上下限溢出和有效數字丟失的問題，使實驗頻散的有效頻率達到5 000 Hz以上，能夠提供準確的頻散函數曲線，尤其是清晰準確的高階頻散圖，為反演的復雜性提供準確的理論參考，并為瑞雷波的實際勘探提供強有力的理論基礎和依據，滿足實際的勘探需要，具有很大的發展潛力。

[1]Xia J H，Miller R D，Xu Y X，et al.High-frequency Rayleigh-wave method[J].Journal of Earth Science，2009，20(3):563－579.

[2]Lu J Q，Li S Y，Li W，et al.A hybrid inversion method of damped least squares with simulated annealing used for Rayleigh wave dispersion curve inversion[J].Earthquake Engineering and Engineering Vibration，2014，13(1):13－21.

[3]童立元，陳征宙，方磊，等.瞬態瑞雷波法在工程檢測中的應用[J].桂林工學院學報，1999，19(4):334－338.

[4]周文宗，熊章強，張大洲.基于軟弱夾層模型的瑞雷面波各模式相速度敏感性分析[J].桂林理工大學學報，2013，33(4):629－633.

[5]Xia J H，Xu Y X，Miller R D，et al.New developments in analysis of high-frequency Rayleigh waves[C]//Goephysical Solutions for Environment and Engineering:Proceedings of the 2nd International Conference on Environmental and Engineering Geophysics，2006，1:10－18.

[6]劉雪峰.層狀介質中瑞利面波多模式性質及其在正反演中的應用[D].哈爾濱:哈爾濱工業大學，2011.

[7]肖柏勛，李長征.瑞雷面波勘探技術研究述評[J].工程地球物理學報，2004，1(1):38－47.

[8]邵廣周.多階模式瑞利波頻散特征與反演研究[D].西安:長安大學，2009.

[9]鄧樂翔.瑞雷波場正演模擬及頻散曲線的提取[D].西安:長安大學，2010.

[10]焦健.基于空間自相關法的微動勘探技術的研究[D].長春:吉林大學，2012.

[11]張碧星.半無界分層介質表面任意面源激發的彈性波場[J].聲學學報，2007，32(3):193－204.

[12]楊威.瑞雷波在巖溶勘查中的應用研究[D].長沙:中南大學，2012.

[13]張凱.粘彈性介質瑞雷波模擬及其頻散曲線反演[D].武漢:中國地質大學，2011.

[14]曾校豐，崔偉群，許維進，等.瑞雷波頻散函數的快速算法[J].地球物理學報，2002，45(S1):261－267.