基于非常快速模擬退火算法的地震譜反演方法①

張廣智，李 寧，李 超，陳懷震，陳 雷

（中國石油大學地球科學與技術學院，山東 青島 266555）

0 引言

隨著石油勘探開發的不斷深入，地震勘探已經從簡單的構造識別轉向復雜地表、復雜構造、裂縫、薄儲層和老油田剩余油的勘探。勘探難度的加大也對能夠得到具有更高信噪比和分辨率的地震資料提出了挑戰。地震子波大約在30Hz左右，對于25m以下的薄層很難區分它的頂部和底部［1］。但是工業的目標地質體的厚度一般在10m甚至更小。反演薄層在儲層識別和油藏描述中起到了重要的作用［2－3］。

Charles I．Puryear和 John P．Castagna［4］在2008年提出了譜反演的方法。譜反演是一種基于譜分解技術，在時頻域建立目標函數的一種地震反演技術，它能夠反演出小于調諧厚度的薄層。在過去的幾十年里人們一直利用 Widess模型［5］理論得到的調諧厚度分析方法指導薄層的識別。Partyka等［6］以及 Marfurt和 Krilin［7］采用基于離散傅里葉變換得到的頻譜分解結果來求取薄層厚度的方法。但是如果地震頻帶寬度不能夠清晰識別譜峰和陷波變化規律時，譜分解對于分辨薄層還是存在困難。這 也 推 動 了 Partyka、Portniaguine、Puryear 和Castagna等人對譜反演的探索和發現。

由于地球物理問題中普遍存在的非線性和不適定性。本文采用全局尋優的非常快速的模擬退火算法對上述的時頻域目標函數求最小值。該方法不依賴初始模型，能夠有效的跳出局部極值，并且對先驗信息獲知較少的情況下也能得到較好的結果。

1 譜反演原理

譜反演研究思路如圖1所示。

圖1 譜反演流程圖Fig．1 The flowchart of spectral inversion．

譜反演方法是一種不同于傳統地震反演的新方法，可以用于薄層成像。譜反演是利用時頻分析和譜分解技術獲取地震資料的局部頻譜信息，首先在時域將反射系數分為奇偶分量，然后在頻域建立目標函數，再利用改進的模擬退火算法求解。其中關鍵技術是目標函數的建立和求解。在理想的情況下可以識別任意薄層，并且精確的反演出反射系數的值。該方法不需要任何先驗模型、任何反射系數假設、任何層位約束、也不需要井強制約束等優點。下面重點介紹目標函數的建立和模擬退火的尋優算法。

1．1 目標函數的建立

目標函數是由頻率域的地震數據、地震子波及時域的反射系數的奇、偶分量組成。如果反演過程中反射系數正確，目標函數能夠達到最小值，作為判斷反演結果正確與否的依據。

基于圖2的多層模型建立目標函數，經過復雜推導，最后整理得目標函數表達式為［4］

式中S（t，f）、W（t，f）分別為短時傅里葉變換地震記錄的振幅譜和子波的振幅譜；re（t），ro（t）為圖2模型中反射系數對的偶分量和奇分量；ae，ao為權衡系數，依據信噪比的高低來確定二者的比值；fl，fh分別為反演頻帶的低截頻和高截頻；tw為反演的半個時窗長度。Ti為反射系數對之間的時間厚度。紅色三角表示分析點的位置。

圖2 多層模型Fig．2 Multi－layers model．

構建目標函數理論是基于任何一個反射系數序列都可以分解成奇偶分量的和。奇分量具有相等的大小和不同的符號，偶分量具有相等的大小和相同的符號；奇分量不利于檢測薄層，而偶分量可以提高分辨薄層的能力。譜反演的實質就是利用偶分量在厚度趨于零時的有效干涉來提高地震資料的分辨率，同時偶分量也有利于抗噪。

目標函數是通過結合小時窗內地震資料時頻分析的結果和反射系數的奇偶分量，類似于頻率域褶積原理建立起來的。譜反演的目標函數與時間域的褶積殘差目標函數相比具有更好的收斂性和約束能力，能夠有效的減少地震反演中的多解性問題。在反演的過程中因為地震資料是帶限的，所以不能在全頻帶進行反演，只能選取信噪比較高的優勢頻帶進行反演。因此不能得到類似于脈沖信號的反射系數，只能得到一個提高分辨率的剖面，但它代表一定的反射系數的信息。

1．2 模擬退火算法

模擬退火法（Simulated Annealing Method，簡稱SA），最早是由 Metropolis等人［8］在1953年提出，1983年 Kirkpatrick等［9］成功將它應用于組合優化問題中。Rothmann最早將模擬退火算法引入到地球物理問題中，解決地震資料的剩余靜校正問題。隨后人們開始利用其進行一維地震波形反演，一維電測深反演以及疊前速度分析［10］。由于模擬退火的計算過程簡單、通用，魯棒性強，適用于并行運算，不過多的依賴初始模型，能夠有效的跳出局部極小，所以它成為人們解決非線性優化問題常用的算法［11］。

Kirkpatrick提出的算法是以Gaussian概率分布的準平衡Boltzmann－Gibbs統計為基礎的，為經典的模擬退火算法（MSA）；Szu和 Hartley用Cauchy－Lorentz分布代替Gaussian分布，形成快速模擬退火算法（FSA）；Stariolo和Tsallis泛化FSA方法，提出廣義模擬退火算法（GSA），包含了CSA和FSA，計算速度比FSA 快；Ingber［12］則采用了依賴于溫度的似Cauchy分布產生新的擾動模型，提出非常快速的模擬退火算法（VFSA）；張霖斌等［13］以廣義Boltzmann－Gibbs統計理論為基礎，采用依賴于溫度的似Cauchy分布產生新的擾動模型，建立一種新的快速模擬退火算法。

陳華根等［14］深入分析VFSA的核心技術，發現VFSA存在的缺陷，將VFSA的退火計劃與模型擾動相匹配，提出了改進的模擬退火算法（MVFSA）。思路是：算法開始在降溫計劃的控制下，模型作全局擾動以搜索并鎖定最優解區間；在鎖定最優解空間后，在新的降溫計劃下，模型在已被接受模型周圍做局部擾動，以近似局部搜索的模擬退火方式，逐步逼近最優解；新的退火計劃作適當的回火升溫，給予已尋找到的最優解空間再一次跳出局部極小值的機會，使得找到的最優解更可靠。同等條件下，MVFSA所花費的時間要比VFSA少20%～30%，算法顯得更穩健。

模擬退火算法之所以能夠稱為全局尋優算法，是因為它在搜索過程中以一定的概率接受“惡化解”，溫度較高時，接受的概率近似為1，隨著溫度的降低，接受“惡化解”的概率逐漸降低，最后變為局部尋優。

圖3 全局搜索與局部搜索Fig．3 Global search and loacal search．

隨著模擬退火算法應用的廣泛，人們對其進行了多種改進，現在應用較多為非常快速的模擬退火算法。本文采用了VFSA算法并將其和經典的模擬退火算法在機理上進行了比較研究［13－14］。

經典的模擬退火算法多采用均勻的模型擾動方式即：mj＝mi＋step×（rand－0．5）。這里mj表示新解；mi表示前一個解；step表示步長（通常與溫度無關）；rand表示0～1之間的隨機數。而VFSA的模型擾動方式采用依賴于溫度的似Cauchy分布法［12］：

由Boltzmann－Gibbs［14］分布，可得出新的接受概率計算公式為

本文中h＝－5；ΔE擾動前后的能量差；T為當前的溫度。T的計算公式為

式中T0為初始溫度，本文中設為T0＝100；α＝0．96；c＝1；N＝2。

圖4 模型擾動量和接受概率隨迭代次數的變化Fig．4 Variation of the model turbulent and acceptance probability with iteration．

經典的模擬退火算法與VFSA的模型擾動量、接受概率隨著迭代次數的變化關系如圖4所示。可以看出VFSA在高溫時搜索范圍比MSA大，在低溫時搜索范圍比MSA小，這樣更有利于跳出局部極值，加快了收斂的速度。VFSA接受概率隨迭代次數的變化較MSA要快，并且曲線的尾端較MSA略有上翹，這是由于在低溫時VFSA的擾動量變小，使得反演解更靠近全局最小解，反演精度更高。

利用好模擬退火這個全局尋優的優化算法需要模型擾動、接受概率、降溫函數3者之間的有效配合。

2 模型試算

2．1 單道多層模型試算

主頻為30Hz，4ms采樣的雷克子波與反射系數模型（圖5（b）1）褶積合成地震記錄，加上5%的高斯隨機噪聲，得到觀測值（圖5（a）1）。利用非常快速模擬退火算法對依式（1）建立的目標函數進行求解。由于噪聲的存在，且從噪聲分布的頻帶上看，噪音使得高頻帶的地震記錄信噪比變低，使用5～90 Hz頻帶信號參與反演。

圖5 觀測與反演的地震記錄和反射系數比較Fig．5 The seismic waves and reflectivities from record and the inversion．

反演地震記錄（圖5（a）2）與觀測地震記錄的相對誤差為0．246%，反演反射系數（圖5（b）2）與真實反射系數的相對誤差為3．814%。通過圖5可以看出利用譜反演方法能夠準確的反演出反射系數的大小和位置。不能夠在地震記錄中分辨的薄層，能夠在反射系數中得到分辨，提高了分辨率。

2．2 楔形模型試算

主頻為30Hz，4ms采樣的雷克子波與反射系數模型褶積，加上5%的高斯隨機噪聲，得到圖6（a）的楔形層狀地震記錄，作為反演輸入。

圖中黑色實豎線的位置即利用λ／4確定的調諧厚度的位置。從反演結果（圖6（b））與輸入地震記錄的比較，可以看出譜反演可以識別出小于調諧厚度的薄層。在250～350ms（矩形框）的位置存在3個反射面，在地震剖面上不能清楚的分辨，但是從反演的結果能夠清楚的識別。理論上譜反演可以反演出任意厚度的薄層，實際由于噪聲的存在只能得到類似反射系數的剖面。

圖6 模型數據的地震記錄與譜反演的結果Fig．6 The model of seismic record and spectral inversion result．

針對上述楔形層狀模型的地震記錄、反演結果、真實反射系數做出頻譜圖（圖7）。曲線1為真實地震記錄頻譜，可以看出原始地震記錄帶限，缺失高頻信息。曲線3為真實反射系數的頻譜頻帶最寬，是最理想的反演結果。由于噪聲的存在和地震記錄的帶限性，譜反演沒有恢復出全頻帶的反射系數，但是對于地震記錄的拓頻作用明顯，極大的提升了60～90Hz頻帶范圍內的有效信號。

圖7 地震記錄的頻譜、反演反射系數的頻譜與真實反射系數的頻譜的比較Fig．7 Comparison of spectrums of seismic record、inversion reflectivity and true reflectivity．

2．3 marmous模型

下面截取部分marmous模型（縱波速度）求取反射系數與4ms采樣的雷克子波褶積合成地震記錄。如（圖8（a））所示，在150道至350道間有一透鏡狀含氣砂體，但是由于子波的影響，其邊界不是很明顯。在已知地震記錄和準確子波的情況下，對地震道進行譜反演試算。

通過譜反演試算，得到結果（圖8（b））。可以看出譜反演能夠有效的去除子波的影響，清楚地識別出砂體及其他反射體的分界面；能夠準確的反演出反射系數的位置和大小；極大的提高了數據體的分辨率，有利于精細的地震解釋。

3 實際數據處理

該工區為加拿大Alberta白堊系淺層的部分地震數據，2ms采樣，在第60～80道0．02～0．03處存在一個氣層。應用測井方法從實際數據中提取較為準確的子波，然后對其進行譜反演（圖9）。

圖8 原始疊后地震剖面和反演結果Fig．8 The original seismic record and spectral inversion result．

圖9 加拿大某工區原始疊后地震剖面和譜反演結果Fig．9 The original seismic record and spectral inversion result at one site in Canada．

通過反演前后數據的比較，可以看出數據的分辨率得到了顯著的提高，同相軸較反演前變細，邊界更加明顯。地震解釋結果比原始數據更加可靠。

4 結論

通過對譜反演方法的推導以及模型數據和實際資料的處理，可以得出以下幾點結論：（1）全局搜索的模擬退火算法可以有效的跳出局部極值，收斂到全局極小，但是收斂速度慢。本文采用的改進的模擬退火算法VFSA提高了搜索過程中的收斂速度，并且保證了解的正確性和精度。（2）譜反演作為一種精細識別薄層、反演反射系數的方法，可以準確的反演出反射系數的大小和位置，其目標函數與時間域的目標函數相比具有更好的收斂性。（3）雖然使用部分頻帶信號參與反演，但是從反演的效果看反演有效地去除了子波的影響，明顯的拓寬了原始資料的頻帶，極大地提高了資料的分辨率，對于精細的地震解釋起到了作用。

在反演過程中，要反復的對地震數據和子波進行傅里葉變換，因此譜反演相對時域的其它反演方法更費時；因為要變換到頻率域處理，所以該方法對于噪音較敏感。譜反演方法在國內雖然還不十分成熟，但是對于薄層的識別，地震資料分辨率的提高具有深刻意義。

［1］ Chopra S，J P Castagna，O Portniaguine．Seismic resolution and thin－bed reflectivity inversion［J］．Canadian Society of Exploration Geophysicists Recorder，2006，31：19－25．

［2］ 秦德文．基于譜反演的薄層預測與反演方法研究［D］．青島：中國石油大學（華東），2009．

［3］ 陳科．基于模擬退火的譜反演方法研究［D］．青島：中國石油大學（華東），2010．

［4］ Charles，Puryear，Castagna John P．Layer－thickness determination and stratigraphic interpretation using spectral inversion：Theory and application［J］．Geophysics，2008，73（2）：37－48．

［5］ Widess M．How thin is a thin bed？［J］．Geophysics，1973，38：1176－1180．

［6］ Partyka G A，Gridley J A，Lopez J A．Interpretational aspects of spectral decomposition in reservoir characterization［J］．The Leading Edge，1999，18（3）：353－360．

［7］ Marfurt K J，Kirlin R L．Narrow－band spectral analysis and thin－bed tuning［J］．Geophysics，2001，66：1274－1283．

［8］ Metropolis N，Rosenbluth A，Rosenbluth M，et al．Equation of state calculations by fast computing machines［J］．J．Chem．Phys，1953，21：1087－1092．

［9］ Rothman D H．Nonlinear inversion statistical mechanics and residual statics estimation［J］．Geophysi－cs，1985，50（12）：2784－2796．

［10］ 師學明，王家映．模擬退火法［J］．工程地球物理學報，2007，4（3）：165－174．

［11］ 康立山，謝云，尤矢勇，等．非數值并行算法－模擬退火算法［M］．科學出版社，2003．

［12］ Inger L．Very fast simulated annealing［J］．Math Conput Modeling，1989，12（8）：967－973．

［13］ 張霖斌，姚振興，紀晨，等．快速模擬退火算法及應用［J］．石油地球物理勘探，1997，32（5）：654－661．

［14］ 陳華根，吳健生，王家林，等．模擬退火算法機理研究［J］．同濟大學學報（自然科學版），2004，32（6）：802－805．