畢業(yè)論文 -一種新的頻率選取策略及其在頻域波形反演中的應(yīng)用

1、PAGE 4一種新的頻率選取策略及其在頻域波形反演中的應(yīng)用A new strategy for selecting temporal frequencies and its application in frequency-domain seismic waveform inversion目 錄摘要ABSTRACT1緒論1.1研究背景和意義1.2國內(nèi)外發(fā)展現(xiàn)狀1.3主要研究內(nèi)容2頻域全波形反演算法2.1正演模擬2.2反演算法3頻域波形反演中的波數(shù)覆蓋原理3.1散射場的積分表示3.2散射場的一階Born近似3.3高頻近似3.4常數(shù)背景和變化背景下散射場Born近似的漸近公式3.5波數(shù)覆蓋4線性遞

2、增速度介質(zhì)中地震波的傳播5頻率點的選取方法摘要頻域全波形反演是一種能夠獲得高分辨率反演結(jié)果的反演方法，但在模型較大時進(jìn)行全頻點反演的計算量很大。Sirgue和Pratt (2004)提出一種波形反演中的頻率點選取的方法，在不影響反演精度的前提下，可以減少頻域波形反演所需計算的頻率點數(shù)目，從而節(jié)省計算量。目前，大多數(shù)的頻率點選取方法都是基于地下介質(zhì)速度為均勻分布這一假設(shè)前提的，為了研究在更加接近實際地下介質(zhì)的背景模型下的情形，本文在線性遞增速度模型的基礎(chǔ)上，利用散射波場的Born近似和WKBJ近似，根據(jù)模型修正量中豎直方向上的波數(shù)覆蓋，提出一種新的頻率點選取策略。結(jié)果表明，在接近實際的地下介質(zhì)速

3、度分布情況下，可以采用比Sirgue和Pratt方法更少的頻點來進(jìn)行頻域反演，而不會影響反演的精度，從而可以進(jìn)一步提高反演計算的效率。本文針對一維速度模型、Marmousi模型和Overthrust三個模型進(jìn)行反演實驗，模擬算例結(jié)果驗證了該優(yōu)化的頻率點選取方法的有效性。關(guān)鍵詞 頻率點選取 ；波數(shù)覆蓋 ； 線性遞增模型 ； 頻域波形反演AbstractThe method of frequency-domain full-waveform inversion is capable of obtaining high resolution results, but it suffered a he

4、avy computation when inverting large-scale models with all frequencies. Sirgue and Pratt (2004) proposed a discretization strategy that reduces the input frequencies required for inversion and increases efficiency without loss accuracy of the inversion results. At present, the most frequency selecti

5、on methods are based on the assumption that the velocity of underground medium follows uniform distribution, in order to study the situation that the background is more close to the actual underground medium, in this paper we use the velocity of medium linearly increases with depth, apply the Born a

6、pproximation of scattered field and the WKBJ approximation, according to the coverage of vertical direction wave number in model correction, and proposed a new frequency selection strategy. Experiment results showed that under the more realistic assumption, our method could obtain the same level of

7、accuracy with fewer frequencies than that required by Sirgue and Pratts strategy. In this paper we conducted experiments with 1D model, Marmousi model and Overthrust model, and the results demonstrated the validity of our method.Key words: Frequency selection; Wave-number coverage; Linear increasing

8、 velocity; Frequency domain waveform inversion1緒論隨著人們對地下非均勻介質(zhì)認(rèn)識的不斷深入以及計算機性能的大幅度提升，地震波反演成像技術(shù)也相應(yīng)地得到了非常迅猛地發(fā)展。然而，在使用傳統(tǒng)頻率點選取策略所得頻率點進(jìn)行頻率從低到高的迭代頻域波形反演的前提下，在勘探的目標(biāo)區(qū)域尺度非常大的情況下，計算效率就會成為此種反演方法所面臨的挑戰(zhàn)。此時，基于一種新的梯度圖像中波數(shù)覆蓋的頻率點選取策略可能會為地震波形反演帶來新的發(fā)展。本章將回顧相關(guān)的研究背景及現(xiàn)狀，簡述傳統(tǒng)的頻率點選取策略的原理，闡述它在實際應(yīng)用中的局限性及其改進(jìn)的可能性，給出本文的總體研究思路和主要研究

9、內(nèi)容。1.1研究背景和意義地震勘探(seismic exploration)是地球物理勘探中的一種重要方法，其主要目標(biāo)是使用人工震源激發(fā)的地震波來定位地下介質(zhì)中的礦藏(如石油，天然氣，煤，礦石，礦物，地?zé)崮茉?，水資源等)，獲得地下的巖石特性、碳?xì)浠衔飳傩缘鹊刭|(zhì)信息。在陸地地震勘探中，由地表或近地表的人工震源激發(fā)的地震波經(jīng)由地下介質(zhì)向地球內(nèi)部傳播，當(dāng)遇到地下介質(zhì)的物性參數(shù)快速變化區(qū)域時地震波便會產(chǎn)生散射波(包括反射波，折射波，繞射波等)，一部分返回到設(shè)置在地表的高靈敏度檢波器的散射波會被此檢波器記錄下來，所有的這些被不同檢波器(hydrophone)記錄的地震記錄構(gòu)成地震數(shù)據(jù)集；同理，在海洋地

10、震勘探中，返回到水面下方沿拖曳電纜方向設(shè)置的高靈敏度檢波器處的散射波會被此檢波器記錄下來，所有這些不同檢波器的散射波記錄構(gòu)成了地震數(shù)據(jù)集。對于上述兩種勘探情形下所得的地震數(shù)據(jù)集，其中主要記錄了地震波傳播的旅行時信息以及振幅響應(yīng)信息。通過分析地震數(shù)據(jù)集，人們一方面想建立一種從地震數(shù)據(jù)集到地下結(jié)構(gòu)概貌的映射，這一過程通常被稱作是地震成像或地震偏移成像；另一方面想通過地震數(shù)據(jù)集來對目標(biāo)區(qū)域中未知的物理參數(shù)(如速度，彈性模量，密度，阻抗等)進(jìn)行適當(dāng)估算，這一過程通常被稱作是地震反演；當(dāng)然，人們更希望可以通過分析地震數(shù)據(jù)集一次性地獲得目標(biāo)區(qū)域的結(jié)構(gòu)概貌和物性參數(shù)。由于地震數(shù)據(jù)集的復(fù)雜性，要想完全正確地實

11、現(xiàn)上述的想法依然是目標(biāo)我們所面臨的困難。地震數(shù)據(jù)集中包括了所有檢波器所能接收的散射波的能量信息，每一個能量信息代表了從震源到檢波器的地震波所特有的散射傳播方式，根據(jù)散射波在地下介質(zhì)中的具體傳播方式可以將其分為單次反射波和多次反射波。根據(jù)微擾理論，可以將地震波的單次散射波表示為目標(biāo)區(qū)域物性參數(shù)擾動的線性映射，相應(yīng)地多次反射則是一個非線性映射?？紤]到非線性反問題相對于線性問題的復(fù)雜性，目前大多數(shù)地震成像和地震反演算法都是基于數(shù)據(jù)中只包含地震波單次散射波的假設(shè)。由于地震數(shù)據(jù)在頻域中存在信息冗余(Mulder and Plessix 2004)，因此可以使用所有頻率點的一個子集來表示整體信息，考慮到頻

12、域中聲波方程的有限差分的計算量與頻率點數(shù)目成正比(Plessix et al. 2001)，利用此結(jié)論就可以將時間域中的波形反演轉(zhuǎn)化到頻率域中從而提高反演的計算效率。就從低頻率點到高頻率點的頻率迭代反演來說，其計算用時線性正比于波形反演時的頻率點數(shù)，因此如何高效地去除地震數(shù)據(jù)集中的冗余頻率點，也即如何合理地選擇最好(或最少)的頻率點組合，就成了提高此時計算效率的關(guān)鍵。現(xiàn)在常用的方法有，(1)基于均勻介質(zhì)假設(shè)的前提下采用固定間隔頻率點；(2)基于均勻介質(zhì)假設(shè)的前提下使用間隔隨頻率值線性遞增的頻率點；(3)通過速度值隨深度線性遞增模型中地震波的最大傳播深度以及多尺度反演中頻率點對局部極點的影響來對

13、方法(2)中的線性系數(shù)進(jìn)行簡單的修正。鑒于實際的地下介質(zhì)往往是非均勻的復(fù)雜情形，以及方法(3)只是對方法(2)的簡單修正，因此我們有必要研究一種新的頻率點選取策略。本論文旨在充分借鑒現(xiàn)有的頻率點選取策略，發(fā)現(xiàn)它們的不足之處，提出一種完全基于速度值隨深度遞增的線性速度模型下的新的頻率點選取方法，為下一步的研究打下良好的基礎(chǔ)。1.2國內(nèi)外發(fā)展現(xiàn)狀人們對于地震數(shù)據(jù)的處理方法一般分為地震偏移成像和地震反演，其中偏移成像的主要目的是通過恢復(fù)光滑背景模型上的不連續(xù)區(qū)域來構(gòu)建地下結(jié)構(gòu)概貌，而地震反演的目的則是定量地恢復(fù)地下介質(zhì)的物理參數(shù)。如果基于最小二乘法的失配函數(shù)(error functional)的成像

14、方法所得的黑塞(Hessian)矩陣是一個對角(或近對角)陣，那么此失配函數(shù)的最小值可以通過對梯度矩陣前乘黑塞矩陣(或偽(pseudo-)黑塞矩陣)的逆矩陣來獲得，此時的成像方法可以用于表示真振幅偏移成像(Mulder and Plessix 2004; Plessix and Mulder 2004)，考慮到基于最小二乘法的頻域波形反演的情形，我們可以發(fā)現(xiàn)地震成像和地震反演有時并沒有明確的界線(Tarantola 1986)。由于地震數(shù)據(jù)集在頻域表示時存在信息冗余(Mulder and Plessix 2004; Plessix et al. 2001)，因此頻域波形反演中可以使用全頻率點的

15、一個子集來提高反演的效率。Mulder和Plessix (2004)指出在頻率點選擇時需要考慮兩個主要準(zhǔn)則：(1)由等價海森堡測不準(zhǔn)原理(Heisenbergs principle)所確定的偏移圖像中的反射體分辨率；(2)頻域奈奎斯特抽樣定理(Nyquists theorem)，考慮到偏移包括在給定背景模型上的地震波傳播以及多個震源檢波器的求和，因此盡管有時不符合奈奎斯特抽樣定理但依然可以獲得一個沒有混疊的結(jié)果。在均勻背景的假設(shè)下，利用上述的兩準(zhǔn)則，得到了一個頻率點間隔為固定值的頻率點選取策略。Sirgue和Pratt (2004)根據(jù)頻域最速下降反演方法，在均勻介質(zhì)假設(shè)的情況下，利用偏移成像

16、和頻域波形反演的等價性(Tarantola 1986)，借助平面波近似及波恩近似(Born approximation)原理，揭示了利用多震源檢波距地震數(shù)據(jù)進(jìn)行頻域反演時在垂直方向上波數(shù)覆蓋存在冗余的現(xiàn)象，據(jù)此提出了一種適用此情形下的頻率點間隔隨頻率線性遞增的選取策略。同年，Yokota和Matsushima (2004)在均勻介質(zhì)假設(shè)基礎(chǔ)上，根據(jù)觀測數(shù)據(jù)與正演數(shù)據(jù)之間的相位誤差大于時就會在結(jié)果上出現(xiàn)跳周現(xiàn)象(Cycle Skipping，即一個異常高的渡越時間記錄)的情形，進(jìn)而提出一種相對應(yīng)的頻率點選取策略。Mulder和Plessix (2008)利用速度值隨深度線性遞增模型中地震波傳播的

17、最大深度以及多尺度反演中頻率點對局部極點的影響(Bunks et al. 1995)，對Sirgue和Pratt (2004)所提出的頻率點選取策略中的線性系數(shù)進(jìn)行了簡單修正。Anagaw和Sacchi (2014)以及Kim等 (2011)對Sirgue和Pratt (2004)的方法所得的頻率點進(jìn)行簡單地分組，從而驗證不同的頻率點組合對反演的影響。目前來說，對于頻率點選取策略的研究還不是很多，前面所述的方法基本上都是基于均勻背景(或只是對均勻背景簡單修正)的情形。隨著反演理論的發(fā)展以及勘探的需要，人們自然要問是否存在基于相對于均勻背景而言更一般的頻率點選取策略。目前，該課題的研究還處于初級

18、階段，對于從事反演研究的學(xué)者來說它是一個具有一定挑戰(zhàn)性的研究方向。1.3主要研究內(nèi)容在認(rèn)真分析總結(jié)現(xiàn)有頻率點選取策略的基于基礎(chǔ)上，本文結(jié)合傅里葉變換的意義，首先建立基于一階Born近似的散射波場積分表達(dá)式；然后在WKBJ近似的情形下獲得Green函數(shù)的解析表達(dá)式；在此基礎(chǔ)上結(jié)合頻域波形反演算法通過梯度圖像中的波數(shù)覆蓋情況推導(dǎo)出線性遞增速度背景模型下相鄰兩個頻率點之間的關(guān)系式，進(jìn)而得到新的頻率點選取策略。本論文各章節(jié)的內(nèi)容安排如下：第一章 介紹了本論文的研究目的及意義，對相關(guān)研究現(xiàn)狀進(jìn)行總結(jié)，并簡要介紹了本論文的研究內(nèi)容及安排。第二章 回顧了頻域波形反演的基于原理和方法，并對預(yù)處理梯度法進(jìn)行了簡

19、要介紹。第三章 針對線性遞增速度背景模型，提出一種基于散射波場Born近似及WKBJ近似的新的頻率點選取策略，并在不同源檢距組合情形下，推導(dǎo)了反演中相鄰兩個頻率點之間的關(guān)系。第四章 對于本文提出的新的頻率點選取策略，在不同的背景速度模型下與Sirgue和Pratt (2004)提出的方法進(jìn)行模型測試。詳細(xì)對比兩種頻率點選取策略下反演結(jié)果的分辨率及計算用時。第五章 總結(jié)本論文的相關(guān)研究結(jié)果，分析研究中存在的問題，并展望下一步研究的工作重點。2頻域波形反演算法簡介地震波形反演是一個非線性問題，為了降低問題的求解難度，通常會使用基于梯度的線性處理方法來對其進(jìn)行求解。常用了求解方法有：全牛頓法，高斯牛

20、頓法，最速下降法，預(yù)處理梯度法等。在實際使用中，由于黑塞矩陣的尺寸及其較高的計算消耗，一般很難對其進(jìn)行直接計算(Mulder and Plessix 2004; Pratt et al. 1998; Tarantola 1984)。通常，最速下降法具有較低的收斂速度(Tarantola 1986)，考慮到黑塞矩陣的求解難度，我們可以通過使用偽黑塞矩陣預(yù)處理來加速梯度法的收斂過程(Mulder and Plessix 2004; Plessix and Mulder 2004)。本章以頻域反演為前提，主要介紹預(yù)處理梯度反演方法，為后續(xù)的研究奠定基礎(chǔ)。2.1 波形反演的基本步驟目前地震成像中絕大多

21、數(shù)的偏移成像方法都是采用線性化簡化，其思想為對散射波場進(jìn)行線性化處理(也即Born近似)，然后將Born近似的積分形式自然地與廣義的Radon變換結(jié)合起來，進(jìn)而通過逆廣義Radon變換來求取近似地反演解，進(jìn)而恢復(fù)地下介質(zhì)結(jié)構(gòu)。就其中的數(shù)學(xué)部分而言，Radon已于1917年首次證明了利用投影函數(shù)可以恢復(fù)出滿足條件的唯一的圖像函數(shù)，因此利用上述的偏移成像方法可以唯一地重建地下介質(zhì)結(jié)構(gòu)，以達(dá)到地震勘探的最終目的。而我們這里要介紹的波形反演方法也遵從上面所述的理論，它們都是通過建立勘探所得的實際數(shù)據(jù)與用理論計算得到的模擬數(shù)據(jù)的數(shù)據(jù)殘差，并將此殘差與實際地下模型與先驗地下模型之間的擾動之間形成映射，通過

22、不斷修正先驗地下模型來達(dá)到最終獲得一個在其上的模擬數(shù)據(jù)與勘探數(shù)據(jù)擬合度最好的地下模型(Anagaw and Sacchi 2014; Gauthier et al. 1986; Mora 1987; Pratt 1999; Tarantola 1984)。原理雖然相同，但是對于具體的實現(xiàn)方法之間則有較大的區(qū)分，為了論文的需要這里僅就波形反演的方法進(jìn)行討論，其實現(xiàn)過程基本可分為以下幾個部分：(1)在某一地理區(qū)域設(shè)置勘探的震源及相應(yīng)的檢波器(也即接收點)。(2)采集各個震源所對應(yīng)的檢波器上的波形記錄，形成實測(觀測)地震記錄。由于實測地震記錄是波形反演中的表征地下結(jié)構(gòu)的最重要的數(shù)據(jù)，故實測地震記錄

23、的質(zhì)量顯著地影響了波形反演的效果及可信程度。(3)根據(jù)實測地震記錄中震源及檢波器的坐標(biāo)位置以及檢波器中的地震道數(shù)據(jù)來建立初始介質(zhì)模型，也即先驗?zāi)Ｐ?。在通常的基于梯度的局部最?yōu)化反演方法中，如果先驗?zāi)Ｐ团c實際勘探區(qū)域的實際模型差別過大，則在反演過程中陷入局部極小點的可能性大大增加，在實際的應(yīng)用過程中一般可以使用走時反演的反演結(jié)果作為初始模型，而對于基于全局優(yōu)化的方法則可以減弱對初始模型的先驗限制。(4)在初始介質(zhì)模型上按實際勘探時的震源及檢波器的設(shè)置情形進(jìn)行其上的正演模擬，并記錄檢波器處的波形信息，即預(yù)測(正演)地震記錄。(5)把正演地震記錄與觀測地震記錄按每個檢波器做差得到地震記錄殘差，由波恩

24、近似可以看出地震記錄殘差為先驗?zāi)Ｐ团c目標(biāo)模型之間的模型擾動的映射，由Radon變換可知可用此地震記錄殘差來解得模型的擾動也即恢復(fù)地下介質(zhì)結(jié)構(gòu)。將地震記錄殘差表示成列矩陣的形式，其中的每個元素代表相應(yīng)的檢波器上的地震記錄的殘差，把地震殘差矩陣按一定方法映射到一個標(biāo)量目標(biāo)函數(shù)(如F范數(shù))。(6)在目標(biāo)函數(shù)的形式上，使用泛函分析的相關(guān)理論獲得目標(biāo)函數(shù)極小時所對應(yīng)的相對于初始介質(zhì)模型的模型修正量，用此修正量修正初始模型。把修正后的初始模型當(dāng)作先驗介質(zhì)模型重復(fù)步驟(2)-(6)，直到修正后的初始模型上的目標(biāo)函數(shù)滿足設(shè)定的閾值時停止對介質(zhì)模型的迭代更新。概括起來波形反演的基本步驟可分為：獲取觀測地震記錄，

25、設(shè)置初始介質(zhì)模型，在初始介質(zhì)模型上正演獲得正演地震記錄，對正演地震記錄與觀測地震記錄做差并將結(jié)果表示成波場殘差矩陣，然后將波場殘差矩陣映射成標(biāo)量目標(biāo)函數(shù)，再通過泛函方法獲得介質(zhì)模型修正量并按此對介質(zhì)模型進(jìn)行更新，具體的流程如圖2.1所示。圖2.1 波形反演流程圖2.2 頻域正演模擬地震波的速度是指地震波在巖層中的傳播速度，簡稱地震速度，有時又叫做巖石速度。地震速度是地震勘探中最重要的一個參數(shù)，也是一個復(fù)雜的參數(shù)，其不僅與巖石的物理性質(zhì)(包括：巖性、孔隙度、孔隙中的充填物、巖石的密度、埋藏深度和地質(zhì)年代等)有關(guān)，而且還與介質(zhì)的結(jié)構(gòu)和求取速度的方法緊密相關(guān)。本文中我們選用基于地震速度的函數(shù)作為介質(zhì)

26、模型的參數(shù)。2.2.1 二維頻域波動方程非均勻介質(zhì)中的地震波的傳播是一個十分復(fù)雜的過程，在實際的傳播過程中一般可以將地震波分為縱波和橫波兩種形式，在傳播過程中縱波的傳播速度大概是橫波的二倍左右。考慮到在實際野外勘探記錄中的記錄時長一般不太長，因此我們完全可以在數(shù)據(jù)集中濾掉橫波然后用只能描述縱波的聲波方程來對地震波縱波的傳播過程進(jìn)行模擬，這樣可以提高波形正反演的精度。這里我們只考慮地震波在二維恒定密度速度介質(zhì)中傳播的情形，其具體的聲波方程可以表示如下：, (2.1)其中x表示水平方向的坐標(biāo)，z表示豎直方向的坐標(biāo)，t為時間，s(x, z, t)表示在坐標(biāo)(x, z)處的時域震源，p(x, z, t

27、)表示由震源s(x, z, t)激發(fā)的地震波在坐標(biāo)(x, z)處的接收點所接收到的時域波形記錄。表示拉普拉斯算子，其具體形式可表示如下. (2.2)對式(2.1)的兩邊同時進(jìn)行時域傅里葉變換，也即把聲波方程從時間域變換到頻率域，可得, (2.3)其中為角頻率，P(x, z, )為時域接收點波場p(x, z, t)的傅里葉變換，表示頻域接收點波場，S(x, z, )為時域震源s(x, z, t)的傅里葉變換，表示頻域震源。在對式(2.3)所表示的微分方程進(jìn)行有限差分求解時，可以看到拉普拉斯算子離散時所使用到的所有坐標(biāo)點構(gòu)成一個十字線，也即只使用了左、右、上、下四個方向上的離散點，為了在有限差分的

28、具體實現(xiàn)中包含左上、左下、右上、右下方向上的離散點，建立一個兩條坐標(biāo)軸分別為離散網(wǎng)格不同對角線的新的坐標(biāo)系(仿射坐標(biāo)系)，其具體情形可由圖2.2給出。若記原來XZ坐標(biāo)系的x軸正方向單位矢量為i，z軸正方向單位矢量為j, 在新坐標(biāo)系的軸正方向單位矢量為，軸正方向單位適量為。并記離散的網(wǎng)格的水平間距為，豎直間距為，對角線間距為，可看出此三者滿足下式. (2.4)從圖2.2可以看出新舊兩個坐標(biāo)系軸向上的單位向量之間應(yīng)滿足如下關(guān)系式, (2.5). (2.6)假設(shè)一個向量在原坐標(biāo)系中表示為(x, z)，而在新坐標(biāo)系中則相應(yīng)地表示為，由于這兩個表示均代表同一個向量，故對于同一向量而言其在新舊坐標(biāo)系下的坐

29、標(biāo)之間應(yīng)滿足如下關(guān)系：, (2.7). (2.8)由式(2.7)和式(2.8)可知，在原坐標(biāo)系下的拉普拉斯算子在新坐標(biāo)系中可以表示為如下形式. (2.9)由式(2.9)可以看出拉普拉斯算子在新坐標(biāo)系中的形式相對于式(2.2)的形式而言，出現(xiàn)了混合偏導(dǎo)形式，為了使之具有與式(2.2)一樣的簡潔的形式，我們在此人為地令網(wǎng)格的水平間距和豎直間距相同，并記, (2.10)則式(2.9)可以簡潔地表示為如下形式. (2.11)由式(2.11)可以看出其與式(2.2)具有非常類似的表達(dá)形式，而相應(yīng)地如果用新坐標(biāo)來表示聲波方程在地下速度介質(zhì)中的傳播規(guī)律則可將聲波方程表示為如下形式：, (2.12)其中為拉普

30、拉斯算子在新坐標(biāo)系下的表示，也即. (2.13)由前面討論可知，式(2.3)和式(2.12)分別為頻域聲波方程在原坐標(biāo)系及新坐標(biāo)系下的表達(dá)式，若和表示同一個接收點的波形記錄，也即坐標(biāo)(x, z)和坐標(biāo)滿足式(2.7)和式(2.8)，那么我們可以將式(2.3)和式(2.12)通過加權(quán)的方式將其合并為下式. (2.14)其中表示在合成的表達(dá)式中新舊坐標(biāo)系下聲波方程所占的比例，在實際應(yīng)用中，可以通過最優(yōu)化方法來求得此參量的使用值。在后文中可以看到，這時在聲波方程的有限差分中不再只是使用十字線上的點同時使用離散網(wǎng)格對角線上的點(而相應(yīng)的不同位置的點所點的權(quán)重則由反映)的差分式，并且當(dāng)時則恢復(fù)到式所表示

31、的情形。圖2.2 兩種坐標(biāo)系的表示2.2.2 有限差分表示有限差分法是使用差分原理來解決微分方程的一種常用的數(shù)值求解方法，其中心思想為：將需要進(jìn)行數(shù)值計算的區(qū)域離散成很多的點的集合，離散方法通常為離散為正方形或矩形網(wǎng)格(此時后續(xù)的求解在形式上較為簡潔)，然后對方程進(jìn)行微分或偏微分向差商形式的代換，在最終的差分表達(dá)形式相對于原來的微分或偏微分方程而言是相容及收斂的，并且差分形式本身是穩(wěn)定的時候那么就可以使用此差分形式來對原微分或偏微分方程進(jìn)行數(shù)值求解。對于頻域的聲波方程而言通常有多種離散方法(Gu et al. 2013; Jo et al. 1996; Liu 2014)，本文中我們選用一種新

32、型的9點離散策略。首先把速度模型按水平及豎直方向間隔均為的正方形網(wǎng)格進(jìn)行離散，并把此對角線的長度記為(明顯有)，離散后水平方向的總的網(wǎng)格點數(shù)目記為，相應(yīng)的豎直方向的總的網(wǎng)格點數(shù)目記為。為了用有限差分法求解式(2.3)、式(2.12)及式(2.14)所表示的聲波方程，需要使用泰勒級數(shù)來完成微分形式到差分形式的轉(zhuǎn)換，對于一元函數(shù)由高數(shù)知識可知若函數(shù)在的某個領(lǐng)域中具有直到(n+1)階的導(dǎo)數(shù)，則在該領(lǐng)域內(nèi)的n階泰勒公式為：. (2.15)其中為拉格朗日余項，在不影響計算或影響很小的情況下可以將其直接計為0。為了求解的方便，我們把簡記為，同理把簡記為。把波場及波場按式(2.15)的形式表示為關(guān)于坐標(biāo)點(

33、m, n)的函數(shù)，可得：, (2.16). (2.17)把式(2.16)和式(2.17)等號兩端分別相加，再經(jīng)過簡單的化簡可得：, (2.18)同理，可得：. (2.19)把式(2.18)與式(2.19)在等號兩邊分別求和可得：. (2.20)同理可求得在新坐標(biāo)系下拉普拉斯算子的相應(yīng)的表達(dá)為：. (2.21)由式(2.7)及式(2.8)，可將式(2.21)中的在新坐標(biāo)系下表示的四階偏導(dǎo)轉(zhuǎn)化為在原來坐標(biāo)系下的情形，可得：. (2.22)如果和表示同一個接收點的數(shù)據(jù)，則可利用式(2.20)和式(2.22)把在混合新舊兩種坐標(biāo)系的波動方程式(2.14)離散為如下所示形式：, (2.33)其中的R為式

34、(2.14)離散時的截斷誤差，其具體形式可根據(jù)前面的表達(dá)表示如下：. (2.34)由圖2.3中離散點的分布情形，以及離散時的離散網(wǎng)格為正方形，則式(2.33)及式(2.34)可化為, (2.35). (2.36)由式(2.36)可知，當(dāng)時，R的值也趨于零，故可忽略式(2.35)中的R進(jìn)而將式(2.35)表示為：. (2.37)為了減少聲波方程離散化的頻散，進(jìn)一步提升正演的精度，可對式(2.37)中所表示的質(zhì)量加速項進(jìn)行優(yōu)化(Jo et al. 1996)，具體的優(yōu)化可以有多種，這里選用上、下、左、右及其本身五個位置上的波場值在一階近似的情況下來近似地表示，也即：, (2.38)其中的參量表示在

35、近似時中心點所占的比例，在實際應(yīng)用中，可以通過最優(yōu)化方法來求得此參量的使用值。把式(2.38)代入式(2.37)可得：, (2.39)顯然上式相對于式(2.14)而言是相容且收斂的，當(dāng)離散網(wǎng)格間距、模型的速度值以及相應(yīng)的頻率點之間滿足一定關(guān)系時(簡單情形下一個波長覆蓋4個或者更多的網(wǎng)格點時)，式(2.39)所表示的離散方程本身是穩(wěn)定的(Boonyasiriwat et al. 2009; Chen 2014; Chen and Cao 2014; Huang and Schuster 2014; Jo et al. 1996; 劉璐等 2013)，故可以將式(2.39)作為聲波方程的有限差分形

36、式來對其進(jìn)行數(shù)值求解。對于離散后的每一個網(wǎng)格點而言，均存在一個如式(2.39)所示的線性方程，就整體而言可形成一個由個方程組成的線性方程組。這樣在速度模型上的正演模擬過程可以表示為如下的線性方程組形式：, (2.40)其中為表示模型空間的壓強場的列向量，速度模型中沿X軸正方向變動過程中沿Z軸正方向從上到下的網(wǎng)格點分別與中從上到下的元素一一對應(yīng)；為離散尺度、頻率點和差分格式等有關(guān)的復(fù)數(shù)方陣，一般情況下，矩陣是高度稀疏、非對稱及非正定的；為震源向量。通過對(2.40)進(jìn)行求解可以獲得此時的壓強場的列向量，也即正演結(jié)果，而求解方法一般分為兩類：迭代法和直接求解法(Pratt 1999)，本文中則通過

37、調(diào)用mumps庫來進(jìn)行直接求解。對于解出的結(jié)果，由于我們只對設(shè)置了檢波器的位置上的波場值感興趣，故我們需要對列矩陣進(jìn)行處理，即保留檢波器位置所對應(yīng)的元素值，而把其余的所未設(shè)置檢波器位置的網(wǎng)格點所對應(yīng)的元素值設(shè)置為0，經(jīng)此處理后的列矩陣記錄了給定頻率點上的波場記錄，對于所有頻率點的波場數(shù)據(jù)進(jìn)行反傅里葉變換即可獲得給定先驗速度模型上的正演波場記錄。圖2.3 在不同坐標(biāo)系中的離散網(wǎng)格點2.2.3 PML層邊界處理方法對于真實的地下介質(zhì)模型而言，一般而言是一個無限區(qū)域，但是由于使用計算機計算的有限差分法只能對有限的計算的計算量進(jìn)行計算，因此在離散時只能考慮我們感興趣的區(qū)域，這樣由于人為邊界所產(chǎn)生邊界反

38、射就不可避免有限差分所得的結(jié)果與與真實情況之間出現(xiàn)差異，就離散后的模型而言越靠近邊界區(qū)域這種差異表現(xiàn)的越明顯。為了盡量地消除邊界反射，在這里我們采用一種常用的也比較有代表性的方法完全匹配層邊界吸收法，這種方法最開始主要用于電磁波傳播規(guī)律的有限差分方程的邊界處理中，后來逐漸被用在地震層析成像中，相對于其它的邊界處理方法而言，PML方法(Berenger 1994; Hustedt et al. 2003; Liu et al. 2011; 郭念民 和 吳國忱 2012; 劉璐等 2013; 張毅 2007)是目前來說最有效并且相對簡潔的方法，這一點在相關(guān)的有很多實驗中已經(jīng)得到了驗證。由于PML層

39、是對有限差分法離散后的區(qū)域的邊界之外的區(qū)域進(jìn)行處理，因此可知此區(qū)域中并沒有任何的震源，故此時的時域地震波方程相對于式(2.1)而言只是使后者的震源變?yōu)榱阒?，也即? (2.41)為了使上式降階，也即由兩個一階偏微分方程來表示PML層中地震波傳播規(guī)律，在這里我們引入兩個中間變量和，中間變量中的下標(biāo)x和z只是一個記號標(biāo)記并不指代任何其它的方向信息，且這兩個中間變量滿足以下方程：, (2.42); (2.43)對于時域的波場響應(yīng)而言，我們在這里將其拆分成和兩部分，相應(yīng)分量的下標(biāo)x和z只是一個記號標(biāo)記并不指代任何其它的方向信息，也即：, (2.44)由數(shù)學(xué)關(guān)系可知一定存在上述的波場響應(yīng)的拆分方式，使得

40、下兩式成立，也即：, (2.45). (2.46)式(2.45)及式(2.46)即為我們用來表示PML層中聲波方程的兩個一階偏微分方程，由形式上我們可以看出上兩式所描述的地震波在傳播過程中并沒有能量的衰減，也即能量是守恒的，這一點使得相對于沒有邊界處理時的情形而言只是相當(dāng)于擴大離散區(qū)域的大小，當(dāng)然如果我們能把區(qū)域擴展到無窮大那么在物理上來說也是對真實問題的較好的近似，但是前面我們也提過這一點對于實際計算而言是沒有任何意義的，因此式(2.45)和式(2.46)的邊界處理對于邊界處的地震波反射而言并沒有得到根本地解決。因此為了消除邊界處的波形反射，PML區(qū)域中傳播的地震波必須在能量上出現(xiàn)一定程度上

41、的衰減，根據(jù)一階常微分方程的解的一般規(guī)律，在這里我們對式(2.42)、式(2.43)、式(2.45)和式(2.46)所表示的PML區(qū)域中的地震波傳播的規(guī)律中加入衰減項(劉璐等 2013)，也即：, (2.47), (2.48), (2.49), (2.50)其中和分別表示在坐標(biāo)(x, z)處下標(biāo)標(biāo)記為x和z的波場的衰減系數(shù)，在實際應(yīng)用中有多種方法來對不同空間坐標(biāo)點的衰減系數(shù)進(jìn)行選取，本文中所使用的衰減系數(shù)規(guī)律為，其中為最大的衰減系數(shù)， L為所添加的PML層的厚度，d為坐標(biāo)點離已經(jīng)離散的區(qū)域的邊界的距離，對于不同的下標(biāo)標(biāo)記的衰減系數(shù)其相應(yīng)的和L可以選用不同的值，本文中的值我們?nèi)?0。對式(2.

42、47)、式(2.48)、式(2.49)及式(2.50)進(jìn)行傅里葉變換將它們從時間域變換到頻域可得：, (2.51), (2.52), (2.53), (2.54)其中, (2.55). (2.56)對于新坐標(biāo)系下的情況而言可以得到類似于式(2.51)、式(2.52)、式(2.53)及式(2.54)所表示的地震波方程，對所得到的新舊坐標(biāo)系下的方程進(jìn)行有限差分處理，并進(jìn)行相應(yīng)的簡化可得類似于式(2.37)的形式，也即：. (2.57)式(2.38)代回式(2.57)可得類似于式(2.39)的形式，也即：. (2.58)如果上式中的所有衰減系數(shù)就限定為0，那么式(2.58)可化為式(2.39)中的所

43、有震源都為零的情形。如果我們把計算分為普通區(qū)域和PML層區(qū)域，并在不同的區(qū)域中分別使用式(2.39)和式(2.58)所表示的地震波傳播方程來進(jìn)行相應(yīng)的正演計算，那么在形式上而在無形中就給計算增加了麻煩，這里為了計算的簡單及形式上的簡潔，我們把普通區(qū)域和PML層區(qū)域整合成一個區(qū)域計算區(qū)域，結(jié)合式(2.39)和式(2.58)的形式，我們可以得到描述此的方程可以表示為：. (2.59)對于普通區(qū)域及PML層區(qū)域中的每個離散網(wǎng)格點都有一個與式(2.59)相對應(yīng)的式子，它們的總體可以組成一個方程的形式，可以用式(2.40)的形式來表示這個矩陣方程。對于添加的PML層而言，一般除了表層上層區(qū)域外都會添加不

44、同厚度的PML層，而對于表層上層的區(qū)域而言，如果是對陸地進(jìn)行勘探則不需要添加任何PML層，但如果是對海底進(jìn)行勘探則需要添加相應(yīng)厚度的PML層。而對于實際的介質(zhì)地表而言，一般都是不規(guī)則或有較多起伏的，這時就需要有相應(yīng)的方法對起伏的地表進(jìn)行相應(yīng)地處理來使之用于正反演過程，本文由于使用的是平地表模型，故而在此不對此問題加以論述。2.3 頻域波形反演上文中我們已經(jīng)介紹了給定區(qū)域中地震波的正演方法，本節(jié)中則主要對本文中將要使用的頻域波形反演進(jìn)行簡要介紹。由正演的差分形式可知，由于正演數(shù)據(jù)相對于介質(zhì)參數(shù)(速度值)而言并不是一個簡單的線性關(guān)系，因此波形反演問題是一個比較復(fù)雜的非線性問題。反演問題的最終目的是

45、使我們所得的最終模型與真實的地下介質(zhì)模型相同或非常接近，但由于我們并不知道真實地下結(jié)構(gòu)(否則就不需要反演了)，而僅僅知道在相應(yīng)地理區(qū)域中由給定位置震源在相應(yīng)檢波器點所激發(fā)的波形響應(yīng)，也即地震數(shù)據(jù)集，因此我們通過實際數(shù)據(jù)集與正演數(shù)據(jù)集的差異目標(biāo)函數(shù)來衡量反演結(jié)果的優(yōu)劣，目標(biāo)函數(shù)值越小則認(rèn)為反演結(jié)果越好。首先對正演模擬波場及勘探所得波場之間做差，并把結(jié)果記為波場殘差，也即：, (2.59)其中的P和均為列矩陣，其中從上到下的每個元素都與離散后并添加了PML層后的網(wǎng)格點上的波場值一一對應(yīng)，但由于勘探波場只是設(shè)置了檢波器的網(wǎng)格點上的波場響應(yīng)，故對于沒有設(shè)置檢波器的網(wǎng)格點所對應(yīng)的元素值我們將其設(shè)置為0。

46、我們可以使用波場殘差的F范數(shù)來表示相應(yīng)的目標(biāo)函數(shù)(Brossier et al. 2010; Shin et al. 2007)，也即：. (2.60)其中M表示模型參數(shù)，在聲波方程中一般表示聲波網(wǎng)格點聲波速度值或聲波慢度或者為前兩者的函數(shù)。由前面的討論可知，反演的過程可以等價于式(2.60)所表示的目標(biāo)函數(shù)收斂到全局極點值的過程。假設(shè)先驗?zāi)Ｐ臀挥谌謽O小點附近，在此我們可以使用基于梯度的方法來完成對原模型的修正，首先應(yīng)先求得目標(biāo)函數(shù)對模型的整體參數(shù)的梯度，也即：, (2.61)其中Re表示取其后得數(shù)的實部，上標(biāo)T表示矩陣的轉(zhuǎn)置，J表示雅可比矩陣，其具體形式如下：, (2.62)可以看出J是二

47、維矩陣，其第i行第j列的元素值可以表示為. (2.62)如果使用式(2.61)式來直接求取梯度矩陣，由于當(dāng)模型規(guī)模增大時求取雅可比矩陣的時耗顯著變大，因此不能有效地用于大規(guī)模模型情形，為了解決這個問題，這里首先對式(2.40)所表示的正演方程兩邊同時對模型參數(shù)M求偏導(dǎo)，考慮到其中震源函數(shù)本身與模型參數(shù)無關(guān)，可得, (2.63)其中A-1表示矩陣A的逆矩陣，可以看出上式實際上是雅可比矩陣的另一種表示形式，將其代回式(2.61)并考慮到矩陣A的近似對稱性(Pratt et al. 1998)，可得. (2.64)其中的上標(biāo)*表示矩陣的共軛。可以看出目標(biāo)函數(shù)關(guān)于模型參數(shù)的梯度可以表示為正演波場與回傳

48、波場的零偏移相關(guān)函數(shù)表示(Gauthier et al. 1986; Mora 1987; Pratt et al. 1998; Ravaut et al. 2004; Tarantola 1984)。目標(biāo)函數(shù)可以在模型參數(shù)M處按泰勒級數(shù)進(jìn)行展開，我們忽略級數(shù)中二階以上的高階小項可得下式：, (2.65)其中H為黑塞矩陣，其具體形式為, (2.66)如果模型參數(shù)M表示真實地下介質(zhì)模型，那么若對式(2.65)兩端同時對求一階導(dǎo)數(shù)，也即, (2.67)由于全局極小點存在且為M，式(2.67)的左端值為0，故可得. (2.68)上式即為所得的實際修正量，用此式所表示的波形反演方法通常被稱作是牛頓法(

49、或全牛頓法)(Pratt et al. 1998)。相關(guān)的大量實驗表明式(2.66)所示的矩陣其值主要集中在第一項上，我們可以把此部分記為Ha，也即, (2.67)其中E表示單位矩陣，上標(biāo)H表示矩陣的共軛轉(zhuǎn)置，此時可稱Ha為近似黑塞矩陣(Pratt et al. 1998)，將其代回式(2.68)時所得到的即為高斯牛頓法的表達(dá)式。若我們用式(2.67)所代表矩陣的對角陣來近似地表示完整的黑塞矩陣(此時的矩陣稱為偽黑塞矩陣)(Plessix and Mulder 2004)，則模型參數(shù)的修正量可以表示為(Operto et al. 2006; Ravaut et al. 2004; Sourbi

50、er et al. 2009), (2.68)其中表示修正步長，為正則化因子，diag表示矩陣對角線上元素保持不變而其它位置元素置零?？梢允褂命S金分割搜索法、拋物線法、二分搜索法等多種進(jìn)行求解，本文中采用拋物線法?？梢钥醋魇菍μ荻鹊念A(yù)處理。式(2.68)表示的是對于單震源多檢波器點的情形，對于通常使用的多震源多檢波器的情形可以通過對單震源的情形的加權(quán)求和來確定其使用的模型參數(shù)的修正量。式(2.68)所表示的預(yù)處理梯度法，雖然仍然需要計算黑塞矩陣，但由于只是計算對角元上的元素，因此相對于式(2.68)所表示的牛頓法而言具體較小的計算費用，由于對梯度矩陣進(jìn)行了預(yù)處理，故相對于通常的最速下降法而言有

51、較好的收斂效率。由式(2.64)所表示的梯度形式及式(2.67)所表示的近似黑塞矩陣的形式可以看出，預(yù)處理梯度法的反演過程可以表示成多次的正演過程，考慮到反演時設(shè)置的震源與實際情形下的震源通常是不同的，故應(yīng)對震源進(jìn)行源估計(Koo et al. 2011; Pratt 1999; Operto et al. 2006; Ravaut et al. 2004)。3頻率點的選取策略基于均勻背景模型，根據(jù)不同的原理或假設(shè)已經(jīng)提出了多種頻率點選取策略(或只是使用其它背景模型下的一些結(jié)果對所得策略做簡單修正)，但就真實地下非均勻介質(zhì)結(jié)構(gòu)而言，此模型假設(shè)存在較大的區(qū)別。為了更好地用于真實地下介質(zhì)反演，本章

52、基于線性遞增速度背景模型，根據(jù)散射波場Born近似及WKBJ近似及反演時梯度圖像上豎直方向上的波數(shù)覆蓋情形，提出一種新的頻率點選取策略。3.1 散射波場的Born近似及Green函數(shù)求解在前面波形反演的算法介紹中，我們只是在二維情形中對正演模擬及反演算法的具體實現(xiàn)進(jìn)行了相關(guān)討論，為了問題的一般性，我們將在n維歐幾里得空間中討論散射波場的線性化(Born近似)和格林函數(shù)(脈沖點震源在空間的響應(yīng))的求解問題。對于空間中的任意點可用笛卡爾坐標(biāo)向量表示為，當(dāng)n=3時就是我們所熟知的三維空間，當(dāng)n=2時就是前面討論中所使用的二維面空間。不同于前面的討論中的面震源，這里我們假設(shè)震源為點源(面震源可以看作多

53、個點震源的疊加情形)，對于震源在處的點源，在速度模型中由此點源產(chǎn)生的地震波會被一個或一組在處的檢波器檢測。反演中速度模型所在的區(qū)域的全體記作D，把區(qū)域的D邊界標(biāo)記為，若點源的波形函數(shù)在頻域中可表示為，故可以使用克羅內(nèi)克函數(shù)把震源的空間及波形信息表示為的形式。此時，常數(shù)密度介質(zhì)中的標(biāo)量聲波方程可以以類似式(2.3)的形式表示為：. (3.1)對于速度模型來說，其模型參量是每個位置上的速度值，若把其上速度值的倒數(shù)記為，通常稱作介質(zhì)在此點的折射率(可以用此來代表整個模型)，也即. (3.2)根據(jù)微擾原理，速度模型中模型參量可以表示成一個已知的背景速度模型及在其上的一個微小擾動。為了計算的簡便性，這里

54、將其如下形式表示(Beylkin et al. 1985)：, (3.3)其中背景模型不一定是均勻不變的，是一個未知的微小擾動。由上式可看出實際上反映了預(yù)測速度模型與真實地下結(jié)構(gòu)之間的差異程度，如果能準(zhǔn)確的求出此值，那么就能完全的還原地下介質(zhì)的具體速度分布。為了使Sommerfeld輻射條件有效(Colon 1988)，應(yīng)限制在全空間的一個小子集中。將式(3.3)代入式(3.1)可得：, (3.4)其中可看作由式子右端的原脈沖源和形式比較復(fù)雜的散射源聯(lián)合作用下的波場，故可將其記全波場。上式顯然是一個Helmholz方程，而其Helmholz操作算子則為。就形式上可以看出，散射源是一個面源而脈沖

55、源只是一個點源，顯然前者含有更多的模型信息特別是速度模型深層區(qū)域的參量信息，在反演中人們非常重視此部分信息?？紤]到脈沖源和散射源在能量強度上有明顯區(qū)別，參考式(3.3)中模型參數(shù)的拆分方式，我們可以把全波場由參考波場和微擾波場兩部分組成(Beylkin et al. 1985; 陳生昌等 2001)，其中參考波場對應(yīng)的震源為脈沖源，而微擾波場所對應(yīng)的震源為散射源。上述的對全波場的劃分及劃分后每部分所對應(yīng)的條件可以表示為：, (3.5), (3.6). (3.7)為了求出式(3.7)中的微擾波場的顯式表達(dá)形式，引入格林函數(shù)和格林函數(shù)，其分別對應(yīng)在源點及檢波器點設(shè)置的脈沖點源在空間中產(chǎn)生的響應(yīng)波形

56、，也即：, (3.8). (3.9)為了求解微擾波場，可構(gòu)建如下方程：. (3.10)由Sommerfeld輻射條件可知，上式的等式右端的面積積分項的值為0，在此基礎(chǔ)上利用式(3.7)式(3.9)可得. (3.11)由于上式中等式右邊中含有微擾與散射波場的乘積項，因此式(3.11)是非線性的，通常我們會把難以直接求解的非線性問題通常一定的近似或其它方法將問題近似地轉(zhuǎn)換我們比較容易求解的線性問題，然后把相應(yīng)的解當(dāng)作是滿足一次近似情形下的非線性問題的解，然后再通過其它方法對此近似解進(jìn)行必要的校正。若能對式(3.11)進(jìn)行線性化處理，也即能忽略式中的項，這意味著需要滿足一定的先驗條件，例如當(dāng)和在數(shù)值

57、上基本相當(dāng)時顯然不能忽略上述項。由式(3.6)和式(3.9)可知參考波場和格林函數(shù)均為脈沖源的波場響應(yīng)，因此當(dāng)震源函數(shù)的絕對值不是特別大或特別小時上述波場是大小相當(dāng)?shù)?。由?3.7)可知，若微擾，則顯然可得散射波場的值為0，相似地當(dāng)?shù)臑橹行狞c為0半徑很小的鄰域內(nèi)的值，若脈沖震源的響應(yīng)隨著震源的強度的變化速率是有限的，則可以相應(yīng)的推知此時的值也是很小的，此時散射波場相對于背景波場而言可以忽略，也即可以對式(3.11)進(jìn)行線性化處理。一般通常把式(3.11)的線性化處理過程稱作為波恩近似，可表示為(Esmersoy and Levy 1986; Huang and Schuster 2014; L

58、evy and Esmersoy 1988; Wei et al. 2014)：, (3.12)當(dāng)?shù)卣鸩ㄔ诘叵滤俣冉橘|(zhì)的分界面處，對于大于臨界反射角之外的反射波而言，由于此時的反射系數(shù)為1，這就意味著此時相對于界面的反射波振幅可以與入射到此界面的入射波相比較。因此，在高維情形下應(yīng)對是否能進(jìn)行波恩近似時行驗算，而不是直接進(jìn)行線性化處理。對于高頻情形，由式(3.6)及式(3.7)可知當(dāng)角頻率大于一定值的時候，參考波場與散射波場在數(shù)值程度上是可以相互比較的，也即此時相對于參考波場而言不能再一味的完全忽略掉散射波場，故在考慮能否用波恩近似來完成線性化時應(yīng)充分驗證角頻率的選擇是否能滿足進(jìn)行如此的簡化。在

59、高頻地震波及足夠小的微小擾動的情形下，利用WKBJ近似可以把方程(3.4)的解析解表示成如下的試驗解(Bleistein et al. 1987; Engquist and Runborg 2003)：, (3.13)其中為實數(shù)，是地震波傳播用時(由下文可以得到反映)，為待定函數(shù)，且其中的為非0量。對于式(3.4)中離開震源所在點之外的任意點的波場解，我們可以把式(3.13)所給出的試驗解看作是一個或多個波場解的線性疊加，對于離開震源點的波場的解析解可以通過把式(3.13)代回式(3.4)求得，也即：, (3.14)對上式進(jìn)行化簡，可得：. (3.15)若式(3.13)所表示的試驗解是式(3.

60、4)的解析解，則可得式(3.15)一定恒成立。若要使式(3.15)恒成立，由于其本身是關(guān)于的多項式級數(shù)，若使得每一項級數(shù)的值均為0，則式(3.15)能恒成立，也即可得(Bleistein et al. 1987; Cohen et al. 1986)：, (3.16) , (3.17). (3.18)式(3.16)是在折射率為n的介質(zhì)中的程函方程，由折射率的定義式(3.2)及式(3.16)的具體形式可以看出可以表示介質(zhì)中地震波的傳播走時(Brown 1982)，同時由式(3.4)可以看出乘積因子保證了在背景上具微擾的介質(zhì)模型中的程函方程能正確地表示介質(zhì)的全波速。式(3.18)表示在折射率為n的