NREE：基於地震波場的跨尺度地球內部成像
地震方法是如今人類認識地球內部結構的重要方法之一，也是人類了解地球深部構造的唯一方法。藉助地震成像方法，地球物理學家將複雜的地震波波形信息轉化為易懂的地球模型參數，如地球內部的速度結構。因此地震成像方法的精度直接決定了人類對地球內部結構的認識程度。
現今精度最高的地震成像方法是誕生於上世紀70年代的地震全波形反演方法。與傳統的射線方法不同，全波形反演方法是基於地震波波場的成像方法。它的終極目標是利用接收到的地震波全部信息來推斷地球的結構特徵和動力學信息，例如獲得油氣藏儲層參數、定位地下熱點以及監測地震活動和板塊活動中的應力場變化等。在過去10年間，得益於計算能力的提高和數值算法的進步，全波形反演方法在解決諸多實際問題中取得了成功。從小微尺度的超聲波醫學成像，到中等尺度的勘探地震成像再到全尺度的地球內部成像，全波形反演方法的應用場景跨越了足有9個量級的頻率及波長。全波形反演方法成功預測了Valhall油田氣層下部的新近紀冰川印記，位於庇里牛斯山脈西部的伊比利亞地殼推覆體，全球範圍內俯衝帶的俯衝結構以及非洲超級地幔柱的結構特徵。展望未來，隨著計算資源的提高及算法效率的提升，全波形反演方法將會呈現出更為精確的地球內部結構以及物理和化學模型參數。普林斯頓大學布萊爾講席教授Jeroen Tromp在Nature Reviews Earth & Environment期刊發表綜述文章，回顧了全波形反演方法的發展歷史、基本方法及若干跨尺度的應用實例並展望了該方法的發展方向及挑戰。
圖1 地震全波形反演方法主要發展線（Tromp, 2019）
一、基本原理及關鍵因素
全波形反演方法的基本原理是通過不斷比較模擬地震波形數據與觀測地震數據來迭代更新地球參數模型（一般為縱、橫波速度和密度）。當預測數據與觀測數據足夠接近時，我們認為此時的地球參數模型也就足以近似地球內部的真實情況。為了生成預測數據，我們需要考慮震源、初始模型及所用的正演模擬方法；為了比較預測數據和觀測數據的不同，我們需要選取合適的目標函數；最終數據域的殘差通過數據反傳映射到模型域並通過優化算法更新參數模型。這些也是決定全波形反演方法應用效果的關鍵因素。
1. 震源
一般而言，數值模擬所用的震源需要接近觀測數據的真實震源以確保數據殘差是由於模型參數的不精確引起的。在求解區域及全球尺度問題時，所用的震源一般為天然地震事件，相應的震源信息需要在反演之初求解（如CMT解）並在迭代過程中適時更新（震源和模型聯合反演）。勘探應用一般選取主動源數據，震源信息較易獲得。當震源信息缺失或者不夠精確時，也可考慮之前提到的震源與模型聯合反演。由于震源反演與模型反演存在一定程度的串擾，現在也發展出了一些不依賴于震源的反演方法（如雙差反演）來避免反演震源信息。
2. 初始模型
對於經典全波形反演算法而言，初始模型需要足夠接近真實模型以避免反演陷入局部極值。即基於初始模型得到的預測數據與觀測數據的波形（走時）差異要控制在半個波長範圍內。區域及全球尺度的應用往往有較為精確的一維甚至三維參考模型，這些模型對常用的長周期觀測數據來說已經是足夠精確的初始模型了。相反的，勘探應用中並沒有足夠精確的初始模型。因此，對勘探地震來說，建立初始模型是首要任務。一些針對全波形反演的長炮檢距、低頻（～1.5 Hz）數據採集工作已經在一些石油公司展開；同時一些可以突破半波長準則限制的目標函數也被陸續提出。
3. 正演模擬
對於給定的震源和初始模型，若要得到預測地震數據，我們需要模擬地震波的傳播過程即正演模擬。在這裡我們考慮兩個問題：選用的波動方程和求解該方程的數值算法。在很長的一段時間內，石油工業界主要採用時間域有限差分方法求解聲波波動方程。這主要是基於計算效率、算法複雜度及觀測數據的特點考量的。隨著勘探目標的轉移，應用於陸地資料處理的彈性波波動方程也逐漸得以應用。與工業界不同，地震學家偏好（黏）彈性方程來處理觀測到的橫波及面波數據。時間域的譜元法因其在處理起伏地表、流固耦合界面及核幔邊界（CMB）的高精度優勢也得到了地震學家的青睞。
4. 目標函數
作為衡量模擬數據與觀測數據近似程度的指標和模型更新方向的指示（伴隨震源），目標函數在某種程度上決定了反演的成敗。經典的全波形反演方法選用模擬數據與觀測數據的最小二乘差異作為目標函數，該目標函數只有在模擬數據與觀測數據較為接近（半波長準則）時才有唯一局部極小值，因此在應用局部優化算法時對初始模型和數據頻率範圍有較高要求。一般來說，我們需要根據需要匹配的數據特點來合理選取目標函數。例如，當觀測數據有明顯的離群點時，L1範數要優於L2範數（最小二乘）。在解決區域及全球的反演問題時，地震學家往往選取那些與觀測數據較為接近的模擬數據參與反演並且優先匹配地震數據中的相位（走時）信息。這是由於地震波振幅信息受震源的不確定性、記錄儀器的影響往往難以匹配。由於初始模型較為精確，全波形反演問題中常見的周波跳躍（非唯一局部極值）問題在區域及全球反演問題中並不突出。但該問題卻是自全波形反演方法誕生以來一直困擾勘探地球物理學家的棘手問題。相應地，勘探地球物理學家可以在反演之初選擇較為簡單的、與觀測數據接近的模擬數據參與反演，如部分初至數據並逐漸增加數據選擇範圍。同時，勘探地震學家致力於尋找新的目標函數來突破半波長準則的限制。如自適應波形反演方法、震源-檢波點拓展方法、時差拓展方法及最優輸運方法等。嚴格意義上講，現有的全波形反演算法並未使用全部的波形信息，而是使用儘量多的信息來約束反演模型。
5. 數據反傳
模型梯度為更新當前參數模型提供了更新方向。全波形反演方法可以使用弗雷謝特微分或者伴隨狀態法兩種方法求取模型梯度。在當前數以億計（大型三維模型）的模型參數條件下，計算每一模型參數的弗雷謝特微分並不可行。正是由於伴隨狀態法的提出，使得全波形反演算法實用化。藉此，模型梯度可以由正傳波場與反傳波場的互相關得到。所謂反傳波場即由目標函數計算得到的伴隨震源在檢波點處沿時間反傳得到。反傳所用的方程為正傳方程的伴隨形式，所用的數值算法為之前提到的時間域有限差分或者譜元方法。在實際應用中，由於正傳與反傳波場存在傳播時間上的差異，一般需要對正傳波場做存儲（硬碟）或者重建處理。
6. 優化算法
現有的全波形反演方法使用基於梯度的局部優化算法，比如最速下降法、共軛梯度法、牛頓法、高斯-牛頓法以及擬牛頓法。其中L-BFGS方法被認為是最為實用的擬牛頓方法。然而這些局部優化方法都面臨著同樣的問題：收斂到局部極值。當該反演問題有唯一局部極值時，該局部極值即為全局最優解，反之反演結果將會是其中某一個局部最優解。反演問題局部極值的多少與上文談到的初始模型、數據頻帶範圍、正演算子與實際物理過程的近似程度及所用的目標函數等因素直接相關。
全波形反演方法與傳統的基於射線理論的層析成像方法有著本質的區別。射線方法假設地球是簡單的球形對稱結構或者是近似的層狀結構，地震波速度的異常是基於一階擾動理論給出的。而全波形反演方法則可以適用於複雜的非均質模型，並考慮到了由此帶來的高階擾動。
圖2 全波形反演方法的基本流程（Tromp, 2019）
二、應用實例
地震全波形反演方法已經用於解決多種尺度的結構成像問題，如醫學超聲成像、無損探傷、地球近地表成像、勘探地球物理和區域及全球成像。所用的震源主要有主動源、天然地震以及環境噪聲等。在這些應用中，全波形反演的基本方法並沒有顯著的不同，但是根據所接收的數據特點及反演目標的不同，上文所述的六大要素需要進行相應的調整。
1. 勘探地震應用實例
過去十幾年間，全波形反演方法在勘探地震學領域產生了深遠影響。2007年盲測數據的成功應用給了勘探地球物理學家信心，自那時起對全波形反演方法的研究一直是勘探地球物理界的熱點領域。值得一提的是全波形反演方法在北海地區Valhall油田的成功應用。該油田位於北海70m深水區，由於淺部氣層的存在，常規的射線方法難以對氣層以下區域進行成像。石油公司在該地區採集了寬頻帶、寬方位角的海底電纜數據並應用多參數三維粘滯聲波全波形反演。反演所用的頻率範圍為3.5-10Hz，模型參數為各向異性縱波速度、密度及品質因子（Q）。如下圖所示，全波形反演方法精確刻畫出了淺層氣藏、古河道沉積以及古海底的冰川移動劃痕。
圖3 Valhall油田的多參數粘滯聲波全波形反演實例。左側為初始垂向縱波速度模型（由反射走時層析成像獲得）；右側為相應的更新模型。a-c為不同深度處的橫向速度切片；在反演結果的175m深度（a）切片處可見古河道沉積及生產平臺產生的散射點（X=6km，Y=11km）；在500m深度（b）切片處可見冰川移動在古海底留下的劃痕；在1km深度（c）切片處可見清晰的氣藏結構以及裂縫構造。相應的構造特徵在兩處縱向速度剖面（d、e）中也有體現（Operto et al.，2018）
2. 全球反演實例
首例全球尺度的全波形反演模型為GLAD-M15各向異性地幔模型。在2019年該模型被最新的GLAD-M25模型代替。該反演過程選取了1480個地震事件，在GLAD-M15模型的基礎上增加了10次擬牛頓迭代更新。GLAD-M25模型刻畫出了更為精細的地球內部結構：俯衝帶地區的細節更加豐富、地幔柱的橫波速度異常更加明顯。
圖4 哈特拉斯（Hatteras）俯衝帶縱波速度擾動模型切片。紅綠白圓圈標記了GLAD-M25模型及兩個參考模型（GAP-P4及UU-P07）速度切片的位置（Lei et al.，2020）
三、機遇與挑戰
在過去的四十多年間，儘管全波形反演方法取得了長足的發展，但其面臨的最大挑戰仍來自實際應用。一些實際應用中的棘手問題仍然存在，例如勘探問題中面臨的初始速度模型問題、低頻數據缺失問題以及全球反演中觀測數據分布不均等問題。作者總結了全波形反演方法在地球科學應用中面臨的機遇和挑戰。
1. 多參數反演
早期的全波形反演應用側重於反演單一的速度模型。然而地震波在地下的傳播過程還受到除速度以外的其他模型參數的影響，常見的有密度、傳播方向（各向異性）以及吸收衰減（Q）等。考慮多參數的全波形反演方法在模擬地震波傳播時更加接近真實的物理過程，這也就降低了由於傳播算子不精確引起的反演不確定性。在勘探地球物理界，聲波各向異性全波形反演已經得以應用，（黏）彈性全波形反演仍在探索中。對地震波的聲波近似並不適用於全球反演問題，因此在區域及全球反演問題發展之初就將（黏）彈性全波形反演作為主要發展方向，並逐漸將各向異性參數納入反演目標。多參數反演加劇了反演問題的病態程度。在有限的觀測範圍內（地表接收），模型參數間的串擾一直存在（如反射地震中的速度和密度），部分模型參數約束不足（如透射波對密度變化不敏感）。
2. 計算效率
有限的計算資源仍然是制約全波形反演應用的主要因素。計算效率的提升主要考慮計算機硬體系統的進步和算法的優化兩個方面。地球物理學家主要通過優化全波形反演算法及流程以提升計算效率，如平衡數值計算與輸入/輸出（I/O）負擔、圖形卡計算等。此外，反演結果的可視化對於解釋反演結果、科學發現有重要作用。
3. 反演結果的不確定性
受制於有限的觀測系統和地震波的傳播規律，反演模型的不同區域（如淺層和深層）、不同參數（如縱波速度和密度）的精確程度是不一樣的。對反演結果的不確定性分析將會幫助解釋反演結果、排除反演假象。目前對反演結果的不確定性分析仍處於探索階段。常用的數值分析方法是對棋盤模型的算例分析，但是該方法相當於進行一次全波形反演，耗費大量計算資源。基於貝葉斯理論的不確定性分析方法有望在未來得以發展並應用。
4.震源編碼
震源編碼方法可以極大地降低全波形反演方法的計算量，但是也容易引入不同震源波場和觀測數據波場之間的串擾問題。最新研究表明當地震波傳播至穩態狀態時，震源波場和反傳波場可以準確解碼，互相關得到的梯度更加準確。
5.全局搜索
如前文所述，現有的全波形反演方法依賴於基於梯度的局部優化算法。這也就不可避免地使反演過程收斂於某一個局部極值。基於貝葉斯理論的全局搜索算法有望跳出局部極值，找到全局最優解。然而這類方法的最大問題是需要對每一模型參數進行採樣（如正演模擬），計算量巨大。漢密爾頓馬科夫鏈蒙特卡洛方法因其可用伴隨狀態法求解有望降低計算量，使得全局搜索變得可行。
四、結論及展望
全波形反演方法經過四十多年的發展已經在科學界和工業界產生了巨大的影響。儘管該方法的成功應用仍然面臨挑戰，但它帶來的回報吸引了來自科學界和工業界的持續投入和研究。對地震學家而言，全波形反演提高了地球內部成像精度、提高了人類對地球內部物理和化學過程的認識。對工業界而言，全波形反演幫助開發現有油藏、發現新的油藏，為石油公司帶來豐厚的收入。展望未來，無串擾的震源編碼方法、基於漢密爾頓蒙特卡洛方法的貝葉斯反演框架、機器學習及數據分析方法將會對全波形方法產生深遠影響，進一步幫助人類認識地球內部結構。
原文：Tromp J. Seismic wavefield imaging of Earth’s interior across scales[J]. Nature ReviewsEarth & Environment, 2019: 1-14.
（編譯：張振東，劉伊克/油氣室）
校對：張騰飛