多尺度全波形反演的正則化思想及速度模型的多尺度分解

肖金梅，曹靜杰，楊賀龍

河北地質大學 河北省戰略性關鍵礦產資源重點實驗室，河北 石家莊 050031

0 前言

自從Tarantola[1]提出基于廣義最小二乘反演理論的時間域全波形反演方法以來，全波形反演逐漸成為獲得高精度速度模型的重要手段，為改善成像效果、區域深部構造及演化分析、淺表層環境調查、宏觀速度場建模與成像、巖性參數反演提供了有力工具。全波形反演不僅僅是為逆時偏移等偏移方法提供高精度的波場，其更深遠的目的是直接由疊前波場數據獲得速度，密度等多種參數，用于儲層預測和描述，因此目前的研究還有很大的發展空間。全波形反演方法基于疊前地震波場的運動學和動力學信息來重建地層結構，具有揭示復雜地質背景下構造與儲層物性的潛力，比常規的層析成像方法獲得速度參數更加精細和準確。全波形反演典型的病態地球物理反問題。反演問題一般存在解的非唯一性、對噪聲敏感性、對初始模型強依賴、易陷入局部極值、計算量大等問題。全波形反演求解困難從理論上說主要是因為其高度非線性，存在多個極值和觀測數據的不完備性，對低頻信息的依賴性較大，需要含有低頻信息的速度模型作為初始解。由于對初始模型的高度依賴性，學者們提出了很多方法來得到好的速度模型，比如旅行時反演，偏移速度分析，Laplace域反演[2]，包絡反演[3]，正則化的積分方法[4]等；另一種獲得理想解的方法是多尺度全波形反演[5]，Sirgue 和 Pratt[6]證明了可以從低頻波場開始進行反演，首先恢復與低頻波場信息對應的大尺度速度結構，然后逐漸提高頻率成分，用低頻信息得到的速度模型作為下一個高頻反演的初始模型，反演出更加精細的速度結構。目前的采集環境噪聲的能量隨著頻率的減小而呈指數增加，真實數據的低頻分量被認為是不可靠的，因此很難產生低頻震源和采集到低頻成分，地球物理學家們處在一個兩難的境地。隨著以小波變換為代表的多尺度變換的興起，基于多尺度變換的反演方法是實現全波形反演的好的策略，由于低頻波場信息對解的敏感度較低，因此可以將觀測數據分成不同頻率的分量，首先反演出大尺度的速度模型，然后逐漸降低尺度的大小，得到小尺度的速度結構，隨著尺度的增多就可以獲得較好的反演結果[7]。

本文首先介紹了現有的主要正則化方法，然后給出了多尺度反演的理論公式，分析其正則化思想，指出多尺度反演是一種正則化方法；然后介紹了多尺度全波形反演的實現思路。速度模型的多尺度分解是多尺度全波形反演中重要組成部分，Curvelet變換[8]和Dreamlet變換[9]是在地震勘探領域應用較多的多尺度變換，耿瑜等[9]對比分析了Curvelet變換和Dreamlet變換對波場數據的分解，表明Dreamlet對波場數據有更稀疏的表達。本文對鹽丘速度模型進行Curvelet多尺度分解和重建，表明Curvelet分解的大尺度模型對應于速度的大尺度構造，小尺度模型對應于速度的精細結構，因此是速度模型多尺度分解的一個很好的選擇。本文對不同尺度鹽丘模型的正演進行結果頻譜的比較，發現在近偏移距大尺度模型正演數據的頻譜成分和真實模型正演數據的頻譜成分類似，在遠偏移距，真實速度模型正演數據的低頻成分缺失。研究為多尺度全波形反演的速度模型分解和實現提供了有價值的研究思路。

1 多尺度反演方法的正則化思想

假設正演問題為L（m）=d，其中L是一個非線性算子，m為模型參數，d為數據參數。反演問題一般是由數據參數d和算子L來獲得模型參數m。反演問題的解可能是非唯一的，不存在的或者不穩定的，地球物理反問題的上述特性使得反演問題的求解比一般的正演問題的求解更加困難，因此反演問題的求解一直是地球物理學家們的研究熱點。求解反演問題的一個重要思想是正則化策略，通過對解的先驗信息的探索，利用合理的先驗信息來求解，或者改變正演算子的性質。正則化方法可以分為很多種，在地球物理領域，經典的方法是Backus-Gilbert方法[10]，在數學領域，則以Tikhonov正則化最為著名[11]。下面對一些重要的正則化思想進行簡單闡述。早期的反演理論主要針對線性反演問題的求解。通過對算子L線性化，先對L算子進行泰勒展開，

然后采用非線性優化方法求解。實際上由于L算子可能是無窮次可微的，往往采用高斯牛頓法，共軛梯度法或速下降法來求解，然而這些算法屬于局部優化方法，對于多極值問題，該方法往往容易陷入局部解， 這就是著名的Tikhonov正則化方法的思想[11]。除了對解的二范數約束，最近幾年在信號領域興起了以解在變換域的稀疏性為先驗約束的壓縮感知理論[12]，使得1范數正則化得到了前所未有的應用和研究[13]。全變差也是一種稀疏約束的正則算子，此外，無窮范數在某些應用問題中也具有重要的物理意義[12]。

另一種常見的正則化方法是截斷奇異值分解方法，該方法對近似線性算子L的奇異值分解，考慮到大的條件數對解的穩定性的影響，以及小奇異值對解的影響，通過強制小奇異值為零來獲得穩定的解。截斷奇異值分解方法是解線性反演問題的一種重要的正則化方法[12]。

除了以上兩種正則化方法，預條件方法也是一種正則化策略。通過對算子的預處理來改變算子L的性質，使得求解變得更加穩定。預條件共軛梯度法是預條件正則化的代表[12]。

正則化方法的求解可以分為迭代正則化和直接正則化，比如共軛梯度法屬于一種迭代正則化方法，而截斷奇異值分解則屬于直接正則化[10]。多尺度反演是類似于預條件方法的一種正則化方法，通過模型和數據的多尺度分解，首先求解出大尺度模型，然后逐漸更新大尺度速度模型來逐漸得到準確的模型參數，下面主要分析多尺度反演方法的正則化思想。

多尺度反演首先對測量數據d進行多尺度分解，

其中Φ1，Φ2，…，ΦN是不同尺度的變換，N表示分解成的尺度的個數。假設Φ1，Φ2，…，ΦN之間是正交的，對于線性問題 Ψm=d，變成

基于傅里葉變換的從低頻到高頻的求解方法[6]和基于小波變換的多尺度反演[7]都屬于這個思路。這個思路依賴于第一方程組中算子的病態性的程度。即使Ψ是病態的，但是復合算子ΦT1Ψ從i=N到i=1的病態性變弱，因此可以從i=1開始逐步求解，這就是是多尺度反演的正則化思想和原理。多尺度反演減弱了對初始解的依賴，是解病態問題的重要方法和策略。

將上述過程寫成矩陣向量乘積的形式

因而

因此方程（6）可以變為,

上式寫成

從形式上說，正定或欠定的問題變成了超定的問題。為了得到穩定的解，先算第一組，得到好的解，作為第二組的初始解，然后依次類推。

對于非線性問題，同樣可以構造以下函數

結合非線性問題的線性化，可以得到好的解。該策略關鍵是使得條件數減小，改變算子的性態，減小對初始解的依賴性。

對于多尺度反演來說，如果變換后的d1，d2，…，dN都是稀疏的，有利于加快反演的計算。實際上，多尺度全波形反演不僅僅可以對觀測數據分解，速度模型和正演算子同樣也可以采用多尺度分解的方法。從物理意義上說，對模型的多尺度分解更有意義，因為速度結構可以看作是由宏觀構造和精細構造組成的。多尺度反演方法依賴于多尺度分解后數據的稀疏性，如果能用少量的系數表達原始空間的數據，則在變換域中進行反演能夠提高計算效率和精度。

假設模型按照如下方式分解

其中s1，s2，…，sN代表了不同尺度的模型在變換域的系數，M為速度模型多尺度分解的個數，Φ1，Φ2，…，ΦN組成多尺度變換。則問題（6）變為

上式等價于

由于si代表不同尺度的信息，對解的貢獻量最大，因此可以先求出大尺度的解，然后逐漸逼近精細結構解。多尺度反演更進一步的研究是將數據d，模型m和算子L都進行分解，然后利用分解后數據的稀疏性和解的多尺度性進行求解，加快收斂速度。

2 多尺度全波形反演的理論

在最小二乘框架下，全波形反演的質量依賴于初始模型的質量，即好的初始模型能夠避免陷入局部最優解。但是基于廣義最小二乘模型的全波形反演是多極值的，如果不改變這個模型，則只能依賴好的初始模型來突破。因此全波形反演的另一個可行的思路就是建立關于速度模型的單極值最小化模型，這樣就會避免對初始解的強烈依賴。此外，正則化方法通過改變算子的性質來使得速度模型容易求解，是目前全波形反演的主要研究思路之一。全波形反演的速度模型m和地震數據d之間是非線性關系：

其中的算子L代表基于速度和子波的波動方程正演過程。根據上面的分析，多尺度全波形反演可以表示為如下的模型：

這里模型和數據之間的關系是非線性的，并且只對模型進行分解。

多尺度全波形反演的關鍵是求出大尺度的速度模型，由于大尺度的速度模型可以在某個變換域用很少的系數表示，因此s1是稀疏的，在利用多尺度全波形反演求解大尺度模型時，只要求出少量的系數s1就能夠表達速度模型的大尺度信息，這對于提高計算速度非常有利。因此在求解時，先求出大尺度模型的系數s1，然后逐漸更新速度模型逼近真實解。由于大尺度模型對解的貢獻大，小尺度模型對應了模型的精細結構，對解的貢獻小，從大尺度的速度模型開始計算可以避免陷入局部最優解。實現大尺度系數的求解是多尺度全波形反演的難點，如果對正演算子線性化，將問題（15）變為線性問題，降低求解難度。從理論上，多尺度全波形反演可以直接獲得大尺度的速度模型，降低對速度模型的依賴，因此是一種很有前景的全波形反演方法。

3 速度模型的多尺度分解

根據上面的介紹，速度模型和地震數據的多尺度分解是多尺度全波形反演的重要組成部分。目前, 在地震勘探領域采用的多尺度變換有Curvelet變換和Dreamlet變換，耿瑜等[9]對地震數據的多尺度分解的研究表明Dreamlet變換比Curvelet變換更加適合數據的分解。當對觀測數據和速度模型都進行多尺度分解時，可以對觀測數據采用Dreamlet變換。本文對速度模型采用Curvelet變換分解，探討Curvelet變換分解速度模型的合理性。 采用的模型為SEG_EAGE 鹽丘模型，模型大小為nx=700, nz=150, 空間網格大小為12米，原始模型見圖1a。Curvelet變換將速度模型劃分成6個尺度，其中第一個精細尺度劃分成8個方向，圖1b給出了不同尺度的系數的對數形式，可見大尺度模型的系數很少，而小尺度模型的系數很多。Curvelet的原子見圖2。全波形反演的關鍵是大尺度速度模型的建立，因此圖3給出了最大尺度的模型及其所對應的系數，大尺度模型在變換域只要用203個系數就能表達，這對于模型的求解是很有意義的。圖4給出了兩個最大尺度模型的疊加和兩個最小尺度模型的疊加，可見兩個大尺度的模型的疊加是真實速度模型的好的近似，而小尺度的模型能夠刻畫模型的細節，類似于反射系數。因此在求解時，為了避免小尺度模型的干擾，應先求出大尺度模型，然后逐漸逼近真實的速度模型。由此可見，對速度模型進行多尺度分解是合理的，大尺度的模型表示了速度的宏觀結構，小尺度模型對應了速度的精細結構，因此可以采用Curvelet變換進行速度模型的多尺度分解。

圖1 a為SEG-EAGE 鹽丘模型；b為雙對數坐標下Curvelet系數的向量形式。Fig. 1 a SEG-EAGE model; b Vector form of curvelet coefficient in double logarithmic coordinates.

圖2 a-f為時空域下最大尺度到最小尺度的Curvelet 原子。 Fig. 2 a-f Curvelet atoms from maximum scale to minimum scale in spatiotemporal domain.

圖3 最大尺度的鹽丘模型及其對應的Curvelet系數（203個）。Fig. 3 Maximum scale salt-dome model and its corresponding Curvelet coefficients （two hundred and three）.

圖4 a為兩個最大尺度組成的速度模型； b為兩個中間尺度組成的速度模型。Fig. 4 a A velocity model composed of two maximum scales; b Velocity model composed of two intermediate scales.

4 大尺度速度模型的正演

由于全波形反演的關鍵在于大尺度模型的獲取，為了更好的研究大尺度模型的地震響應，下面對上一節的原始鹽丘速度模型和最大尺度的模型分別進行正演，分析正演結果的特點。圖5a是對原始速度模型正演的一炮記錄，炮點在模型中間，震源采用主頻為15 Hz的雷克子波，時間采樣率為1毫秒。圖5（b）給出了該炮記錄的F-K譜。圖6給出了最大尺度模型正演的炮記錄和F-K譜。由F-K譜可以看到，大尺度模型正演的頻率成分較低。圖7畫出了上述兩種炮記錄的遠偏移距的地震道及其振幅譜。在遠炮檢距的第50道以及650道，大尺度模擬的數據在遠偏移距低頻成分相比而言要更強，真實模型模擬的波場數據中低頻成分很弱，這是地震數據低頻成分缺失的表現之一。圖8給出了近炮檢距的地震道的振幅譜，對于近偏移距的300和450道，兩種炮記錄的頻譜很類似，不存在明顯的低頻缺失。因此，首先采用近偏移距數據反演大尺度模型可能獲得較好的效果。

圖5 原始鹽丘模型正演的炮記錄及其對應的F-K譜。Fig. 5 Gun records of forward modeling of original salt dome model and its corresponding F-K spectrum.

圖6 大尺度模型正演的炮記錄及其對應的頻譜。Fig. 6 Shot records of large scale forward modeling and their corresponding spectrum.

圖7 a為原始鹽丘模型正演炮記錄遠炮檢距的地震道；b為a的振幅譜；c為大尺度模型正演炮記錄遠炮檢距的地震道；d為c的振幅譜。Fig. 7 a Seismic trace of far offset recorded by forward shot of original salt dome model; b Amplitude spectrum of a; c Seismic traces of long offset recorded by large scale forward modeling; d Amplitude spectrum of c.

圖8 a為原始鹽丘模型正演炮記錄近炮檢距的地震道；b為a的振幅譜；c為大尺度模型正演炮記錄近炮檢距的地震道；d為c的振幅譜。Fig. 8 a Near offset seismic traces recorded by forward shot of original salt dome model; b Amplitude spectrum of a; c Seismic trace of near offset recorded by forward shot of large scale model; d Amplitude spectrum of c.

5 結論

本文簡單介紹了現有的正則化方法，重點分析了多尺度反演的正則化思想，指出多尺度反演可以看作一種正則化方法，通過不同尺度的變換來改變原始算子的病態性，降低對初始解的依賴，多尺度全波形反演作為多尺度反演的特例自然也是一種正則化方法。本文還簡單介紹了多尺度全波形反演的實現原理和優越性。通過對鹽丘速度模型的Curvelet多尺度分解和重建，表明對速度模型的多尺度分解是合理的，因此Curvelet變換可以作為速度模型的多尺度變換。通過對大尺度鹽丘模型和真實鹽丘模型的正演模擬表明，在近偏移檢距，大尺度模型正演數據和真實模型正演數據的頻譜成分相似，兩者的低頻分量都存在，在遠偏移距，大尺度模型正演數據的低頻信息豐富，但是原始速度模型正演的數據的低頻信息缺失。因此多尺度全波形反演可以從近偏移距波場數據開始多尺度反演，提取大尺度的速度模型信息，然后逐漸增大偏移距，利用大偏移距波場信息提高模型的精度。本文的研究對速度模型的多尺度變換的選擇和多尺度全波形反演的實現提供了有價值的研究思路。