祝賀君，劉沁雅，楊繼東

1 美國得克薩斯州立大學(xué) 達拉斯分校 地球科學(xué)系，美國得克薩斯州理查德森 75080

2 加拿大多倫多大學(xué) 物理與地球科學(xué)系，加拿大安大略省多倫多市 M5H 2N2

3 中國石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院，青島 266580

0 引言

研究地球內(nèi)部物質(zhì)結(jié)構(gòu)及其與板塊運動的相互作用是地球科學(xué)的一個重要研究內(nèi)容.目前由于鉆探技術(shù)的限制，無法直接測量大部分地球深部的物理參數(shù)，包括地震波傳播速度、電導(dǎo)率、溫度和巖石成分等.所以我們需要使用間接的物理和化學(xué)方法來探測地球甚至于其他行星的內(nèi)部構(gòu)造.由于地震波傳播的穿透性，地震成像給我們提供了一個很好的間接探測工具（Dziewonski and Romanowicz,2015;Liu and Gu,2012;Nolet,2011;Rawlinson et al.,2010;Ritsema and Lekic,2020;Romanowicz,2003;Tromp,2020）.從1970 年代開始，地球物理學(xué)家們就嘗試使用地表記錄的地震信號來推測地球內(nèi)部的溫度以及物質(zhì)組成的三維分布，以用于研究地幔對流、巖石圈演化和板塊俯沖等重要的地球科學(xué)問題（Aki et al.,1977;Dziewonski and Woodhouse,1987;Grand et al.,1997;Masters et al.,1996;Su et al.,1994;Van der Hilst et al.,1997;Woodhouse and Dziewonski,1984;Zhao,2004）.這種使用地表記錄的地震波信號來反演地球內(nèi)部三維結(jié)構(gòu)的方法被稱為地震層析成像，該方法與醫(yī)學(xué)和聲學(xué)成像有很多的共同點.它能夠為地球科學(xué)家們提供關(guān)于地球內(nèi)部的地震波速度、各向異性、密度乃至于地震波衰減的三維分布信息.通過研究這些三維模型以及結(jié)合巖石物理實驗和地球動力學(xué)模擬結(jié)果，我們可以進一步推測地球內(nèi)部的溫度、部分熔融、物質(zhì)組成以及水和揮發(fā)成分的分布.這些結(jié)果能夠進一步地幫助我們了解地球內(nèi)部的結(jié)構(gòu)和演化，進而更好地了解地幔對流和板塊運動.

在天然地震學(xué)中，所使用的地表記錄大多來自于天然地震，又被稱為被動震源.也可以使用人工震源激發(fā)的彈性波或者聲波來進行成像，比如能源勘探中所使用的氣槍和可控震源車、炸藥以及核爆等，這些震源被稱之為主動源.與主動震源和醫(yī)學(xué)成像不同，在被動源地震學(xué)研究當(dāng)中，由于地質(zhì)構(gòu)造以及研究區(qū)域的限制，無法自由地選擇天然震源以及接收臺站的位置.這就直接導(dǎo)致了在天然地震成像中，很難獲得均勻的成像覆蓋率.比如，目前全球地震臺網(wǎng)以及區(qū)域臺網(wǎng)主要分布在北半球的陸地上，所以相對而言，在南半球以及海洋區(qū)域的成像分辨率還比較低，這對于全球尺度的地幔成像來說是一個很大的問題.目前正在廣泛應(yīng)用和發(fā)展中的海底地震儀和海洋流動地震儀（Nolet et al.,2019），以及利用海底光纖電纜提取震動和應(yīng)變信號（Lindsey et al.,2019;Zhan et al.,2021），在未來會很好地幫助我們解決這一問題.

早期的地震層析成像主要是基于高頻射線理論和采用P 波或者S 波的走時信息.通過運動學(xué)和動力學(xué)射線追蹤，可以計算地震體波在假設(shè)的一維或三維地球模型中的走時和振幅（Cerveny,2001）.通過測量實際地震信號的走時和振幅，并且與模擬結(jié)果進行對比，試圖找到一個最優(yōu)的三維地球模型來解釋所觀測得到的地震信號.該方法可以幫助我們研究地殼、地幔以及地核的地震學(xué)結(jié)構(gòu).另外，相對于體波到時，擬合面波頻散能夠幫助我們更好地研究地殼、巖石圈以及上地幔的結(jié)構(gòu)（Barmin et al.,2001;Ekstrom et al.,1997）.另一種被用于地震成像的信號是地球的自由振蕩.通過擬合觀測和模擬的自由振蕩頻率，可以推測地幔乃至于地核的地震學(xué)參數(shù)（Dziewonski,1971;Masters et al.,1982;Ishii and Tromp,1999）.在早期的走時層析成像方法當(dāng)中，通常使用人工或者自動相位識別方法，比如長時平均/短時平均（STA/LTA）算法（Allen,1978），去測量實際的地震震相到時.這些測量方法主要提取的是高頻的到時信息，所以與高頻射線理論相匹配.目前，通過人工智能的方法也可以很好地提取高頻到時信息（Ross et al.,2018）.

隨著寬頻帶地震儀精度和穩(wěn)定性的提高，地震學(xué)家們開始使用實際觀察到的地震波波形和理論計算所得到的地震圖之間的互相關(guān)來測量特定震相的到時差.該方法所得到的到時差與信號的頻率有關(guān)，而射線理論是基于無限高頻假設(shè)，所以出現(xiàn)了理論與實際觀測之間的不匹配.從2000 年左右開始，Dahlen 等提出有限頻率層析成像（finite-frequency tomography）的概念（Dahlen et al.,2000）.這一理論結(jié)合了射線理論和散射波的波恩近似，相關(guān)的結(jié)論顯示對于一維模型，互相關(guān)所提取的到時只受到傳統(tǒng)射線周圍的速度結(jié)構(gòu)的影響，反而與射線路徑本身上的速度擾動無關(guān)，并且所影響的范圍與菲涅爾帶（Fresnel zone）的寬度有關(guān)（Hung et al.,2000;Zhou et al.,2004）.幾乎同時，Zhao 等（2000）使用自由振蕩疊加理論，同樣得到了有限頻率敏感核.相對于前者，使用自由振蕩疊加的計算量較大，但是由于使用了波動方程的嚴(yán)格解，所以它適用于一切震相，包括臨界反射/折射波，核幔邊界衍射以及面波等.這一有限頻率敏感核（finite-frequency sensitivity kernel）也被形象地稱為香蕉—甜圈敏感核（banana-doughnut kernel）.該有限頻率敏感核的概念把地震波散射融入到層析成像當(dāng)中，并考慮了波前彌合（wavefront healing）的問題（Nolet and Dahlen,2000）.在這一理論框架下，成像結(jié)果與所使用的地震波的波長有關(guān).從理論上而言，它能夠提高成像的分辨率和精確度（Montelli et al.,2004a），尤其是能夠更好地研究低速帶分布，比如地幔柱（mantle plume）等.值得注意的是，在早期的有限頻成像當(dāng)中，仍然使用高頻射線或者自由振蕩疊加理論來計算一維速度模型下的理論地震圖以及敏感核.另外，相關(guān)的敏感核在勘探地震學(xué)中也有所討論和應(yīng)用（Luo and Schuster,1991;Woodward,1992），但是直到2000 年之后這兩個方向的理論才得以統(tǒng)一.

另一方面，在勘探地震學(xué)領(lǐng)域，Lailly（1984）及Tarantola（1984）革命性地提出在聲波介質(zhì)當(dāng)中，通過直接擬合觀測和理論計算所得到的地震圖來反演地球內(nèi)部的波速結(jié)構(gòu).這一理論隨后被拓展到二維彈性/黏彈性介質(zhì)反演當(dāng)中（Gauthier et al.,1986;Mora,1987;Tarantola,1986,1988）.在這一全波形反演（full waveform inversion）理論框架下，Tarantola 使用伴隨方法（adjoint method）來計算目標(biāo)函數(shù)的梯度，然后通過局部最優(yōu)化算法，比如最速下降法或者共軛梯度法（Fletcher and Reeves,1964），來估計局部最優(yōu)的三維速度模型，進而擬合實際所觀察到的波形.值得注意的是，在使用伴隨方法來數(shù)值計算目標(biāo)函數(shù)的梯度時，需要使用逆時間傳播來計算伴隨波場.這一過程與勘探地震學(xué)中所使用的逆時偏移成像（reverse time migration）是相似的（Baysal et al.,1983;Claerbout,1971;McMechan,1983）.所不同的是在逆時偏移成像當(dāng)中，需要計算正傳以及反射波場，其中后者使用的是地表記錄所得到的首次反射波（primary reflections）波形信號，然后使用互相關(guān)成像條件（cross correlation imaging condition）來得到地殼中的反射面位置和振幅信息.在全波形反演當(dāng)中，同樣需要計算正傳波場.然而在計算反傳（伴隨）波場中，使用的是實際記錄和理論計算所得到的波形殘差來構(gòu)建伴隨波場的激發(fā)震源.可以使用多種震相，包括直達波、反射波和面波等，來構(gòu)建地球的三維結(jié)構(gòu)模型，而不僅僅像逆時偏移成像一樣來尋找地下反射界面.另外，傳統(tǒng)的逆時偏移成像不需要使用非線性迭代算法來求解一個局部最優(yōu)化問題，它可以被看成是全波形反演的首次迭代結(jié)果（Tarantola,1984）.

Tromp 等（2005）將以上的有限頻成像、伴隨方法、逆時偏移和香蕉—甜甜圈（banana-doughnut）敏感核統(tǒng)一到一個伴隨層析成像的框架中.在這一框架下，可以自由地選用各種不同的反演目標(biāo)函數(shù)，比如全波形殘差、有限頻到時或者振幅等.在伴隨成像中，可以使用相同的算法去構(gòu)建這些目標(biāo)函數(shù)的梯度.對于不同的目標(biāo)函數(shù)，只需要修改相應(yīng)的激發(fā)伴隨波場的震源時間函數(shù)即可.同時，Tromp等（2005）也提出了使用伴隨成像的框架來反演地震波衰減、各向異性以及震源參數(shù)（包括地震位置以及震源機制解）.在伴隨成像的每次迭代當(dāng)中，對于每一個主動或者被動震源，只需要計算兩次地震波場以獲得目標(biāo)函數(shù)的梯度.一次是用原有的震源信息來計算正傳波場，另一次則是使用與目標(biāo)函數(shù)有關(guān)的伴隨震源來計算反傳的伴隨波場（Liu and Tromp,2006;Tape et al.,2007），所以它的計算量與震源的數(shù)目成正比，而與所使用的臺站和震相的數(shù)量無關(guān).這種方法可以幫助我們有效地計算目標(biāo)函數(shù)的梯度.但是仍然無法得到相應(yīng)的二階導(dǎo)數(shù)，也就是Hessian 矩陣的信息（Plessix,2006）.另一種相似的方法則是由Zhao 等（2005）提出的散射積分方法（scattering-integral）.利用地震波動方程格林函數(shù)所具有的源與觀測點之間的互易性，通過使用硬盤存儲接收點格林函數(shù)，不僅可以計算每個源和臺站對之間的敏感核，而且能夠同時構(gòu)建目標(biāo)函數(shù)的梯度及其二階導(dǎo)數(shù).這樣能夠幫助我們分析反演結(jié)果的分辨率和不確定性信息.但是由于需要存儲四維地震波場庫，而且每一次迭代時都需要重新計算所有接收點格林函數(shù)，所以這種方法相對于伴隨成像需要更大的存儲空間（Chen et al.,2007a）.而在有些特定情況下，其計算量可能與伴隨成像接近.另外，由于伴隨層析成像和全波形反演具有許多相似性.所以在本文中，我們有時會相互替代使用這兩種名稱.全波形反演理論上需要使用記錄的全部波形信號.但是考慮到實際記錄當(dāng)中的噪聲成分以及反演當(dāng)中的非線性，很少有研究將全部記錄所得的波形直接應(yīng)用于反演當(dāng)中.

自從全波形反演理論的提出，這種新方法已經(jīng)被應(yīng)用到各種尺度的地球內(nèi)部結(jié)構(gòu)成像研究當(dāng)中.比如，Tape 等和Chen 等使用伴隨/全波形成像來研究南加州地殼的三維速度結(jié)構(gòu)（Chen et al.,2007b;Lee et al.,2014;Tape et al.,2009,2010）.他們的模型能夠幫助我們很好地研究這一區(qū)域的地質(zhì)構(gòu)造，并且已被美國南加州地震中心（Southern California Earthquake Center）作為區(qū)域標(biāo)準(zhǔn)的三維公共速度模型（比如最新的CVM-H15.1.0 和CVMS4.26）.同時這些高分辨率高精度的三維速度模型能夠幫助我們更好地研究南加州地震的震源信息（Lee et al.,2011;Liu et al.,2004;Wang and Zhan,2020;Zhao et al.,2006）.另一方面，全波形反演方法也被廣泛地應(yīng)用到區(qū)域上地幔成像，包括歐洲（Fichtner et al.,2013;Fichtner and Villasenor,2015;Zhu et al.,2012,2015;）、美洲（Krischer et al.,2018;Zhu et al.,2017,2020b）、大西洋（Colli et al.,2013）、亞洲（Chen et al.,2015;Simut? et al.,2016;Tao et al.,2018）、澳大利亞（Fichtner et al.,2009,2010）、新西蘭（Chow et al.,2022）、南極洲（Lioyd et al.,2019）等.目前也被應(yīng)用于全球地幔研究 當(dāng)中（Bozda? et al.,2016;French et al.,2013;French and Romanowicz,2015;Lei et al.,2020）.同時，全波形反演也在油氣資源開發(fā)過程中發(fā)揮越來越重要的作用（Operto et al.,2015;楊積忠等,2014）.我們將會在2.8—2.10 節(jié)中詳細討論這些具體的例子.迄今為止，已經(jīng)有一些全波形反演的開源軟件可以供研究人員下載和使用，比如Seisflow（https://github.com/rmodrak/seisflows）、LASIF（https://dirkphilip.github.io/LASIF_2.0/）、PySIT（https://pysit.readthedocs.io/）和IFOS3D（https://git.scc.kit.edu/GPIAG-Software/IFOS3D/）等.

1 理論框架

本文只是回顧全波形反演的主要理論框架和計算過程，具體的推導(dǎo)可以參考其他文獻和專著（Chen and Lee,2015;Fichtner,2010;蔣夢凡等,2021;Liu and Gu,2012;Tromp et al.,2005;Virieux and Operto,2009;楊午陽等,2013;楊勤勇等,2014）.

1.1 目標(biāo)函數(shù)

全波形反演在數(shù)學(xué)上是求解一個非線性最優(yōu)化問題（Liu and Gu,2012;Plessix,2006;Tromp et al.,2005）.由于所涉及的模型參數(shù)量以及正演數(shù)值計算量十分巨大，目前在實際應(yīng)用當(dāng)中，只能通過使用目標(biāo)函數(shù)的導(dǎo)數(shù)與局部最優(yōu)化算法來求解這一問題.常用的方法有共軛梯度（Fletcher and Reeves,1964）或者L-BFGS 算法（Nocedal and Wright,2006）.在目前階段尚無法使用全局優(yōu)化算法（Sen and Stoffa,2013;Tarantola,2005），比如模擬退火法或者蒙特卡羅法，來解決大型的三維反演問題.這一局限性致使目前大部分研究只能給出某一個問題的局部最優(yōu)解，而無法準(zhǔn)確地提供關(guān)于這一問題的后驗概率分布（posteriori probability density distribution）以及不確定性分析（uncertainty quantification），這仍然是當(dāng)前研究的前沿問題（詳見2.5 節(jié)）.

由于以上的問題，所以選取一個最佳的目標(biāo)函數(shù)/殘差函數(shù)來求解這一局部最優(yōu)化問題至關(guān)重要.這一目標(biāo)函數(shù)是我們想要最小化的標(biāo)量，它也被用來衡量最優(yōu)化過程是否收斂.地震學(xué)家們主要使用的觀測數(shù)據(jù)是地表位移記錄，也就是地震圖.它本質(zhì)上是記錄地表震動的時間和空間的函數(shù).目前大部分區(qū)域尚沒有密集的臺站，只有少數(shù)區(qū)域具有覆蓋面積大且分布密集的臺網(wǎng)，比如USArray 和Hinet 等.所以不失一般性，我們假設(shè)所采集的地震圖是一維的時間函數(shù).使用目標(biāo)函數(shù)來衡量實際觀測與理論模擬的地震記錄之間的差別.在伴隨層析成像和全波形反演中，最直接也是最早使用的目標(biāo)函數(shù)是波形之間的最小二乘殘差，即：

式中，目標(biāo)函數(shù)J是一個標(biāo)量.在公式（1）中它的單位是 m2.sw和dw分別表示模擬和觀測所得到的某個分量w上的地震圖.e和r分別表示某一震源和臺站.Ns、Nr、Nw分別表示地震、臺站和時間窗的數(shù)目.由于觀測所得到的地震圖中難免有噪聲，所以我們可以使用特定的時間窗來測量某一個時間段中的波形差（Maggi et al.,2009）.時間窗的選取可以根據(jù)波形之間的相似度以及所需要研究的特定震相來決定.其中值得注意的是，時間窗中的信號可以是單一的傳統(tǒng)震相，也可以是多個不同路徑的震相疊加.這一優(yōu)勢使得全波形反演和伴隨成像能夠更加充分地利用地震記錄信息來約束地球內(nèi)部結(jié)構(gòu).由于陸地地震臺可以記錄三分量地表震動，所以目標(biāo)函數(shù)可以根據(jù)三分量位移來計算求得.

全波形目標(biāo)函數(shù)的好處在于它能夠最直接地衡量觀測與理論計算之間的細微差別.所以它能夠幫助我們更好地約束地球模型，以及提供更高分辨率的反演結(jié)果（Virieux and Operto,2009）.但是它的缺點包括：（1）如果記錄地震圖的信噪比較低，則反演過程主要是去擬合不必要的噪聲，這些噪聲最終會被投射到反演模型當(dāng)中，導(dǎo)致錯誤的結(jié)果和解釋；（2）如果觀測和模擬波形之間的到時差大于信號的半周期，則會出現(xiàn)所謂的周期跳躍（cycle skipping）現(xiàn)象（Virieux and Operto,2009）.從而使得全波形目標(biāo)函數(shù)中出現(xiàn)許多局部最小值（董良國等,2013;Engquist and Yang,2019;Métivier et al.,2016;Yang et al.,2018）.因為目前我們大多數(shù)時候使用的是目標(biāo)函數(shù)的梯度來求解這一局部最優(yōu)化問題，所以如果起始模型遠離全局最優(yōu)解，只能獲得起始模型周圍的局部最優(yōu)解，而非全局最優(yōu)解.這一問題在高頻全波形反演中尤為嚴(yán)重.

考慮到以上全波形目標(biāo)函數(shù)所存在的諸多問題，另一種被廣泛使用的目標(biāo)函數(shù)是基于到時信息.這一點與傳統(tǒng)的P 波或S 波射線走時層析成像的概念相同.可以設(shè)計如下的目標(biāo)函數(shù)：

式中， ΔTw(e,r;m)是觀測以及預(yù)測的某個震相的到時差.相對于全波形最小二乘目標(biāo)函數(shù)[公式（1）]，使用到時殘差目標(biāo)函數(shù)可以更好地幫助我們在反演過程中避免陷入局部最小值（Luo and Schuster,1991）.

另一類目標(biāo)函數(shù)被稱之為雙差函數(shù)，它使用兩個臺站接收到的同一震源的到時差來反演地球結(jié)構(gòu).這與傳統(tǒng)的雙差層析成像（Zhang and Thurber,2003）和地震定位（Waldhauser and Ellsworth,2000）的概念是相似的，該方法目前也被用于伴隨層析成像當(dāng)中（Yuan et al.,2016）.它可以減少震源對地球結(jié)構(gòu)反演的影響，同時主要反映臺站之間的速度結(jié)構(gòu).

1.2 模型參數(shù)化

在設(shè)立了目標(biāo)函數(shù)之后，第二個問題就是要確定我們所要求解問題的模型參數(shù).它涉及到以下幾個方面的問題：

（1）如何離散化所關(guān)心的研究目標(biāo).其中最直接的方法就是使用正演模擬當(dāng)中的數(shù)值網(wǎng)格來參數(shù)化連續(xù)性的地球模型，這也是目前大多數(shù)研究中所采用的方法.它的優(yōu)點是正演和反演使用的是同一套網(wǎng)格點，易于操作.缺點則是在正演模擬當(dāng)中，為了滿足數(shù)值計算的精度，穩(wěn)定性以及頻散要求（Komatitsch and Tromp,1999;Virieux et al.,2011），網(wǎng)格點的數(shù)目一般十分巨大.比如，在地殼或地幔反演當(dāng)中，所使用的正演網(wǎng)格點可以輕松到達百萬級別（Tape et al.,2010;Zhu et al.,2015）.另一種方法就是用全球?qū)游龀上裰兴鶑V泛使用的基函數(shù)來離散化所要研究的目標(biāo)（Nolet,2011;Rawlinson et al.,2010）.常用的基函數(shù)包括水平方向上的球諧函數(shù)以及深度方向上的樣條函數(shù)（Dziewonski and Romanowicz,2015），使用這些基函數(shù)可以大大地減少模型參數(shù)量.另一個優(yōu)點則是當(dāng)?shù)厍蚰Ｐ捅旧砭哂杏邢薜牟〝?shù)分布時，這樣也可以大大地減少模型參數(shù)量.這點在全球地幔成像當(dāng)中尤為突出，比如，地幔的主要結(jié)構(gòu)用6～8 階的球諧函數(shù)就可以很好地描述（Su and Dziewonski,1992）.另外使用基函數(shù)還可以起到對反演模型進行平滑或者是低通濾波的效果.

（2）傳統(tǒng)地震成像最關(guān)心的地球模型參數(shù)是各向同性的P 波和S 波速度（Grand et al.,1997;Van der Hilst et al.,1997;Zhao,2004）.通過反演求得的這些波速信息可以幫助我們分析地球內(nèi)部的溫度、物質(zhì)組成以及部分熔融和揮發(fā)成分等.我們可以通過各種方法來構(gòu)造目標(biāo)函數(shù)與模型參數(shù)之間的關(guān)系如下：

式中，α、β 和ρ 分別表示P 波和S 波的速度以及密度.Kα、Kβ和Kρ分別是目標(biāo)函數(shù)相對于以上三個模型參數(shù)的梯度.如果我們使用到時差[公式（2）]作為目標(biāo)函數(shù)，同時用高頻射線或者有限頻方法來計算這些導(dǎo)數(shù)，則以上的公式與傳統(tǒng)的射線走時層析成像是一致的.與全波形反演所不同的是，我們使用三維地球模型以及數(shù)值伴隨算法來計算上述目標(biāo)函數(shù)的梯度（見1.3 節(jié)）.

使用伴隨方法的另一個好處是可以考慮多個不同的模型參數(shù)，包括各向異性波速度以及地震波衰減.即可以在一個統(tǒng)一的理論框架下來求解地球的多參數(shù)模型（Fichtner and Driel,2014;Sieminski et al.,2007a,2007b;Tromp et al.,2005）.這也是目前的一個研究重點，本文將會在2.2、2.3 和2.6 節(jié)中具體討論.

1.3 使用伴隨方法計算目標(biāo)函數(shù)的梯度

從公式（3）中可以看到，伴隨成像/全波形反演與高頻射線和有限頻成像的一個重要區(qū)別在于如何計算目標(biāo)函數(shù)的梯度.在伴隨層析成像和全波形反演當(dāng)中，使用的是伴隨方法來計算這一三維空間函數(shù)（Plessix,2006）.根據(jù)不同的目標(biāo)函數(shù)，這些梯度的單位也有所不同.如果目標(biāo)函數(shù)是全波形最小二乘殘差，梯度的單位是而如果我們選取到時殘差作為目標(biāo)函數(shù)，那么它們的單位則是同時值得注意的是，這里的目標(biāo)函數(shù)的梯度與前面所描述的敏感核是不一樣的，后者不包括實際觀測的數(shù)據(jù)和殘差.使用伴隨方法來推導(dǎo)目標(biāo)函數(shù)的梯度大體上可以分為以下兩類：一種是基于地震波場的波恩近似展開來推導(dǎo)（Tarantola,1984;Tromp et al.,2005）.另一種則是使用拉格朗日乘數(shù)（Lagrange multiplier）的方法來推導(dǎo)（Liu and Tromp,2006;Virieux and Operto,2009）.以上兩種方法所得到的結(jié)果是相同的.具體的數(shù)學(xué)推導(dǎo)過程可以參考相關(guān)的文獻和專著（Chen and Lee,2015;Fichtner,2010）.

基于Tromp 等（2005）的推導(dǎo)，可以使用正傳和伴隨波場的卷積來構(gòu)建目標(biāo)函數(shù)的梯度，其中伴隨波場可以通過計算反傳波場得到（如圖1 所示）.用包含測量信息的伴隨時間函數(shù)在接收臺站上作為震源來激發(fā)伴隨波場.所以我們可以使用正演的數(shù)值算法和程序來求解伴隨波場.對于二階波動方程，正傳和伴隨波場之間的差別只在于波傳播的時間方向以及震源時間函數(shù).根據(jù)不同的模型參數(shù)，可以使用正傳和伴隨波場的時間和空間導(dǎo)數(shù)來計算它們的梯度.如果選取密度ρ、體積模量 κ和剪切模量 μ作為模型參數(shù)，它們的梯度可以通過以下的公式求得：

圖1 二維均勻速度模型下構(gòu)建SH 波的敏感核.從左到右分別是正傳波場、伴隨波場、相互作用波場和剪切波速度的敏感核.五角星和方塊分別表示震源和接收臺站的位置.時間從下往上分別是8 s、16 s、24 s、32 s 和44 s.地表使用的是自由表面邊界條件（修改自Tromp et al.,2005 中的圖3）Fig.1 Construction of an SH sensitivity kernel in a 2D homogeneous velocity model.From left to right are forward wavefield,adjoint wavefield,interaction wavefield and shear wave velocity sensitivity kernel.The star and rectangular represent source and receiver.The time steps from bottom to top are 8,16,24,32 and 44 seconds.A traction free boundary condition is applied to the Earth's surface (modified from Figure 3 in Tromp et al.,2005)

式中，s和s+分別是正傳和伴隨波場，D和D+則是正傳和伴隨應(yīng)變偏量（strain deviator），T是記錄時間.

如果我們選取波阻抗Zα，以及P 波和S 波速度作為模型參數(shù)，則可以通過以下的公式來計算它們的梯度：

值得注意的兩點是：（1）在有的文獻中可以看到使用的是卷積[如公式（4）]，而在另一些論文中使用的則是互相關(guān)（主要在勘探地震學(xué)中）.它們本質(zhì)上是一樣的，區(qū)別只在于對伴隨波場的時間傳播方向的定義不同.（2）公式（4）中的密度梯度與勘探地震學(xué)中的逆時偏移成像非常相似，這一點在Tarantola（1984）中已經(jīng)指出.Claerbout（1971）得到的成像條件是基于物理的直觀考慮，就是地下反射面存在于下行（正傳）波場與反射（伴隨）波場相位一致的時候.在這一點上數(shù)學(xué)推導(dǎo)與物理直觀達到了高度統(tǒng)一.只不過Claerbout 的成像條件是目標(biāo)函數(shù)相對于密度的梯度.而后通過數(shù)值計算，Zhu 等（2009）發(fā)現(xiàn)相對于密度梯度而言，波阻抗梯度能夠更好地反映模型的高波數(shù)信息（Wu and Aki,1985），也就是反射率.這一點已被工業(yè)界用于實際勘探開發(fā)當(dāng)中，被稱為逆反射成像（inverse scattering）（Douma et al.,2010）或者是能量模式成像條件（Rocha et al.,2016）.

1.4 最優(yōu)化迭代

通過以上的伴隨方法可以求得相關(guān)模型的目標(biāo)函數(shù)梯度.然后通過基于導(dǎo)數(shù)的局部最優(yōu)算法，比如共軛梯度或者是L-BFGS 方法，可以獲得初始模型附近的局部最優(yōu)解（Akcelik et al.,2002,2003）.其中，共軛梯度法使用的是本次迭代的梯度以及上次迭代的方向之間的加權(quán)平均來求得本次迭代的方向.相對于最簡單的最速下降法，它具有更好的收斂性.另一類方法被稱為擬牛頓迭代法，其中一個典型代表是L-BFGS 方法（Nocedal and Wright,2006）.由于使用了前幾次迭代的方向信息，LBFGS 方法能夠幫助我們進一步構(gòu)建近似的二階導(dǎo)數(shù)信息.所以相對于共軛梯度和最速下降算法，它具有更好的收斂性（圖2），進而加快整個反演的進度（Modrak and Tromp,2016）.這一點在實際應(yīng)用中至關(guān)重要，因為相對于射線和有限頻成像，伴隨層析成像和全波形反演的計算量要求仍然較高.所以如果能夠減少迭代的次數(shù)，則能更加有效地反演地下模型，進而更好地研究相關(guān)的科學(xué)問題.

圖2 對比九個不同理論模型下使用三種迭代方法所得到的收斂性，包括最速下降法（黑線），非線性共軛梯度法（綠線）和L-BFGS 法（紅線）.不同的理論模型如每個圖的標(biāo)題所示（修改自 Modrak and Tromp,2016 中的圖2）Fig.2 Convergency comparison of three iterative methods for nine synthetic velocity models,including the steepest descent method(black lines),nonlinear conjugate gradient method (green lines) and L-BFGS method (red lines).The 2D synthetic velocity models are shown in the titles of each panel (modified from Figure 2 in Modrak and Tromp,2016)

另外兩種可以加速迭代收斂的方法是：（1）使用近似方法來構(gòu)建Hessian 矩陣.因為Hessian 矩陣本身十分巨大，它的大小與模型的網(wǎng)格點數(shù)的平方成正比（Pratt et al.,1998）.所以目前仍然很難直接地計算這一矩陣.但是通過各種近似，可以構(gòu)建它的對角量，使用該信息可以幫助我們更好地平衡模型淺部與深部的梯度大小，并且可以減少在震源以及接收臺站上模型梯度數(shù)值較大的問題.（2）使用高斯-牛頓方法來迭代地考慮Hessian 矩陣的作用（Epanomeritakis et al.,2008;劉璐等,2013;Métivier et al.,2014），這被稱之為截斷高斯牛頓（Truncated Gauss-Newton）方法.在這一方法中，在外部共軛梯度或者L-BFGS 循環(huán)當(dāng)中加入一個內(nèi)循環(huán)，通過求解法方程式（Normal equation）來迭代地考慮Hessian 矩陣的作用.這一方法已被用于勘探地震學(xué)的數(shù)值試驗當(dāng)中，但是相對于一般的共軛梯度或者L-BFGS 方法，它涉及到內(nèi)外兩個迭代循環(huán)，所以計算量仍然十分巨大.

通過以上的迭代，如果目標(biāo)函數(shù)的數(shù)值不斷減小，而且迭代模型趨于穩(wěn)定，則在一定迭代次數(shù)之后，我們可以得到一個最終優(yōu)化的模型.另外，在迭代過程中所廣泛使用的技巧是多尺度反演（Bunks et al.,1995;王毓瑋等,2016;張文生等,2015）.也就是在迭代的早期，我們使用長周期的信號來反演大尺度的結(jié)構(gòu).隨著反演的推進，我們逐漸地加入短周期的信號來進一步約束小尺度的結(jié)構(gòu).使用這一方法能夠更好地幫助我們規(guī)避目標(biāo)函數(shù)在初始模型周圍的局部最小值以及反演當(dāng)中的強非線性問題.

2 研究進展

2.1 目標(biāo)函數(shù)的設(shè)計

如1.1 節(jié)所述，設(shè)計一個好的目標(biāo)函數(shù)對于全波形反演是否成功具有十分重要的意義.而波形最小二乘殘差具有前述的諸多缺點，所以如何設(shè)計一個好的目標(biāo)函數(shù)是目前的一個研究熱點.本文總結(jié)幾個大的目標(biāo)函數(shù)類型.

（1）直接基于時間信號的特征.這一類型包括前面的波形最小二乘差（Tarantola,1984）以及后來提出的波形包絡(luò)差（胡勇等,2017;Wu et al.,2014）.相關(guān)的研究顯示使用波形包絡(luò)可以很好地增加信號中的低頻信息，進而減小反演問題的非線性以及避免陷入局部最小區(qū)域；（2）從地震圖中提取到時和頻率的信息.這一類包括使用互相關(guān)來提取觀測與預(yù)測的時間差（Leeuwen and Mulder,2010;Luo et al.,2016）.可以直接在時間域來測量，也可以通過小波變換，進而在時間—頻率域中測量（Fichtner et al.,2008;Yuan and Simons,2014）.該方法也可以用來測量振幅差；（3）自適應(yīng)反演.Warner 等提出使用維納濾波（Wiener filter）的概念來測量兩個波形之間的差（Warner and Guasch,2016;Zhu and Fomel,2016），這一方法可以不使用時間窗來分割地震信號.相關(guān)的數(shù)值試驗顯示很好的反演效果.其他類似的方法包括動態(tài)擬合（dynamic warping）（Hale,2013;Ma and Hale,2013）和地震圖匹配（seismogram registration）（Baek et al.,2013;Zhu,2018a）；（4）最優(yōu)傳輸方法.這一方法是基于數(shù)學(xué)中的最優(yōu)傳輸路徑來測量兩個時間序列之間的差別（Engquist and Yang,2019;Métivier et al.,2016;Yang et al.,2018）.相對于最小二乘波形差，這一目標(biāo)函數(shù)具有更寬的收斂帶，從而對初始模型的要求更低（圖3）.對于一維時間序列，有很好的數(shù)學(xué)方法來求解這一最優(yōu)化問題，進而得到相關(guān)的伴隨震源信號.相關(guān)的研究也被應(yīng)用于實際地殼反演當(dāng)中（Dong et al.,2022）.但是目前對于二維圖像的對比，需要求解Monge-Ampere 方程（Engquist and Froese,2014）.這是一個非線性方程組，它的穩(wěn)定數(shù)值解還有待進一步研究.同時這一方法要求兩個信號之間具有相同的能量，這在實際的地震記錄中是很難得到滿足的.其他在勘探地震學(xué)中開發(fā)的用于解決周期跳躍的方法還有時間軸擴展（Biondi and Almomin,2014）和波場重構(gòu)（van Leeuwen and Herrmann,2013）等.

圖3 對比不同目標(biāo)函數(shù)的特征.（a）顯示的是子波波形.（b）和（c）分別顯示的是基于最小二乘波形殘差和最優(yōu)化路徑所得到的目標(biāo)函數(shù)隨著不同時間移動（s）的特征（修改自Engquist and Froese,2014 中的圖1）Fig.3 Comparison of misfit functions based on least-square waveform differences and optimal transport distances.Panel (a) shows the input source wavelet.Panels (b) and(c) are the misfits as functions of time shifts for the leastsquares waveform differences and optimal distances,respectively (modified from Figure 1 in Engquist and Froese,2014)

2.2 各向異性反演

我們使用各向異性來描述地震波沿不同方向所具有的不同傳播速度，在地震學(xué)研究中包括徑向各向異性和方位各向異性（Backus,1965;Hess,1964;Montagner and Nataf,1986;Smith and Dahlen,1973）.這些各向異性來自于形狀或是晶格的定向排列[shape-preferred orientation (SPO)/lattice-preferred orientation (LPO)]（Jung and Karato,2001;Karato et al.,2008;Nicolas and Christensen,1987;Zhang and Karato,1995）.地殼和地幔中的各向異性礦物包括橄欖石、輝石以及云母，它們是導(dǎo)致地震波各向異性的主要原因.目前可以使用剪切波分裂（Silver and Chan,1991;Silver,1996）以及面波頻散（Lin et al.,2011;Marone and Romanowicz,2007;Moschetti et al.,2010）來研究地球內(nèi)部的各向異性分布.相對來說，如果有很好的臺站分布，橫波分裂可以提供很好的橫向分辨率.但是由于使用的是垂向入射的遠震體波信號，所以對于垂向上的各向異性分布具有很差的分辨率.相反，因為不同頻率的面波對不同深度的速度具有不同的敏感性,使用面波頻散可以幫助我們提高垂向上的分辨率.但是由于在面波相速度或群速度反演當(dāng)中，使用最小二乘法來求解一個層析成像問題，通常需要加入正則化約束.同時由于面波的橫向傳播，相對來說，它的橫向分辨率要低于垂向分辨率.通過研究地球內(nèi)部的各向異性，可以討論地殼、巖石圈與軟流圈中的應(yīng)力、形變以及流動性特征，同時也可以用來研究地幔中的對流以及地核構(gòu)造（Long and Silver,2008;Long and Becker,2010;Park and Levin,2002）.

如前所述，伴隨層析成像和全波形反演提供了一個統(tǒng)一的框架來反演各種不同模型參數(shù)的組合.所以它也可以用于反演各向異性（Rusmanugroho et al.,2017;Sieminski et al.,2007a,2007b）.比如在地殼成像方面，Wang 等（2020）在Tape 等的南加州地殼模型基礎(chǔ)之上，加入了徑向各向異性.他們發(fā)現(xiàn)在圣安德烈斯斷層西面的中下地殼中，橫向偏振的剪切波速度比縱向偏振要低（負徑向各向異性～-6%）.而跨過斷層，這一比例正好相反.這與上地殼和上地幔中常見的徑向各向異性特征正好相反.為了解釋這一現(xiàn)象，他們提出一種可能性是由于下地殼的斜長石具有與黑云母和橄欖石不一樣的各向異性特征（Wang et al.,2020），導(dǎo)致它們的快速軸與應(yīng)變方向垂直而非平行.在地幔成像方面，Zhu 等（2020c）使用伴隨方法反演上地幔的方位各向異性，進而研究俯沖板塊與地幔對流的關(guān)系.例如，他們發(fā)現(xiàn)在中南美洲和卡斯凱迪亞俯沖帶中，當(dāng)出現(xiàn)板塊斷裂時，上地幔物質(zhì)會流動通過這些板塊的間斷空間（圖4），從而在板塊窗以及板塊周圍產(chǎn)生環(huán)流（Hall et al.,2000;Russo and Silver,1994）.這些觀測與巖石物理實驗（Buttles and Olson,1998;Kincaid and Griffiths,2003）和地球動力學(xué)模擬相吻合（Schellart,2004;Schellart et al.,2007;Stegman et al.,2006）.另外，在歐洲上地幔反演結(jié)果當(dāng)中（Zhu and Tromp,2013），發(fā)現(xiàn)反演所得到的方位各向異性方向與過去三千萬年的大地構(gòu)造過程有很好的相關(guān)性.同時通過分析方位各向異性的大小在垂向上的變化，他們推測出脆性和韌性變形在巖石圈中的分布，這一結(jié)果與實驗和震源深度分布結(jié)果相吻合（Chen and Molnar,1983;Kohlstedt et al.,1995）.另外，通過同時反演Rayleigh 波和Love 波，F(xiàn)ichtner 等（2010）發(fā)現(xiàn)澳大利亞大陸和海洋巖石圈在徑向各向異性上存在差別.比如在海洋巖石圈中，正徑向各向異性（橫向偏振剪切波快于縱向偏振剪切波）主要存在于80～150 km 之間.而在大陸巖石圈中，徑向各向異性比較弱且主要出現(xiàn)在小于100 km 左右的區(qū)域.這與早期的全球面波成像結(jié)果相吻合（Nettles and Dziewonski,2008）.

圖4 三維俯沖板塊導(dǎo)致的地幔對流.（a）和（b）分別顯示中美洲和卡斯凱迪亞俯沖帶的反演結(jié)果.綠色塊體顯示的是剪切波速度擾動大于1.5%的結(jié)果.黃色箭頭表示通過伴隨方位各向異性成像所得到的地幔對流情況[（a）和（b）分別修改自 Zhu et al.,2020a 中的圖7 和 Zhu et al.,2020c 中的圖5 ]Fig.4 Three dimensional subducting slabs and induced mantle flows.Panels (a) and (b) are results for the Middle American and Cascadian subduction zones.Green bodies represent regions with shear wave velocity perturbations greater than 1.5%.Yellow arrows represent horizontal mantle flows constrained by azimuthal anisotropy adjoint tomography [Panels (a) and (b) are modified from Figure 7 in Zhu et al.,2020a and Figure 5 in Zhu et al.,2020c,respectively]

2.3 地震波衰減成像

地震波在地球內(nèi)部傳播過程中經(jīng)受衰減作用，其中動能和彈性勢能被轉(zhuǎn)換成熱能.相對于過去30多年地震波波速層析成像的顯著進步來說（Grand et al.,1997;Meschede and Romanowicz,2015;Schaeffer and Lebedev,2013），地震波衰減成像的進步相對較?。―alton et al.,2008;Gung and Romanowicz,2004;Romanowicz,1995）.相對于以前的結(jié)果，在最近的一些大陸和全球尺度的衰減成像當(dāng)中（圖5），能夠看到相對較好的可比性（Adenis et al.,2017a,2017b;Bao et al.,2016a,2016b;Ma et al.,2016;Karaoglu and Romanowicz,2017）.相對于波速成像，地震波衰減成像有以下的幾個難點：（1）需要有高質(zhì)量的地震波振幅信息，這就需要具有較高信噪比的地震記錄，同時要對地震儀器響應(yīng)進行校正；（2）體波走時和面波頻散主要受地震波波速的影響，然而地震波振幅受到多重因素的影響，其中包括介質(zhì)衰減、波速、震源大小、輻射模式、彈性聚焦與散焦（Durek et al.,1993;Romanowicz,1987;Ruan and Zhou,2010,2012）、散射以及臺站局部的放大效應(yīng)等（Eddy and Ekstrom,2014;Lin et al.,2012）；（3）目前尚難以區(qū)分內(nèi)在衰減和由地震波散射所導(dǎo)致的表象衰減.所以在衰減成像當(dāng)中，需要仔細地考慮這些因素.

目前，全波形反演也被用于研究地球內(nèi)部地震波衰減的分布（Karaoglu and Romanowicz,2017,2018;Zhu et al.,2013）.Tromp 等揭示了在伴隨成像/全波形反演框架下，如果使用標(biāo)準(zhǔn)線性固體模型（standard linear solid model）（Liu et al.,1976）來表示一定頻段下與頻率無關(guān)的品質(zhì)因子Q，那么Q的目標(biāo)函數(shù)的梯度計算公式與剪切模量的公式是一樣的.只需要把計算伴隨波場的震源時間函數(shù)改成新的計算公式即可（Tromp et al.,2005）.所以，在這一框架下，相對于波速成像，需要使用兩個不同的伴隨震源時間函數(shù)來計算伴隨波場，從而同時得到目標(biāo)函數(shù)對于波速和衰減的梯度信息，這也就使得反演的計算量翻倍.這一方法被Zhu 等（2013）應(yīng)用到歐洲地幔反演當(dāng)中.他們發(fā)現(xiàn)在上地幔淺部，品質(zhì)因子Q模型與地震波波速分布具有較好的相關(guān)性，這一點與全球尺度的地震波衰減成像結(jié)果相吻合（圖5）.比如，洋中脊地區(qū)具有較高的衰減和低速度帶，而在大陸克拉通下則是低衰減和高波速帶.另外，在北大西洋的地幔轉(zhuǎn)換帶當(dāng)中，他們發(fā)現(xiàn)Q值比周圍較小，表示這一區(qū)域的衰減較大.為了解釋這一現(xiàn)象，他們認為這一區(qū)域應(yīng)受含水量的影響較大，這與地幔轉(zhuǎn)換帶水過濾的科學(xué)假設(shè)相吻合（Bercovici and Karato,2003;Huang et al.,2005;Schmandt et al.,2014）.另一方面，Komatitsch 等（2016）改進了伴隨波場的存儲方法，從而大大地改進了計算衰減敏感核的穩(wěn)定性和精度問題.

圖5 對比三個全球尺度的上地幔地震波衰減成像結(jié)果.（a-c）分別來自于QRLW8（Gung and Romanowicz,2004）、SEMUCB-UMQ（Karaoglu and Romanowicz,2018）和QRFSI12（Dalton et al.,2008）.每一個圖的下方第一行和第二行分別顯示的是各個深度的剪切模量衰減值及其擾動量，其中SEMUCB-UMQ 是全波形衰減反演的結(jié)果（修改自Karaoglu and Romanowicz,2018 中的圖9 ）Fig.5 Comparison of three global scale upper mantle seismic attenuation tomography models.Panels (a-c) are results from QRLW8(Gung and Romanowicz,2004)，SEMUCB-UMQ (Karaoglu and Romanowicz,2018) and QRFSI12 (Dalton et al.,2008).The first and second lines in each panel represent absolute and relative perturbations of shear modulus attenuation (modified from Figure 9 in Karaoglu and Romanowicz,2018)

標(biāo)準(zhǔn)線性固體模型只是眾多表示地震波衰減理論模型當(dāng)中的一個，它具有以下幾個缺點：（1）品質(zhì)因子Q并不在所推導(dǎo)的地震波公式當(dāng)中，相反，需要求得相應(yīng)的應(yīng)力和應(yīng)變松弛時間因子（Carcione et al.,1988;Robertsson et al.,1994）.這給反演和正演都帶來了一定的困難；（2）如上所述，如果要同時獲得波速和衰減的梯度信息，需要計算兩次伴隨波場，從而加大了反演的計算量.為了解決上述問題，在過去幾年當(dāng)中，一些用以描述地震波在地球內(nèi)部衰減現(xiàn)象的新的理論模型被提出.比如Fichtner 和Driel（2014）提出一個改進的標(biāo)準(zhǔn)線性固體模型，使得Q直接出現(xiàn)在波場公式中，進而簡化了計算目標(biāo)函數(shù)相對于Q的梯度的過程.同時這種方法也可以幫助我們表示與頻率相關(guān)的Q的情況（Lekic et al.,2009）.另一方面，Zhu 等和Yang 等從不同的角度提出了常量Q的聲波方程（Yang and Zhu,2018;Zhu and Harris,2014）.而后這些方程被推廣到黏彈性介質(zhì)（Zhu and Sun,2017），以及品質(zhì)因子Q的全波形反演過程當(dāng)中（Yang et al.,2020）.

另外在地震波衰減反演當(dāng)中，選取好的反演策略和目標(biāo)函數(shù)也十分重要.相關(guān)的研究表明使用兩個反演階段具有更好的效果，即第一階段主要使用到時信息來約束地震波速度，而在第二階段中加入振幅并同時反演速度和衰減（Kamei and Pratt,2013;Zhu et al.,2013）.另外，相對于波形和到時目標(biāo)函數(shù)，基于振幅的目標(biāo)函數(shù)，比如波形包絡(luò)差和振幅頻譜比，對衰減反演具有更好的效果（Pan and Wang,2020）.

2.4 模型正則化

模型正則化是求解反演問題當(dāng)中一個必不可少的重要環(huán)節(jié)，它在射線和有限頻成像中都起到平滑模型以及提供先驗信息的作用（Aster et al.,2018;Backus and Gilbert,1968,1970;Tarantola,2005）.同樣在全波形反演當(dāng)中，正則化也被考慮進來.因為和所有其他的反演問題一樣，我們也要面對問題的非唯一性，數(shù)據(jù)的不準(zhǔn)確性以及采樣的不完備性.目前主要有以下幾種途徑來提供正則化約束：（1）對目標(biāo)函數(shù)相對于模型參數(shù)的梯度進行高斯平滑或者低通濾波.這樣做的效果與使用Tikhonov正則化是一致的（Operto et al.,2006;Tape et al.,2007）；（2）顯性地使用Tikhonov正則化.這樣能夠幫助我們在擬合實際觀測數(shù)據(jù)的同時，提供一個光滑的或者是具有最小能量/振幅的模型；（3）可以在正則化約束當(dāng)中提供相關(guān)的先驗信息，比如Asnaashari 等（2013）將測井所得到的一維速度模型考慮到正則化約束當(dāng)中，這樣能夠幫助我們在擬合地表觀測地震圖的同時也擬合實際得到的測井?dāng)?shù)據(jù)；（4）使用稀疏正則化.相對于Tikhonov 正則化，如果我們的模型在物理空間或者某個特定轉(zhuǎn)換空間下具有系數(shù)稀疏性的特征.那么可以在反演過程中構(gòu)建一個稀疏化約束，從而使用更少的模型數(shù)量來擬合實際觀測數(shù)據(jù)（Lin and Huang,2014）.如圖6 所示，使用改進的全變差正則化方法所得到的彈性波全波形反演模型能夠更好地約束強速度間斷面，同時具有提高模型信噪比的作用.另外，可以使用小波變換或seislet 來將模型轉(zhuǎn)換到相應(yīng)的變換空間，然后用全變差（total variation）來做相關(guān)的稀疏性約束（Askan and Bielak,2008;Xue et al.,2017）.

圖6 不同正則化約束下二維彈性波全波形反演所得到的P 波速度結(jié)果.（a-c）分別表示基于Tikhonov、全變差和改進的全變差正則化的結(jié)果（修改自 Lin and Huang,2014 中的圖10）Fig.6 P wave velocity models from an elastic full waveform inversion with different regularization schemes.Panels(a-c) are results constrained with Tikhonov,Total variation and modified Total variation regularization,respectively (modified from Figure 10 in Lin and Huang,2014)

2.5 分辨率和不確定性

相對于射線層析成像，全波形反演能夠幫助我們通過擬合實際觀察和模擬所得到的波形來更好地提取地球內(nèi)部物質(zhì)信息.在這一過程當(dāng)中，我們能夠考慮復(fù)雜的物理過程，比如波前彌合，有限頻和多路徑傳播等.但是對于伴隨層析成像和全波形反演，目前的一大挑戰(zhàn)是如何評估所得到的模型的分辨率以及不確定性.由于在全波形反演當(dāng)中，我們需要很高的計算量以及較大的模型空間，目前我們尚無法得到傳統(tǒng)的分辨率矩陣的信息（Backus and Gilbert,1968,1970）.同時因為龐大的計算量，我們也很難像射線層析成像一樣進行棋盤（checkerboard）實驗（Leveque et al.,1993）.在過去的一段時間里，我們有以下幾種方法來分析相關(guān)的成像分辨率和不確定性：（1）Fichtner 和Trampert（2011）以及Zhu 等（2016）使用點擴散函數(shù)來分析某個區(qū)域的分辨率，以及多參數(shù)反演下模型參數(shù)之間的互相干擾（trade off）.這些點擴散函數(shù)是Hessian 矩陣與模型在某個區(qū)域的擾動量相乘的結(jié)果.比如，圖7 顯示在一個模擬反演當(dāng)中，不同的地區(qū)由于射線覆蓋、背景速度等的影響，所得到的點擴散函數(shù)具有比較大的差別.這一分析方法的缺點是計算量比較大，尤其是當(dāng)我們要分析多個區(qū)域的分辨率時；（2）Bui-Thanh 等（2013）在貝葉斯反演框架下重新構(gòu)建全波形反演.通過使用low-rank 近似，他們能夠得到后驗協(xié)方差矩陣的特征值和特征向量.這樣能幫助我們不僅得到某個局部最優(yōu)解，而且可以得到在這一特定結(jié)果下的后驗概率分布信息（Zhu et al.,2016）；（3）Fichtner 和van Leeuwen（2015）使用一些隨機樣本的互相關(guān)來得到方向以及位置相關(guān)的分辨率信息，以及各個參數(shù)之間的互相干擾.他們的結(jié)果顯示使用多個隨機樣本就可以提取出相關(guān)的分辨率信息，這一方法被應(yīng)用到歐洲上地幔成像當(dāng)中；（4）Liu 和Peter（2019）使用square-root variable metric 的方法來取代傳統(tǒng)的共軛梯度算法，從而在不需要矩陣參與計算的情況下來構(gòu)建Hessian矩陣，進而分析后驗概率分布的信息；（5）Thurin等（2019）使用卡曼濾波（Karman filtering）來分析全波形反演的不確定性.這些方法都在一定程度上得到全波形反演的分辨率以及不確定信息.

圖7 使用點擴散函數(shù)分析二維全波形反演所得到的不同位置上的分辨率.（a-i）分別顯示點擴散函數(shù)在右下角圖中不同區(qū)域的結(jié)果（修改自 Fichtner and van Leeuwen,2015 中的圖5）Fig.7 Using point-spread functions to analyze resolution at different locations in an 2D full waveform inversion.Panels (a-i) illustrate results at different locations as shown in the bottom right panel (modified from Figure 5 in Fichtner and van Leeuwen,2015)

2.6 多參數(shù)反演

如2.2 和2.3 節(jié)所示，能夠在全波形反演這一統(tǒng)一的框架下構(gòu)建不同地震學(xué)模型參數(shù)的梯度，進而反演它們的三維分布.但是目前大多數(shù)的方法只是基于一階導(dǎo)數(shù)的信息，所以就難以分析并且減少各個不同參數(shù)之間的互相串?dāng)_（孫史磊等，2020）.而要提取這一信息，需要考慮Hessian 矩陣的影響（Pratt et al.,1998）.一種可以幫助我們更好地設(shè)計不同參數(shù)組合的方法是通過理論分析各個參數(shù)組合對于入射地震波的脈沖響應(yīng).希望找到一個參數(shù)組合，使得各個參數(shù)的脈沖響應(yīng)在方位上的重合盡可能地?。╓u and Aki,1985）.如圖8 所示，相對于P 波速度和密度的組合，P 波速度和波阻抗的組合能夠更好地區(qū)分寬角度和窄角度的輻射響應(yīng)（Operto et al.,2013）.所以后一組合在多參數(shù)反演中具有更好的效果.另外，在多參數(shù)反演下，Hessian矩陣的大小由參數(shù)的個數(shù)以及模型網(wǎng)格點數(shù)目所決定，它是一個區(qū)塊化的對稱矩陣.區(qū)塊之間的數(shù)值表示了各個參數(shù)對反演的影響，以及它們之間互相干擾的情況.如果能夠構(gòu)建Hessian 矩陣以及計算它的逆矩陣，進而在迭代中考慮這一因素，則可以更好地進行多參數(shù)聯(lián)合反演（劉璐等,2013;Operto et al.,2013）.但是如前2.5 節(jié)所述，目前我們?nèi)匀缓茈y得到Hessian 矩陣以及它的逆矩陣，所以大部分的研究在于近似Hessian 矩陣，比如：（1）使用理論方法來近似Hessian 矩陣的對角分量，并在共軛梯度迭代當(dāng)中將他們當(dāng)作預(yù)調(diào)節(jié)算子（preconditioner）.這樣能夠在多參數(shù)聯(lián)合反演當(dāng)中減少各個參數(shù)之間的互相干擾，同時能夠更好地平衡深部與淺部的反演結(jié)果；（2）在迭代過程中使用不同的算法來近似地構(gòu)建Hessian 矩陣.這些算法包括L-BFGS、truncated Gauss Newton（Métivier et al.,2014）以及前面所提到的square root metric 方法（Liu and Peter,2019）.它們都可以相對有效地減少各個參數(shù)之間的互相干擾.

圖8 散射波對于不同模型參數(shù)組合的輻射樣式.（a,b）分別表示散射波對于聲波速度和密度組合所得到的輻射樣式.（c,d）分別表示聲波速度和波阻抗組合的結(jié)果.藍色和白色五角星分別表示激發(fā)震源和散射點.綠色曲線表示由射線加上波恩近似所得到的理論波前（修改自 Operto et al.,2013 中的圖2）Fig.8 Radiation patterns of scattered waves for different combinations of model parameters.Panels (a,b) show the radiation patterns for P wave velocity and density,respectively.Panels (c,d) are results for the combination of P wave velocity and impedance.Blue and white stars denote exciting source and scattering point.Green curves represent the amplitudes of scattered waves from Ray+Born approximation (modified from Figure 2 in Operto et al.,2013)

2.7 與背景噪聲記錄的結(jié)合

從2000 年開始，背景噪聲成像已經(jīng)被廣泛地應(yīng)用于研究地殼波速（Lin et al.,2008;Shapiro et al.,2005;Yang et al.,2007）、各向異性（Huang et al.,2010;Lin et al.,2011;Moschetti et al.,2010;Yao et al.,2010），以及觀測地殼速度隨時間變化（Brenguier et al.,2008a,2008b）.在傳統(tǒng)的背景噪聲成像中，首先計算背景噪聲互相關(guān)函數(shù)（Bensen et al.,2007;Yao et al.,2006），然后提取面波的相速度或者群速度的頻散曲線（Barmin et al.,2001）.通過使用一維速度反演（比如馬爾科夫鏈蒙特卡羅MCMC）來擬合這些頻散曲線（Shen et al.,2013）.最后通過插值這些一維速度模型來得到地球的三維速度結(jié)構(gòu).

背景噪聲成像數(shù)據(jù)也被應(yīng)用于全波形反演當(dāng)中.最常見的方法是假設(shè)通過背景噪聲互相關(guān)所得到的記錄能夠很好地近似兩個臺站之間的格林函數(shù).通過擬合觀測的經(jīng)驗格林函數(shù)和理論地震圖，可以用全波形反演來約束三維速度模型.這一方法已經(jīng)被應(yīng)用于研究西藏東南部的地殼速度（Chen et al.,2014）以及卡斯凱迪亞地殼與上地幔結(jié)構(gòu)當(dāng)中（Gao and Shen,2014）.另外，Lee 等（2014）和Wang 等（2020）結(jié)合天然地震記錄以及背景噪聲互相關(guān)結(jié)果，分別反演美國南加州地殼速度和各向異性分布.如圖9 所示，相對于只使用天然地震記錄所得到的模型，加入背景噪聲互相關(guān)能夠更好地約束模型區(qū)域的速度擾動.另外，Zhu（2018b）使用Rayleigh 面波互相關(guān)來研究美國得克薩斯州北部和俄克拉何馬州的地殼速度.相對于天然地震記錄，使用背景噪聲成像能夠研究構(gòu)造活動比較弱、地震震源分布不均勻區(qū)域的地殼波速結(jié)構(gòu).美國得克薩斯州北部和俄克拉何馬州在2008 年以前沒有太多的地震，所以相關(guān)的研究比較少.新的地殼模型能夠更好地研究這一區(qū)域的誘發(fā)地震的震源特征（Ellsworth,2013）.以上的方法都是基于垂直分量Rayleigh 面波的經(jīng)驗格林函數(shù)，最近Wang 等（2019）將這一方法推廣到三分量格林函數(shù)噪聲伴隨成像當(dāng)中，并在南加州各向異性成像中取得了很好的應(yīng)用效果（Wang et al.,2020）.

圖9 通過使用背景噪聲信號改進南加州地殼剪切波速度模型.從左到右分別是M16（只使用天然地震記錄）、M21（結(jié)合天然地震記錄和背景噪聲互相關(guān)函數(shù)），以及兩者之間的差（Diff）（修改自 Wang et al.,2019 中的圖9）Fig.9 Comparison of shear wave velocity models by incorporating ambient noise cross correlation functions.From left to right are model M16 (only constrained by earthquake records),M21(further constrained by ambient noise cross correlation functions)and differences between M16 and M21 (modified from Figure 9 in Wang et al.,2019)

上述研究中的一個重要假設(shè)是使用背景噪聲互相關(guān)所得到的記錄能夠很好地近似兩個臺站之間的經(jīng)驗格林函數(shù).這一假設(shè)只有在以下條件滿足的情況下才成立：（1）激發(fā)背景噪聲的震源必須是均勻并且同勢能分布.（2）所使用的噪聲記錄時間足夠長，進而得到的互相關(guān)函數(shù)能夠穩(wěn)定地收斂到經(jīng)驗格林函數(shù).然而這些假設(shè)在實際情況下是很難得到滿足的.為了解決這一問題，Tromp 等（2010）提出了在伴隨層析成像的框架下，通過對比觀測和計算的集合互相關(guān)函數(shù)，來同時反演地球速度結(jié)構(gòu)以及背景噪聲的震源分布.這樣能夠避免使用均勻且同勢能噪聲來源的假設(shè).這一方法目前被Ermert 等（2016）用來研究激發(fā)全球背景噪聲的震源分布.

2.8 地殼速度反演

相對于傳統(tǒng)的射線走時層析成像，全波形反演具有更高的分辨率以及能夠更好地反演具體的速度擾動大小.這里列舉一些重要的地殼成像方面的應(yīng)用.（1）盆 地 結(jié) 構(gòu).Tape 等（2009）和Lee 等（2014）早期將伴隨層析成像和全波形反演應(yīng)用到研究南加州地殼當(dāng)中.他們發(fā)現(xiàn)相對于以前的走時層析成像所得到的三維地殼模型，通過擬合三分量的地震圖，新的模型中所成像的盆地結(jié)構(gòu)更加清晰，尤其是盆地以及近地表的低速帶相對于初始模型的擾動可以到達20%～30%.這么大的速度擾動是很難通過傳統(tǒng)的射線走時層析成像得到的（圖10）；（2）斷層結(jié)構(gòu).南加州具有很多大型的斷層，比如圣安德烈斯和Garlock 斷層等.在Tape 和Lee 等的三維模型中，我們可以清晰地看到跨過斷層帶時所出現(xiàn)的大的速度變化，這與早期的到時層析成像的結(jié)果相吻合（Lin et al.,2010），但是他們所得到的速度變化更大.相對于Tape 等的模型，Lee 等（2014）結(jié)合了背景噪聲所得到的格林函數(shù)和天然地震數(shù)據(jù)，所以相對來說具有更好的成像覆蓋率.同時上述模型與二維主動源反演速度模型擬合較好（圖10）.以上的研究都是使用各向同性的速度反演.最近，Wang 等（2020）在Tape模型的基礎(chǔ)上加入了背景噪聲數(shù)據(jù)，反演了各向同性及徑向各向異性的南加州地殼結(jié)構(gòu)；（3）下地殼結(jié)構(gòu).通過結(jié)合伴隨層析成像以及背景噪聲所得到的經(jīng)驗格林函數(shù)，Chen 等（2014）反演了西藏東南部的地殼結(jié)構(gòu).他們發(fā)現(xiàn)中下地殼（約20 km 厚）具有很低的剪切波速度（-6%～-12%），這可能來自于印度板塊和歐亞板塊的強烈碰撞所導(dǎo)致的下地殼部分熔融.另外，使用全波形反演所得到的三維地殼模型可以更好地研究地震震源，比如地震定位以及震源機制解（Lee et al.,2011;Liu et al.,2004;Wang and Zhan,2020）.

圖10 對比不同南加州地殼模型在二維LARSE-I（a）和LARSE-II（b）剖面上的P 波速度結(jié)果.CVM-S4.26、CVM-S4 和CVM-H11.9 分別是三個南加州地震研究中心的標(biāo)準(zhǔn)公共地殼模型.Lutter 等（1999,2004）以及Fuis 等（2003）分別是兩個二維主動源地殼速度反演結(jié)果（修改自 Lee et al.,2014 中的圖8）Fig.10 Comparisons of P wave velocities in different Southern Californian velocity models for 2D LARSE-I and LARSE-II cross sections.Models CVM-S4.26,CVM-S4 and CVM-H11.9 are three 3D reference crustal velocity models from the Southern California Earthquake Center (SCEC).Lutter et al.(1999,2004) and Fuis et al.(2003) are 2D crustal velocity models constrained by active source experiments (modified from Figure 8 in Lee et al.,2014)

上述地殼反演的另一個重要要求是，需要設(shè)計一個模擬區(qū)域能夠把盡可能多的震源和臺站包括其中.對于地震活動不頻繁、但有密集地震臺陣的區(qū)域，常常無法只使用本地的震源來進行全波形反演.一種可以在地殼和上地幔反演中同時考慮遠震震源的方法是在三維反演區(qū)域之外，使用一些快速算法和較簡單的背景模型來計算遠震波場（比如FK 或者DSM 方法）.然后將計算所得到的波場當(dāng)作邊界條件，而在所要研究的區(qū)域三維模型中仍然使用數(shù)值算法（比如有限差分或者譜元法）來計算三維波場以及目標(biāo)函數(shù)梯度.這種混合數(shù)值模擬波場的計算方法使得我們可以在三維地殼波速反演以及反射面成像當(dāng)中使用遠震記錄（Beller et al.,2017;Monteiller et al.,2015,2021;Pienkowska et al.,2021;Tong et al.,2014;Wang et al.,2016），從而大大地改進密集臺陣下面的結(jié)構(gòu)分辨率以及成像覆蓋率.例如，Wang 等（2016）展示只要用極少數(shù)的幾個遠震，并使用一個近線性臺陣上的P 波及其散射和折射波的記錄，就可以反演出歐洲西南部比利牛斯山脈（Pyrenees）下方巖石圈的縱波和橫波速度.相似的方法也被應(yīng)用到反演區(qū)域在震源和臺站中間的情況.從震源用快速算法和簡單模型傳到模擬區(qū)域的地震波經(jīng)過三維數(shù)值模擬后，仍需要通過邊界條件用快速算法和簡單地球模型傳回到接收臺站.這類混合數(shù)值計算波場的方法被用于全波形反演中，也被稱為box tomography，并被應(yīng)用于上地幔成像當(dāng)中（Clouzet et al.,2018;Masson and Romanowicz,2017）.

2.9 地幔速度反演

如上所述，伴隨層析成像和全波形反演已經(jīng)被廣泛地應(yīng)用于很多地幔尺度的研究當(dāng)中.本文列舉一些與板塊俯沖帶和地幔柱相關(guān)的典型研究成果.

（1）板塊俯沖帶成像.研究海洋板塊俯沖的幾何特征以及其波速結(jié)構(gòu)對于我們了解板塊構(gòu)造、地球動力學(xué)以及礦物物理具有十分重要的意義（Stern,2002）.目前在大陸尺度上，全波形反演能夠更好地刻畫海洋巖石圈在俯沖過程當(dāng)中所經(jīng)歷的復(fù)雜過程（Tao et al.,2018;Zhu et al.,2012）.相關(guān)的成像結(jié)果與傳統(tǒng)射線走時層析成像結(jié)果相當(dāng)接近.但同時通過提高模型分辨率以及更好地反演具體的波速大小，我們也發(fā)現(xiàn)了一些新的現(xiàn)象.比如Zhu 等（2012）在歐洲區(qū)域的研究中能夠清晰地看到板片滑脫（slab detachment），這一現(xiàn)象與以前走時層析成像的結(jié)果相似（Wortel and Spakman,2000），但在新的模型中我們能夠找到更多的區(qū)域具有此滑脫現(xiàn)象.另外，現(xiàn)實中海洋板塊在俯沖過程當(dāng)中會發(fā)生斷裂，而且這一斷裂面會沿著海溝方向不斷擴大，進而導(dǎo)致滑脫不斷增大.這些斷裂會導(dǎo)致軟流圈物質(zhì)上涌，進而與地表的火山作用相關(guān).最近，Zhu 等（2020c）在中南美洲的全波形反演結(jié)果中發(fā)現(xiàn)Rivera 和Cocos 板塊在俯沖過程中也出現(xiàn)類似的滑脫現(xiàn)象.同時在板塊斷裂區(qū)間填充有大量的低速帶物質(zhì)，其中可能發(fā)生部分熔融現(xiàn)象.另一方面，通過反演方位各向異性，他們發(fā)現(xiàn)板塊的快速后退以及板塊之間的空隙會產(chǎn)生巨大的壓強差從而導(dǎo)致局部的軟流圈物質(zhì)流動.這一流動會使得橄欖石等各向異性礦物的快速軸方向沿著流動/應(yīng)變方向定向排布，從而產(chǎn)生可觀測的地震波各向異性.這對于我們研究俯沖板塊與周圍的地幔物質(zhì)相互對流作用具有十分重要的意義.另一方面，在西太平洋和東北亞，早期的射線走時層析成像已經(jīng)發(fā)現(xiàn)太平洋板塊在俯沖過程中與上地幔，特別是地幔轉(zhuǎn)換帶具有復(fù)雜的關(guān)系（Fukao and Obayashi,2013;Zhao et al.,1992,1994）.比如在東北亞，太平洋板塊在地幔轉(zhuǎn)換帶中可以橫向地延展上千千米.同時Tang 等（2014）的研究顯示在地幔過渡帶中，太平洋板塊發(fā)生斷裂，導(dǎo)致地幔物質(zhì)上涌，進而與地表的長白山五大連池火山相連.Tao 等（2018）通過使用全波形反演更好地描述這一俯沖板塊在上地幔，特別是地幔過渡帶中的特征.如圖11 所示，通過擬合地幔過渡帶的SH 波形，相對于以前的成像結(jié)果，他們的俯沖板塊分辨率更高，成像結(jié)果更加清晰.在全球尺度上，相關(guān)的伴隨層析成像結(jié)果也更好地顯示俯沖板塊在下地幔的具體形態(tài).而這些特征在傳統(tǒng)的走時層析成像結(jié)果當(dāng)中常常并不十分清楚，進而導(dǎo)致模型解釋上的一系列爭論.通過更加清晰地描繪俯沖板塊在地幔過渡帶和下地幔的幾何學(xué)特征，能夠更好地幫助我們了解地幔對流以及地球內(nèi)部的動力學(xué)過程（Lei et al.,2020）.

圖11 使用通過地幔轉(zhuǎn)換帶的SH 波改進俯沖板塊伴隨層析成像結(jié)果.（a）顯示的是二維剖面和所使用的地震臺站的位置.（b）對比初始模型（左）和迭代改進之后的模型（右）對SH 地幔轉(zhuǎn)換波的影響.黑色和紅色地震圖分別是實際觀測和模擬計算所得到的結(jié)果.（c）和（d）分別是初始模型和迭代改進之后的模型在（a）中所示的二維剖面上的剪切波速度擾動（修改自 Tao et al.,2018 中的圖6）Fig.11 Imaging subducting slabs by using SH triplication waveforms.Panel (a) shows the locations of 2D vertical cross sections and stations.Panel (b) compares observed (black) and synthetic (red) SH waveforms from the starting (left) and updated models(right).Panels (c) and (d) show shear wave velocity perturbations for the starting and updated models on the 2D vertical cross section shown in Panel (a) (modified from Figure 6 Tao et al.,2018)

（2）地幔柱成像.相對于海洋板塊的俯沖，地幔柱提供了地幔對流的上升區(qū)域.從地幔柱假說提出至今（Morgan,1971），地震學(xué)界對地幔柱的成像一直十分關(guān)注，而且不同的成像方法常常得到了不同結(jié)論（Montelli et al.,2004b,2006）.比如有的研究顯示冰島的地幔柱來自于核幔邊界（Bijwaard and Spakman,1999），而有的研究顯示這一區(qū)域的低速帶主要局限于上地幔（Ritsema et al.,1999）.相對于高波速的海洋板塊俯沖帶來說，低波速的地幔柱對于成像精度和分辨率的要求更大.由于波前愈合以及成像分辨率的局限，傳統(tǒng)的走時層析成像很難反演幾百千米直徑的上升低速地幔柱.相反，全波形反演考慮了上述的物理過程，從而能夠更好地幫助我們反演尺度更小的低波速帶.比如Richers 等（2013）發(fā)現(xiàn)冰島地幔柱在上地幔中向東部偏移，可能與地幔對流方向有關(guān) .在全球尺度上，新的成像結(jié)果發(fā)現(xiàn)了更多的地幔柱，它們連接核幔邊界和地表的熱點.French 和Romanowicz（2015）還發(fā)現(xiàn)這些地幔柱與軟流圈上垂直于洋中脊的低速帶相連.值得注意的是這一研究使用的是譜元方法來正演模擬地震波傳播，但是使用的是近似的方法（nonlinear asymptotic coupling theory,NACT）來計算敏感核.相對于伴隨層析成像和全波形反演，這樣能夠大大地減少計算量.如圖12 所示，Lei 等（2020）在全球伴隨層析成像模型當(dāng)中發(fā)現(xiàn)南非地幔柱在上升到地幔轉(zhuǎn)換帶時，分成多個更小的地幔柱，然后上升連接到地表的火山和熱點.這些結(jié)果大大地提高了對地幔柱的認識，從而為更好地了解地幔對流提供了巨大的幫助.

圖12 對比三個全球尺度地幔剪切波成像模型在6 個二維剖面上的結(jié)果.從左到右分別是來自于GLAD-M25（Lei et al.,2020）、TX2015（Lu and Grand,2016）和SEMUCB-WM1（French and Romanowicz,2015）.（a-f）分別對應(yīng)于以下的熱點:（a）阿法爾州;（b）百慕大群島和加那利群島;（c）佛得角和哈加爾高原;（d）冰島和艾費爾高原;（e）復(fù)活節(jié)島 和加拉帕戈斯群島；（f）馬里昂縣和凱爾蓋朗群島（修改自 Lei et al.,2020 中的圖15）Fig.12 Comparisons of shear wave velocity perturbations in three global-scale mantle tomography models.From left to right are results from GLAD-M25 (Lei et al.,2020)，TX2015 (Lu and Grand,2016) and SEMUCB-WM1 (French and Romanowicz,2015).Panels (a) to (f) show results for the following hotspots: (a) Afar;(b) Bermuda and Canary;(c) Cape Verde and Hoggar;(d) Iceland and Eifel;(e) Easter and Galapagos and (f) Marion and Kergulen (modified from Figure 15 in Lei et al.,2020)

2.10 油氣儲層反演

雖然我們對全波形反演和伴隨方法早在1980 年代就已經(jīng)在理論和實踐上進行了探討（Tarantola,1984,1988），但是直到2007 年左右，因為其在二維盲測實驗當(dāng)中所取得的成功（Brenders and Pratt,2007），才又一次得到人們的關(guān)注.這來源于過去幾十年中計算能力和算法精度的提高.目前全波形反演已經(jīng)被廣泛地應(yīng)用于能源勘探項目當(dāng)中，為復(fù)雜條件下尋找油氣資源提供了巨大的幫助（Operto et al.,2015;Operto and Miniuss,2018）.比如，圖13 對比在北海Vahall 油田，通過擬合海底地震臺站觀測波形所得到的速度模型與傳統(tǒng)的反射到時層析成像的結(jié)果，其中全波形反演所得到的結(jié)果能夠更好地尋找油氣儲層.相對于經(jīng)典的全波形反演，目前在勘探領(lǐng)域有以下的一些發(fā)展.

圖13 對比全波形反演在北海Valhall 油田勘探當(dāng)中的應(yīng)用.左圖和右圖分別是使用反射波到時層析成像和全波形反演所得到的結(jié)果.（a-c）以及（h-j）分別是在175 m、500 m 和1 000 m 所得到的P 波速度結(jié)果.（d-g）以及（k-n）分別是相應(yīng)虛線位置的二維剖面結(jié)果（修改自 Operto et al.,2015 中的圖4 和圖11）Fig.13 Comparison of velocity models for the ocean bottom node data from the Valhall oilfield.The left panel and right panel are models from reflection travel time tomography and full waveform inversion,respectively.(a-c) and (h-j) in each panel are P wave velocities at 175 m,500 m and 1000 m.(d-g) and (k-n) demonstrate results in corresponding vertical profiles as shown by dashed lines (modified from Figures 4 and 11 in Operto et al.,2015)

（1）時間與頻率域的選擇.相對于早期的時間域正演和反演，使用頻率域、Laplace 域以及Fourier-Laplace 域的反演具有以下幾個優(yōu)點（卞愛飛等,2010;Pratt and Worthington,1990;Pratt,1999;Pratt and Shipp,1999;Shin and Cha,2008,2009）：第一，如果能夠在三維空間下快速而準(zhǔn)確地求解Helmholtz 方程，那么在頻率域以及Laplace 域的計算量將與震源數(shù)目無關(guān).這大大地減少了全波形反演的計算量，對大型三維勘探項目尤其重要（Engquist and Ying,2011;Operto et al.,2007;Poulson et al.,2013）.第二，在頻率域以及Laplace 域中，可以更好地使用多尺度反演策略.相關(guān)的研究表明，只需要使用幾個特定的頻率或者頻率組合就能夠很好地覆蓋所要研究對象的波數(shù)分布（Bunks et al.,1995）.第三，相對于時間域，在頻率域下可以更直接地描述地震波的衰減，進而更好地反演品質(zhì)因子Q的結(jié)構(gòu)（Kamei and Pratt,2013;Prieux et al.,2013）.

（2）反射波全波形反演.早期的二維全波形理論實驗表明，如果只使用透射波，比如井間層析成像，則可以很好地約束低波數(shù)信息，即光滑的速度結(jié)構(gòu).相反，如果只使用反射波（地表接收信號），則主要只能約束高波數(shù)的速度分布（Gauthier et al.,1986;Mora,1987）.如何使用反射波來約束低波數(shù)的速度模型一直以來是一個重要的研究問題.通過構(gòu)建反射波的敏感核，可以使用反射波的到時以及波形來達到約束低波數(shù)的速度模型的效果.相對于直達P 波的敏感核，反射波的敏感核被形象地稱為Rabbit ear（Ma and Hale,2013;姚剛和吳迪,2017;Zhou et al.,2015）.它與傳統(tǒng)的偏移脈沖響應(yīng)（migration impulse response）不同，后者主要描述高波數(shù)的速度結(jié)構(gòu)（反射率）（圖14）.

圖14 在二維均勻速度模型下對比直達波和反射波的梯度.（a）和（b）分別是直達波和反射波的梯度，其中反射界面如圖（b）中的紫色直線所示.（c）直達波+反射波的梯度.（d）改進的聯(lián)合全波形梯度結(jié)果（修改自 Zhou et al.,2015 中的圖3）Fig.14 Gradients for direct and reflected waves in a 2D homogeneous velocity model.Panels (a) and (b) are gradients for direct and reflected P waves.The reflector is shown as the purple line in Panel (b).Panels (c) and (d) are gradients for direct+reflected waves,and joint full waveform inversion gradient (modified from Figure 3 in Zhou et al.,2015)

（3）高頻成像與反演.在傳統(tǒng)勘探開發(fā)當(dāng)中，我們使用層析成像來反演光滑的低波數(shù)速度結(jié)構(gòu).然后使用偏移方法，比如逆時偏移成像（McMechan,1983），來找到相應(yīng)的高波數(shù)反射界面以及油氣儲層位置.隨著全波形反演的發(fā)展，它被逐漸地推廣到高頻段.目前的反演頻率可以達到30～50 Hz.如果仍然能夠使這一非線性反演問題得到收斂，而且避免前述的周期跳躍問題，在未來有希望使用全波形反演來替代以前的層析成像+偏移的勘探開發(fā)步驟.

（4）與全局反演的結(jié)合.如前所述，全波形反演本質(zhì)上是一種求解局部非線性最優(yōu)化的問題.它只能找到目標(biāo)函數(shù)的局部最優(yōu)解，而無法保證得到全局最優(yōu)，同時也很難提供相關(guān)的不確定性信息.目前，有相關(guān)的研究結(jié)合數(shù)值正演以及全局反演方法，比如MCMC 或者Hamiltonian MCMC 方法，來進行速度模型的采樣（Gebraad et al.,2020;Stuart et al.,2019;Zhao and Sen,2021）.后者相對于前者，由于使用了目標(biāo)函數(shù)梯度的信息，能夠使得反演過程更快地收斂（圖15）.這一點對該類方法十分重要.由于這些取樣涉及到上千甚至上萬的正演模擬，在目前的計算能力下做三維規(guī)?；茝V還是十分具有挑戰(zhàn)性的.所以目前大部分類似的研究還是在二維模擬反演當(dāng)中.如果能夠使用基函數(shù)來更好地參數(shù)化模型，則有希望大幅地減少模型空間的大小，進而加快采樣的收斂速率.

圖15 使用基于目標(biāo)函數(shù)梯度的馬爾科夫鏈蒙特卡羅方法反演二維Marmousi 模型.（a-c）分別表示平均模型、標(biāo)準(zhǔn)差模型以及平均與真實模型之間的差.（d）顯示6 個采樣點的一維邊界后驗概率分布情況（修改自 Zhao and Sen,2021 中的圖9 和10）Fig.15 Results for the 2D Marmousi model from a gradient based Markov chain Monte Carlo sampling.Panels (a) to (c) are results for the mean velocity model,standard deviation and differences between the mean and true velocity models.Panel (d) shows the 1D marginal posteriori probability density functions at six different locations (shown in Panel a) (modified from Figures 9 and 10 in Zhao and Sen,2021)

其他的勘探領(lǐng)域的全波形反演問題涉及到前面所述的設(shè)計優(yōu)化的目標(biāo)函數(shù)，多參數(shù)聯(lián)合反演，與機器學(xué)習(xí)的結(jié)合以及在實際數(shù)據(jù)中提取低頻信號等，這里不再做贅述.

3 結(jié)論與展望

全波形反演目前已經(jīng)被廣泛地應(yīng)用于多尺度的地球內(nèi)部成像研究當(dāng)中，它可以更好地尋找油氣資源，描述巖石圈的波速結(jié)構(gòu)進而更好地研究地震震源分布和機制，同時還可以用于研究板塊俯沖過程以及地幔柱的結(jié)構(gòu).這些成果已經(jīng)幫助我們更好地研究地球內(nèi)部的物質(zhì)組成以及動力學(xué)過程.相關(guān)的研究成果對巖石礦物物理和地球動力學(xué)具有十分重要的科學(xué)意義.目前主要的研究問題還包括選取優(yōu)化的目標(biāo)函數(shù)，如何評價反演結(jié)果的分辨率以及不確定性，多參數(shù)聯(lián)合反演以及反演結(jié)果的解釋等.相關(guān)方面的突破在未來可以更好地了解地球以及其他行星的內(nèi)部物理及動力學(xué)過程.

致謝

衷心感謝中國科學(xué)技術(shù)大學(xué)姚華建教授的邀請，以及三位評審專家提出的寶貴修改意見.同時感謝編輯部老師在稿件修改和發(fā)表過程中的大力支持.