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

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

3、度分布情況下，可以采用比Sirgue和Pratt方法更少的頻點來進行頻域反演，而不會影響反演的精度，從而可以進一步提高反演計算的效率。本文針對一維速度模型、Marmousi模型和Overthrust三個模型進行反演實驗，模擬算例結果驗證了該優化的頻率點選取方法的有效性。關鍵詞 頻率點選取 ；波數覆蓋 ； 線性遞增模型 ； 頻域波形反演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緒論隨著人們對地下非均勻介質認識的不斷深入以及計算機性能的大幅度提升，地震波反演成像技術也相應地得到了非常迅猛地發展。然而，在使用傳統頻率點選取策略所得頻率點進行頻率從低到高的迭代頻域波形反演的前提下，在勘探的目標區域尺度非常大的情況下，計算效率就會成為此種反演方法所面臨的挑戰。此時，基于一種新的梯度圖像中波數覆蓋的頻率點選取策略可能會為地震波形反演帶來新的發展。本章將回顧相關的研究背景及現狀，簡述傳統的頻率點選取策略的原理，闡述它在實際應用中的局限性及其改進的可能性，給出本文的總體研究思路和主要研究

9、內容。1.1研究背景和意義地震勘探(seismic exploration)是地球物理勘探中的一種重要方法，其主要目標是使用人工震源激發的地震波來定位地下介質中的礦藏(如石油，天然氣，煤，礦石，礦物，地熱能源，水資源等)，獲得地下的巖石特性、碳氫化合物屬性等地質信息。在陸地地震勘探中，由地表或近地表的人工震源激發的地震波經由地下介質向地球內部傳播，當遇到地下介質的物性參數快速變化區域時地震波便會產生散射波(包括反射波，折射波，繞射波等)，一部分返回到設置在地表的高靈敏度檢波器的散射波會被此檢波器記錄下來，所有的這些被不同檢波器(hydrophone)記錄的地震記錄構成地震數據集；同理，在海洋地

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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