祝賀君，劉沁雅，楊繼東

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

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

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

0 引言

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

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

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

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

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

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

自從全波形反演理論的提出，這種新方法已經(jīng)被應(yīng)用到各種尺度的地球內(nèi)部結(jié)構(gòu)成像研究當(dāng)中.比如，Tape 等和Chen 等使用伴隨/全波形成像來(lái)研究南加州地殼的三維速度結(jié)構(gòu)（Chen et al.,2007b;Lee et al.,2014;Tape et al.,2009,2010）.他們的模型能夠幫助我們很好地研究這一區(qū)域的地質(zhì)構(gòu)造，并且已被美國(guó)南加州地震中心（Southern California Earthquake Center）作為區(qū)域標(biāo)準(zhǔn)的三維公共速度模型（比如最新的CVM-H15.1.0 和CVMS4.26）.同時(shí)這些高分辨率高精度的三維速度模型能夠幫助我們更好地研究南加州地震的震源信息（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）.同時(shí)，全波形反演也在油氣資源開發(fā)過(guò)程中發(fā)揮越來(lái)越重要的作用（Operto et al.,2015;楊積忠等,2014）.我們將會(huì)在2.8—2.10 節(jié)中詳細(xì)討論這些具體的例子.迄今為止，已經(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 理論框架

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

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

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

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

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

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

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

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

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

1.2 模型參數(shù)化

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

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

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

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

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

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

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

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

圖1 二維均勻速度模型下構(gòu)建SH 波的敏感核.從左到右分別是正傳波場(chǎng)、伴隨波場(chǎng)、相互作用波場(chǎng)和剪切波速度的敏感核.五角星和方塊分別表示震源和接收臺(tái)站的位置.時(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+分別是正傳和伴隨波場(chǎng)，D和D+則是正傳和伴隨應(yīng)變偏量（strain deviator），T是記錄時(shí)間.

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

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

1.4 最優(yōu)化迭代

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

圖2 對(duì)比九個(gè)不同理論模型下使用三種迭代方法所得到的收斂性，包括最速下降法（黑線），非線性共軛梯度法（綠線）和L-BFGS 法（紅線）.不同的理論模型如每個(gè)圖的標(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）使用近似方法來(lái)構(gòu)建Hessian 矩陣.因?yàn)镠essian 矩陣本身十分巨大，它的大小與模型的網(wǎng)格點(diǎn)數(shù)的平方成正比（Pratt et al.,1998）.所以目前仍然很難直接地計(jì)算這一矩陣.但是通過(guò)各種近似，可以構(gòu)建它的對(duì)角量，使用該信息可以幫助我們更好地平衡模型淺部與深部的梯度大小，并且可以減少在震源以及接收臺(tái)站上模型梯度數(shù)值較大的問(wèn)題.（2）使用高斯-牛頓方法來(lái)迭代地考慮Hessian 矩陣的作用（Epanomeritakis et al.,2008;劉璐等,2013;Métivier et al.,2014），這被稱之為截?cái)喔咚古ｎD（Truncated Gauss-Newton）方法.在這一方法中，在外部共軛梯度或者L-BFGS 循環(huán)當(dāng)中加入一個(gè)內(nèi)循環(huán)，通過(guò)求解法方程式（Normal equation）來(lái)迭代地考慮Hessian 矩陣的作用.這一方法已被用于勘探地震學(xué)的數(shù)值試驗(yàn)當(dāng)中，但是相對(duì)于一般的共軛梯度或者L-BFGS 方法，它涉及到內(nèi)外兩個(gè)迭代循環(huán)，所以計(jì)算量仍然十分巨大.

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

2 研究進(jìn)展

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

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

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

圖3 對(duì)比不同目標(biāo)函數(shù)的特征.（a）顯示的是子波波形.（b）和（c）分別顯示的是基于最小二乘波形殘差和最優(yōu)化路徑所得到的目標(biāo)函數(shù)隨著不同時(shí)間移動(dòng)（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 各向異性反演

我們使用各向異性來(lái)描述地震波沿不同方向所具有的不同傳播速度，在地震學(xué)研究中包括徑向各向異性和方位各向異性（Backus,1965;Hess,1964;Montagner and Nataf,1986;Smith and Dahlen,1973）.這些各向異性來(lái)自于形狀或是晶格的定向排列[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）來(lái)研究地球內(nèi)部的各向異性分布.相對(duì)來(lái)說(shuō)，如果有很好的臺(tái)站分布，橫波分裂可以提供很好的橫向分辨率.但是由于使用的是垂向入射的遠(yuǎn)震體波信號(hào)，所以對(duì)于垂向上的各向異性分布具有很差的分辨率.相反，因?yàn)椴煌l率的面波對(duì)不同深度的速度具有不同的敏感性,使用面波頻散可以幫助我們提高垂向上的分辨率.但是由于在面波相速度或群速度反演當(dāng)中，使用最小二乘法來(lái)求解一個(gè)層析成像問(wèn)題，通常需要加入正則化約束.同時(shí)由于面波的橫向傳播，相對(duì)來(lái)說(shuō)，它的橫向分辨率要低于垂向分辨率.通過(guò)研究地球內(nèi)部的各向異性，可以討論地殼、巖石圈與軟流圈中的應(yīng)力、形變以及流動(dòng)性特征，同時(shí)也可以用來(lái)研究地幔中的對(duì)流以及地核構(gòu)造（Long and Silver,2008;Long and Becker,2010;Park and Levin,2002）.

如前所述，伴隨層析成像和全波形反演提供了一個(gè)統(tǒng)一的框架來(lái)反演各種不同模型參數(shù)的組合.所以它也可以用于反演各向異性（Rusmanugroho et al.,2017;Sieminski et al.,2007a,2007b）.比如在地殼成像方面，Wang 等（2020）在Tape 等的南加州地殼模型基礎(chǔ)之上，加入了徑向各向異性.他們發(fā)現(xiàn)在圣安德烈斯斷層西面的中下地殼中，橫向偏振的剪切波速度比縱向偏振要低（負(fù)徑向各向異性～-6%）.而跨過(guò)斷層，這一比例正好相反.這與上地殼和上地幔中常見的徑向各向異性特征正好相反.為了解釋這一現(xiàn)象，他們提出一種可能性是由于下地殼的斜長(zhǎng)石具有與黑云母和橄欖石不一樣的各向異性特征（Wang et al.,2020），導(dǎo)致它們的快速軸與應(yīng)變方向垂直而非平行.在地幔成像方面，Zhu 等（2020c）使用伴隨方法反演上地幔的方位各向異性，進(jìn)而研究俯沖板塊與地幔對(duì)流的關(guān)系.例如，他們發(fā)現(xiàn)在中南美洲和卡斯凱迪亞俯沖帶中，當(dāng)出現(xiàn)板塊斷裂時(shí)，上地幔物質(zhì)會(huì)流動(dòng)通過(guò)這些板塊的間斷空間（圖4），從而在板塊窗以及板塊周圍產(chǎn)生環(huán)流（Hall et al.,2000;Russo and Silver,1994）.這些觀測(cè)與巖石物理實(shí)驗(yàn)（Buttles and Olson,1998;Kincaid and Griffiths,2003）和地球動(dòng)力學(xué)模擬相吻合（Schellart,2004;Schellart et al.,2007;Stegman et al.,2006）.另外，在歐洲上地幔反演結(jié)果當(dāng)中（Zhu and Tromp,2013），發(fā)現(xiàn)反演所得到的方位各向異性方向與過(guò)去三千萬(wàn)年的大地構(gòu)造過(guò)程有很好的相關(guān)性.同時(shí)通過(guò)分析方位各向異性的大小在垂向上的變化，他們推測(cè)出脆性和韌性變形在巖石圈中的分布，這一結(jié)果與實(shí)驗(yàn)和震源深度分布結(jié)果相吻合（Chen and Molnar,1983;Kohlstedt et al.,1995）.另外，通過(guò)同時(shí)反演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)致的地幔對(duì)流.（a）和（b）分別顯示中美洲和卡斯凱迪亞俯沖帶的反演結(jié)果.綠色塊體顯示的是剪切波速度擾動(dòng)大于1.5%的結(jié)果.黃色箭頭表示通過(guò)伴隨方位各向異性成像所得到的地幔對(duì)流情況[（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)部傳播過(guò)程中經(jīng)受衰減作用，其中動(dòng)能和彈性勢(shì)能被轉(zhuǎn)換成熱能.相對(duì)于過(guò)去30多年地震波波速層析成像的顯著進(jìn)步來(lái)說(shuō)（Grand et al.,1997;Meschede and Romanowicz,2015;Schaeffer and Lebedev,2013），地震波衰減成像的進(jìn)步相對(duì)較?。―alton et al.,2008;Gung and Romanowicz,2004;Romanowicz,1995）.相對(duì)于以前的結(jié)果，在最近的一些大陸和全球尺度的衰減成像當(dāng)中（圖5），能夠看到相對(duì)較好的可比性（Adenis et al.,2017a,2017b;Bao et al.,2016a,2016b;Ma et al.,2016;Karaoglu and Romanowicz,2017）.相對(duì)于波速成像，地震波衰減成像有以下的幾個(gè)難點(diǎn)：（1）需要有高質(zhì)量的地震波振幅信息，這就需要具有較高信噪比的地震記錄，同時(shí)要對(duì)地震儀器響應(yīng)進(jìn)行校正；（2）體波走時(shí)和面波頻散主要受地震波波速的影響，然而地震波振幅受到多重因素的影響，其中包括介質(zhì)衰減、波速、震源大小、輻射模式、彈性聚焦與散焦（Durek et al.,1993;Romanowicz,1987;Ruan and Zhou,2010,2012）、散射以及臺(tái)站局部的放大效應(yīng)等（Eddy and Ekstrom,2014;Lin et al.,2012）；（3）目前尚難以區(qū)分內(nèi)在衰減和由地震波散射所導(dǎo)致的表象衰減.所以在衰減成像當(dāng)中，需要仔細(xì)地考慮這些因素.

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

圖5 對(duì)比三個(gè)全球尺度的上地幔地震波衰減成像結(jié)果.（a-c）分別來(lái)自于QRLW8（Gung and Romanowicz,2004）、SEMUCB-UMQ（Karaoglu and Romanowicz,2018）和QRFSI12（Dalton et al.,2008）.每一個(gè)圖的下方第一行和第二行分別顯示的是各個(gè)深度的剪切模量衰減值及其擾動(dòng)量，其中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)中的一個(gè)，它具有以下幾個(gè)缺點(diǎn)：（1）品質(zhì)因子Q并不在所推導(dǎo)的地震波公式當(dāng)中，相反，需要求得相應(yīng)的應(yīng)力和應(yīng)變松弛時(shí)間因子（Carcione et al.,1988;Robertsson et al.,1994）.這給反演和正演都帶來(lái)了一定的困難；（2）如上所述，如果要同時(shí)獲得波速和衰減的梯度信息，需要計(jì)算兩次伴隨波場(chǎng)，從而加大了反演的計(jì)算量.為了解決上述問(wèn)題，在過(guò)去幾年當(dāng)中，一些用以描述地震波在地球內(nèi)部衰減現(xiàn)象的新的理論模型被提出.比如Fichtner 和Driel（2014）提出一個(gè)改進(jìn)的標(biāo)準(zhǔn)線性固體模型，使得Q直接出現(xiàn)在波場(chǎng)公式中，進(jìn)而簡(jiǎn)化了計(jì)算目標(biāo)函數(shù)相對(duì)于Q的梯度的過(guò)程.同時(shí)這種方法也可以幫助我們表示與頻率相關(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的全波形反演過(guò)程當(dāng)中（Yang et al.,2020）.

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

2.4 模型正則化

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

圖6 不同正則化約束下二維彈性波全波形反演所得到的P 波速度結(jié)果.（a-c）分別表示基于Tikhonov、全變差和改進(jìn)的全變差正則化的結(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 分辨率和不確定性

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

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

圖8 散射波對(duì)于不同模型參數(shù)組合的輻射樣式.（a,b）分別表示散射波對(duì)于聲波速度和密度組合所得到的輻射樣式.（c,d）分別表示聲波速度和波阻抗組合的結(jié)果.藍(lán)色和白色五角星分別表示激發(fā)震源和散射點(diǎn).綠色曲線表示由射線加上波恩近似所得到的理論波前（修改自 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），以及觀測(cè)地殼速度隨時(shí)間變化（Brenguier et al.,2008a,2008b）.在傳統(tǒng)的背景噪聲成像中，首先計(jì)算背景噪聲互相關(guān)函數(shù)（Bensen et al.,2007;Yao et al.,2006），然后提取面波的相速度或者群速度的頻散曲線（Barmin et al.,2001）.通過(guò)使用一維速度反演（比如馬爾科夫鏈蒙特卡羅MCMC）來(lái)擬合這些頻散曲線（Shen et al.,2013）.最后通過(guò)插值這些一維速度模型來(lái)得到地球的三維速度結(jié)構(gòu).

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

圖9 通過(guò)使用背景噪聲信號(hào)改進(jìn)南加州地殼剪切波速度模型.從左到右分別是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)

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

2.8 地殼速度反演

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

圖10 對(duì)比不同南加州地殼模型在二維LARSE-I（a）和LARSE-II（b）剖面上的P 波速度結(jié)果.CVM-S4.26、CVM-S4 和CVM-H11.9 分別是三個(gè)南加州地震研究中心的標(biāo)準(zhǔn)公共地殼模型.Lutter 等（1999,2004）以及Fuis 等（2003）分別是兩個(gè)二維主動(dòng)源地殼速度反演結(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)

上述地殼反演的另一個(gè)重要要求是，需要設(shè)計(jì)一個(gè)模擬區(qū)域能夠把盡可能多的震源和臺(tái)站包括其中.對(duì)于地震活動(dòng)不頻繁、但有密集地震臺(tái)陣的區(qū)域，常常無(wú)法只使用本地的震源來(lái)進(jìn)行全波形反演.一種可以在地殼和上地幔反演中同時(shí)考慮遠(yuǎn)震震源的方法是在三維反演區(qū)域之外，使用一些快速算法和較簡(jiǎn)單的背景模型來(lái)計(jì)算遠(yuǎn)震波場(chǎng)（比如FK 或者DSM 方法）.然后將計(jì)算所得到的波場(chǎng)當(dāng)作邊界條件，而在所要研究的區(qū)域三維模型中仍然使用數(shù)值算法（比如有限差分或者譜元法）來(lái)計(jì)算三維波場(chǎng)以及目標(biāo)函數(shù)梯度.這種混合數(shù)值模擬波場(chǎng)的計(jì)算方法使得我們可以在三維地殼波速反演以及反射面成像當(dāng)中使用遠(yuǎn)震記錄（Beller et al.,2017;Monteiller et al.,2015,2021;Pienkowska et al.,2021;Tong et al.,2014;Wang et al.,2016），從而大大地改進(jìn)密集臺(tái)陣下面的結(jié)構(gòu)分辨率以及成像覆蓋率.例如，Wang 等（2016）展示只要用極少數(shù)的幾個(gè)遠(yuǎn)震，并使用一個(gè)近線性臺(tái)陣上的P 波及其散射和折射波的記錄，就可以反演出歐洲西南部比利牛斯山脈（Pyrenees）下方巖石圈的縱波和橫波速度.相似的方法也被應(yīng)用到反演區(qū)域在震源和臺(tái)站中間的情況.從震源用快速算法和簡(jiǎn)單模型傳到模擬區(qū)域的地震波經(jīng)過(guò)三維數(shù)值模擬后，仍需要通過(guò)邊界條件用快速算法和簡(jiǎn)單地球模型傳回到接收臺(tái)站.這類混合數(shù)值計(jì)算波場(chǎng)的方法被用于全波形反演中，也被稱為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)對(duì)于我們了解板塊構(gòu)造、地球動(dòng)力學(xué)以及礦物物理具有十分重要的意義（Stern,2002）.目前在大陸尺度上，全波形反演能夠更好地刻畫海洋巖石圈在俯沖過(guò)程當(dāng)中所經(jīng)歷的復(fù)雜過(guò)程（Tao et al.,2018;Zhu et al.,2012）.相關(guān)的成像結(jié)果與傳統(tǒng)射線走時(shí)層析成像結(jié)果相當(dāng)接近.但同時(shí)通過(guò)提高模型分辨率以及更好地反演具體的波速大小，我們也發(fā)現(xiàn)了一些新的現(xiàn)象.比如Zhu 等（2012）在歐洲區(qū)域的研究中能夠清晰地看到板片滑脫（slab detachment），這一現(xiàn)象與以前走時(shí)層析成像的結(jié)果相似（Wortel and Spakman,2000），但在新的模型中我們能夠找到更多的區(qū)域具有此滑脫現(xiàn)象.另外，現(xiàn)實(shí)中海洋板塊在俯沖過(guò)程當(dāng)中會(huì)發(fā)生斷裂，而且這一斷裂面會(huì)沿著海溝方向不斷擴(kuò)大，進(jìn)而導(dǎo)致滑脫不斷增大.這些斷裂會(huì)導(dǎo)致軟流圈物質(zhì)上涌，進(jìn)而與地表的火山作用相關(guān).最近，Zhu 等（2020c）在中南美洲的全波形反演結(jié)果中發(fā)現(xiàn)Rivera 和Cocos 板塊在俯沖過(guò)程中也出現(xiàn)類似的滑脫現(xiàn)象.同時(shí)在板塊斷裂區(qū)間填充有大量的低速帶物質(zhì)，其中可能發(fā)生部分熔融現(xiàn)象.另一方面，通過(guò)反演方位各向異性，他們發(fā)現(xiàn)板塊的快速后退以及板塊之間的空隙會(huì)產(chǎn)生巨大的壓強(qiáng)差從而導(dǎo)致局部的軟流圈物質(zhì)流動(dòng).這一流動(dòng)會(huì)使得橄欖石等各向異性礦物的快速軸方向沿著流動(dòng)/應(yīng)變方向定向排布，從而產(chǎn)生可觀測(cè)的地震波各向異性.這對(duì)于我們研究俯沖板塊與周圍的地幔物質(zhì)相互對(duì)流作用具有十分重要的意義.另一方面，在西太平洋和東北亞，早期的射線走時(shí)層析成像已經(jīng)發(fā)現(xiàn)太平洋板塊在俯沖過(guò)程中與上地幔，特別是地幔轉(zhuǎn)換帶具有復(fù)雜的關(guān)系（Fukao and Obayashi,2013;Zhao et al.,1992,1994）.比如在東北亞，太平洋板塊在地幔轉(zhuǎn)換帶中可以橫向地延展上千千米.同時(shí)Tang 等（2014）的研究顯示在地幔過(guò)渡帶中，太平洋板塊發(fā)生斷裂，導(dǎo)致地幔物質(zhì)上涌，進(jìn)而與地表的長(zhǎng)白山五大連池火山相連.Tao 等（2018）通過(guò)使用全波形反演更好地描述這一俯沖板塊在上地幔，特別是地幔過(guò)渡帶中的特征.如圖11 所示，通過(guò)擬合地幔過(guò)渡帶的SH 波形，相對(duì)于以前的成像結(jié)果，他們的俯沖板塊分辨率更高，成像結(jié)果更加清晰.在全球尺度上，相關(guān)的伴隨層析成像結(jié)果也更好地顯示俯沖板塊在下地幔的具體形態(tài).而這些特征在傳統(tǒng)的走時(shí)層析成像結(jié)果當(dāng)中常常并不十分清楚，進(jìn)而導(dǎo)致模型解釋上的一系列爭(zhēng)論.通過(guò)更加清晰地描繪俯沖板塊在地幔過(guò)渡帶和下地幔的幾何學(xué)特征，能夠更好地幫助我們了解地幔對(duì)流以及地球內(nèi)部的動(dòng)力學(xué)過(guò)程（Lei et al.,2020）.

圖11 使用通過(guò)地幔轉(zhuǎn)換帶的SH 波改進(jìn)俯沖板塊伴隨層析成像結(jié)果.（a）顯示的是二維剖面和所使用的地震臺(tái)站的位置.（b）對(duì)比初始模型（左）和迭代改進(jìn)之后的模型（右）對(duì)SH 地幔轉(zhuǎn)換波的影響.黑色和紅色地震圖分別是實(shí)際觀測(cè)和模擬計(jì)算所得到的結(jié)果.（c）和（d）分別是初始模型和迭代改進(jìn)之后的模型在（a）中所示的二維剖面上的剪切波速度擾動(dòng)（修改自 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）地幔柱成像.相對(duì)于海洋板塊的俯沖，地幔柱提供了地幔對(duì)流的上升區(qū)域.從地幔柱假說(shuō)提出至今（Morgan,1971），地震學(xué)界對(duì)地幔柱的成像一直十分關(guān)注，而且不同的成像方法常常得到了不同結(jié)論（Montelli et al.,2004b,2006）.比如有的研究顯示冰島的地幔柱來(lái)自于核幔邊界（Bijwaard and Spakman,1999），而有的研究顯示這一區(qū)域的低速帶主要局限于上地幔（Ritsema et al.,1999）.相對(duì)于高波速的海洋板塊俯沖帶來(lái)說(shuō)，低波速的地幔柱對(duì)于成像精度和分辨率的要求更大.由于波前愈合以及成像分辨率的局限，傳統(tǒng)的走時(shí)層析成像很難反演幾百千米直徑的上升低速地幔柱.相反，全波形反演考慮了上述的物理過(guò)程，從而能夠更好地幫助我們反演尺度更小的低波速帶.比如Richers 等（2013）發(fā)現(xiàn)冰島地幔柱在上地幔中向東部偏移，可能與地幔對(duì)流方向有關(guān) .在全球尺度上，新的成像結(jié)果發(fā)現(xiàn)了更多的地幔柱，它們連接核幔邊界和地表的熱點(diǎn).French 和Romanowicz（2015）還發(fā)現(xiàn)這些地幔柱與軟流圈上垂直于洋中脊的低速帶相連.值得注意的是這一研究使用的是譜元方法來(lái)正演模擬地震波傳播，但是使用的是近似的方法（nonlinear asymptotic coupling theory,NACT）來(lái)計(jì)算敏感核.相對(duì)于伴隨層析成像和全波形反演，這樣能夠大大地減少計(jì)算量.如圖12 所示，Lei 等（2020）在全球伴隨層析成像模型當(dāng)中發(fā)現(xiàn)南非地幔柱在上升到地幔轉(zhuǎn)換帶時(shí)，分成多個(gè)更小的地幔柱，然后上升連接到地表的火山和熱點(diǎn).這些結(jié)果大大地提高了對(duì)地幔柱的認(rèn)識(shí)，從而為更好地了解地幔對(duì)流提供了巨大的幫助.

圖12 對(duì)比三個(gè)全球尺度地幔剪切波成像模型在6 個(gè)二維剖面上的結(jié)果.從左到右分別是來(lái)自于GLAD-M25（Lei et al.,2020）、TX2015（Lu and Grand,2016）和SEMUCB-WM1（French and Romanowicz,2015）.（a-f）分別對(duì)應(yīng)于以下的熱點(diǎn):（a）阿法爾州;（b）百慕大群島和加那利群島;（c）佛得角和哈加爾高原;（d）冰島和艾費(fèi)爾高原;（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 油氣儲(chǔ)層反演

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

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

（2）反射波全波形反演.早期的二維全波形理論實(shí)驗(yàn)表明，如果只使用透射波，比如井間層析成像，則可以很好地約束低波數(shù)信息，即光滑的速度結(jié)構(gòu).相反，如果只使用反射波（地表接收信號(hào)），則主要只能約束高波數(shù)的速度分布（Gauthier et al.,1986;Mora,1987）.如何使用反射波來(lái)約束低波數(shù)的速度模型一直以來(lái)是一個(gè)重要的研究問(wèn)題.通過(guò)構(gòu)建反射波的敏感核，可以使用反射波的到時(shí)以及波形來(lái)達(dá)到約束低波數(shù)的速度模型的效果.相對(duì)于直達(dá)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 在二維均勻速度模型下對(duì)比直達(dá)波和反射波的梯度.（a）和（b）分別是直達(dá)波和反射波的梯度，其中反射界面如圖（b）中的紫色直線所示.（c）直達(dá)波+反射波的梯度.（d）改進(jìn)的聯(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)中，我們使用層析成像來(lái)反演光滑的低波數(shù)速度結(jié)構(gòu).然后使用偏移方法，比如逆時(shí)偏移成像（McMechan,1983），來(lái)找到相應(yīng)的高波數(shù)反射界面以及油氣儲(chǔ)層位置.隨著全波形反演的發(fā)展，它被逐漸地推廣到高頻段.目前的反演頻率可以達(dá)到30～50 Hz.如果仍然能夠使這一非線性反演問(wèn)題得到收斂，而且避免前述的周期跳躍問(wèn)題，在未來(lái)有希望使用全波形反演來(lái)替代以前的層析成像+偏移的勘探開發(fā)步驟.

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

圖15 使用基于目標(biāo)函數(shù)梯度的馬爾科夫鏈蒙特卡羅方法反演二維Marmousi 模型.（a-c）分別表示平均模型、標(biāo)準(zhǔn)差模型以及平均與真實(shí)模型之間的差.（d）顯示6 個(gè)采樣點(diǎn)的一維邊界后驗(yàn)概率分布情況（修改自 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)域的全波形反演問(wèn)題涉及到前面所述的設(shè)計(jì)優(yōu)化的目標(biāo)函數(shù)，多參數(shù)聯(lián)合反演，與機(jī)器學(xué)習(xí)的結(jié)合以及在實(shí)際數(shù)據(jù)中提取低頻信號(hào)等，這里不再做贅述.

3 結(jié)論與展望

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

致謝

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