許忠淮

（中國地震局地球物理研究所，北京100081）

地震科普

地震學(xué)百科知識（三）
——地震學(xué)反演問題＊

許忠淮

（中國地震局地球物理研究所，北京100081）

引言

求解地震學(xué)反演問題是一般數(shù)學(xué)物理反演問題的分析方法在地震學(xué)研究中的應(yīng)用。一般反演問題的提出可表述如下。

當(dāng)人們研究一個物理對象時，首先要選定一些描述對象特征的參量χ，根據(jù)已知的物理規(guī)律列出χ變化所遵從的數(shù)量聯(lián)系，即遵從的方程

函數(shù)F的具體形式（如地震波的波動方程、地震波的走時方程等）由研究對象運動變化的具體規(guī)律決定。矢量χ＝（χ1，χ2，…，χL）T表示共有L個參量（上標(biāo)T表示轉(zhuǎn)置），函數(shù)矢量F＝（F1，F(xiàn)2，…，F(xiàn)N）T表示式（1）代表含N個方程的方程組：

對許多問題，人們還沒有完全認(rèn)識有關(guān)的規(guī)律，只是在觀測和實驗的基礎(chǔ)上提出表示各參量之間可能物理聯(lián)系的假設(shè)關(guān)系式來試驗性地分析問題，這時F就稱為描述認(rèn)識對象的模型。參量χ中包含兩部分，一部分是未知的構(gòu)建模型用的模型參數(shù)m；另一部分是可觀測的參數(shù)d。觀測量是會有誤差的，但在一定的認(rèn)識條件和認(rèn)識水平下，有些量的誤差可認(rèn)為很小，它們對問題分析結(jié)果的影響可忽略，這些量就被當(dāng)作模型中的常量?？紤]模型參數(shù)和觀測參數(shù)的劃分后，式（1）可改寫為

對有些實際問題，式（3）可以寫成以下觀測方程的形式［1］：

根據(jù)設(shè)定的模型關(guān)系式，對給定的模型參量，通過分析演繹，預(yù)測會有怎樣的觀測量的問題，稱為正演問題；由已知的觀測值，根據(jù)模型推斷可能有怎樣的模型參數(shù)值的問題稱為反演問題。

在求解反演問題之前，人們必須先能求解正演問題，即需要能理解產(chǎn)生觀測結(jié)果的物理過程，以便建立該過程的可靠的數(shù)學(xué)模型。

在確定模型F（或G）和已知數(shù)據(jù)d以后，反演問題就歸結(jié)為根據(jù)式（3）或式（4）求解模型參數(shù)m的數(shù)學(xué)問題。對不少實際問題，函數(shù)F（m，d）或G（m）是非線性函數(shù)，這時要求解的數(shù)學(xué)問題是個非線性反演問題；當(dāng)它們是線性函數(shù)時，則要求解的是線性反演問題。對于許多實際問題，這一求解過程常常是比較復(fù)雜的，目前已發(fā)展了許多求解反演問題的數(shù)學(xué)方法。

在地震學(xué)中，人們已經(jīng)研究了許多“正演問題”的解，即對給定的地球結(jié)構(gòu)和地震震源的模型，人們已計算出了可觀測的地震運動的特性，例如，各種地震波的走時、地震波的頻散曲線、地震引起的強(qiáng)地震動與近場地震動的頻譜、地球自由震蕩周期、遠(yuǎn)場地震波的波形和完整的地震圖，等等。“反演問題”就是將這些計算結(jié)果與實際觀測數(shù)據(jù)進(jìn)行比較，使計算結(jié)果能最好地擬合觀測數(shù)據(jù)，以便確定描述地球結(jié)構(gòu)和地震震源的各種模型參數(shù)。

1 近震定位問題

作為地震學(xué)反演問題的例子，現(xiàn)分析一個最簡單的近地震的定位問題。

設(shè)均勻水平地殼內(nèi)P波傳播速度vP為已知量。選取一原點o在地表、χ軸向北、y軸向東、z軸向下的直角坐標(biāo)系，設(shè)地殼內(nèi)在時刻τ（發(fā)震時刻）、深度為z處發(fā)生一個小地震，其震中位置的坐標(biāo)為（χ，y）。設(shè)地震周圍有方位不同、離地震遠(yuǎn)近不同的N個地震臺接收到地震發(fā)出的直達(dá)P波的到達(dá)時間各臺的空間位置（χsj，ysj，zsj）可認(rèn)為是已知量。根據(jù)所用的均勻地殼模型，各臺站的理論到時Tj應(yīng)等于發(fā)震時刻τ加上P波到各臺的走時：

式中Rj表示地震至臺站j的距離。我們求解的未知地震參數(shù)（τ，χ，y，z）應(yīng)能使各臺的計算到時與各臺的觀測到時相符：

可將上式的N個方程組用矢量形式表達(dá)為

式中m＝（τ，χ，y，z）T是4維的模型參數(shù)矢量，T（m）＝（T1（m），…，TN（m））T是N維的計算到時矢量，是N維的觀測到時矢量。式（7）所示的地震定位問題正是一般反演問題公式（4）的一個具體表達(dá)形式。式（6）或（7）可稱為地震定位問題的觀測方程組。由式（5）可見，各臺的計算走時Tj是模型參數(shù)的非線性函數(shù)，即T（m）代表一組非線性函數(shù)，因此，這里的地震定位問題是個非線性反演問題。

由于定位模型的近似性和觀測存在誤差等因素，一般是不能找到使式（6）或（7）嚴(yán)格成立的模型參數(shù)的。方程組（7）中未知數(shù)是4個，一般臺站數(shù)N都會多于4個，即式（7）是個超定方程組，無法求數(shù)學(xué)上嚴(yán)格成立的解答。求解反演問題通常是求能使殘差矢量

盡量小的模型參數(shù)。r的任一分量rj＝Tj（τ，是第j臺的計算到時與觀測到時的差異，即第j臺的到時殘差。我們所要求的解答是要使各臺的計算到時與觀測到時總體上能盡量符合。

可以有不同的標(biāo)準(zhǔn)來衡量計算到時與觀測到時的符合程度，一種常用方法是使以下的失配函數(shù)Q（亦稱目標(biāo)函數(shù)）取極小值：

選取不同的失配函數(shù)常會使最后的解答有些差異。

2 反演算法

有多種從式（9）求解答m＝（τ，χ，y，z）T的方法，大致可分為兩類。一類是只計算并比較失配函數(shù)值的搜索法，這類方法多用來求解非線性反演問題；另一類是需要計算失配函數(shù)對未知參數(shù)的導(dǎo)數(shù)的方法，或稱求導(dǎo)法，使用該法的前提是存在失配函數(shù)對未知量的一定階數(shù)的導(dǎo)數(shù)。

（1）搜索法也有多種。例如隨機(jī)嘗試法，亦稱蒙特卡羅法，即在（τ，χ，y，z）的可能取值范圍內(nèi)用很多組（數(shù)量可很大）隨機(jī)數(shù)構(gòu)成嘗試解答，由式（9）計算每組嘗試解的Q值，最后將與最小Q值對應(yīng)的那組解作為可接受的解答。又如應(yīng)用較廣的遺傳算法，這是模仿生物進(jìn)化過程而設(shè)計的一種不斷改變和更新一組嘗試模型參數(shù)、逐步找到逼近觀測結(jié)果的優(yōu)化模型參數(shù)的算法。

（2）求導(dǎo)法中簡單的一種是只求一階導(dǎo)數(shù)的方法，現(xiàn)以上述近震定位的例子說明如下。

首先要選定一組近似的初始模型參數(shù)m（0）＝（τ0，χ0，y0，z0）T。例如，在臺站分布較密時，可將初至波到時最早的臺站位置設(shè)為震中的初始位置，根據(jù)以往對該地區(qū)地殼地震的了解設(shè)定一個可能的震源深度值，進(jìn)而由各臺波的到時減去由初始震源位置計算的到時，可估算出一個平均的初始發(fā)震時刻。也可用其他方法（如隨機(jī)嘗試法等）選定初始模型參數(shù)。

由式（5）知，各臺的計算到時Tj是待求模型參數(shù)的非線性函數(shù)，這樣導(dǎo)致式（8）是一個難以直接求解的非線性方程組。實際求解的方法之一是將Tj（m）在初值m（0）處作泰勒展開，并只保留一階導(dǎo)數(shù)項［2］：

式中偏導(dǎo)數(shù)的下標(biāo)0表示取初值m（0）處的導(dǎo)數(shù)值；后式中用Ajk（k＝1，2，3，4）表示了前式中的4個偏導(dǎo)數(shù)值；后式中的δτ＝ττ0，等等。對各個Tj（m）如此作泰勒展開近似后，求解非線性方程組（7）的問題遂可轉(zhuǎn)化為求解近似的線性觀測方程組

的問題，式中χ＝（δτ，δχ，δy，δz）＝δm＝ m-m（0）是待求的對初始模型參數(shù)m（0）的修正量；矩陣

后項矩陣的下標(biāo)0表示矩陣內(nèi)各偏導(dǎo)數(shù)皆取初始模型參數(shù)m（0）處的數(shù)值。式（11）中的數(shù)據(jù)矢量b＝T（m（0））-d。與觀測方程（11）相應(yīng)的殘差矢量r＝Aχ-b。

用最小二乘法求式（11）問題的解答，就是要求使［2］

取極小值的解答χ（1）。分別求Q對χ的4個分量的偏導(dǎo)數(shù)，并令這4個偏導(dǎo)數(shù)的表達(dá)式等于0后，可得到

即

該方程組稱為線性方程組Aχ＝b的正規(guī)方程組，由于未知量是4個，矩陣ATA是個4×4的對稱方陣，其本征值都是非負(fù)實數(shù)。求矩陣ATA的逆矩陣即可得到地震參數(shù)的修正值

將初始模型參數(shù)加上這一修正量，即得待求的地震定位參數(shù)m（1）＝m（0）＋χ（1）。通過迭代法，每次將新求得的參數(shù)再當(dāng)作初始模型參數(shù)，重復(fù)上述求解過程，可求得第k次迭代的結(jié)果m（k）＝m（k-1）＋χ（k），直至m（k）與m（k-1）的差異已在容許誤差之內(nèi)，即得最終的地震定位結(jié)果。

上述所用的求導(dǎo)法只在式（10）的計算到時泰勒展開中保留了一階導(dǎo)數(shù)項，對于非線性較強(qiáng)的復(fù)雜問題，一階導(dǎo)數(shù)近似法常收斂很慢，有時甚至得不到合適解答；于是需要在泰勒展開中保留2階導(dǎo)數(shù)項，相應(yīng)地已發(fā)展了專門算法。

對于實際地震定位問題，需要考慮地殼地震波速的分層結(jié)構(gòu)，或波速隨深度連續(xù)變化的結(jié)構(gòu)，甚至還要考慮波速的橫向變化，這時理論到時的計算方法要遠(yuǎn)比式（5）表達(dá)的公式復(fù)雜得多，因而通常的地震定位問題在數(shù)學(xué)上也是個高度非線性的問題。

3 解的不唯一性

反演問題的解答一般都是不唯一的。這主要是由于觀測數(shù)據(jù)中缺少全部未知模型參數(shù)的足夠信息，有限的觀測結(jié)果造成了數(shù)學(xué)問題的欠定性。以求解式（11）所示的線性反演問題

為例，對上述簡化的地震定位問題，未知量χ含4個未知數(shù)，通常臺站數(shù)會多于4個，或有些臺可有不同震相的到時數(shù)據(jù)，即觀測方程數(shù)多于4個。該方程組是超定的，一般應(yīng)能求出使Aχ≈b的解答來。但實際問題常使形式上的超定方程組實際是欠定的。例如，當(dāng)缺少靠近地震震中的近距離臺站時，單靠震中距較大的臺站P波初動到時觀測數(shù)據(jù)不能將震源深度約束住，即觀測矢量b中缺少未知量深度的信息，在數(shù)學(xué)上表現(xiàn)為式（11）中有些線性方程是近似相關(guān)的，真正獨立的方程少于4個。又如，當(dāng)用分布在一個小區(qū)域中的臺網(wǎng)去測定在臺網(wǎng)外較遠(yuǎn)處的地震的位置時，若只用各臺的P波到時數(shù)據(jù)，考慮觀測誤差后，各臺的到時幾乎相當(dāng)于一個地震臺的觀測結(jié)果，這時待解方程組將是高度欠定的。

對其他未知模型參數(shù)很多的地震學(xué)反演問題，一組觀測值常常對一些參數(shù)是超定的，而對另一些參數(shù)卻是欠定的。例如，反演多臺記錄的多個地震的震相到時數(shù)據(jù)，以求某地區(qū)地震波速度局部變化的地震層析成像問題，需要將研究區(qū)分成許多小的區(qū)塊，未知數(shù)是各區(qū)塊的波速相對于全區(qū)平均速度模型（初始模型）的擾動量。由于臺站和地震位置的局限性，有些區(qū)塊可能沒有地震射線穿過，或只有稀少的、方向比較單一的射線穿過，于是到時觀測數(shù)據(jù)中就不含有這些區(qū)塊的波速信息，造成反演的數(shù)學(xué)問題對這些區(qū)塊的波速擾動量是欠定的。

此外，采用不同的計算方法，有時也會得到不大一樣的解答。例如，選擇不同的失配函數(shù)，有可能使我們所求的最優(yōu)化解答結(jié)果不同。對復(fù)雜的非線性反演問題，失配函數(shù)常常是模型參數(shù)的多極值函數(shù)，當(dāng)用搜索法求解時，有時會將與失配函數(shù)局部極小值（而不是全局極小值）對應(yīng)的模型參數(shù)確定為解答，這也會造成解答的不唯一。

以上所述是指對給定的物理模型，在求解反演的數(shù)學(xué)問題中出現(xiàn)的解答的不唯一性。物理模型的不確定性也會使問題的解答含有不定性。例如，基于初步認(rèn)識，人們對地震定位問題選擇了速度均勻的地殼模型，如果另一些研究人員采用了地殼分層的速度模型，則地震的定位結(jié)果將會有所不同。這種模型的近似性可能會影響與之相關(guān)的反演的數(shù)學(xué)問題是否適定，特別是不符合實際的模型是難以從實際觀測數(shù)據(jù)中找到合理的支持證據(jù)的；但反演問題的不適定性主要是從分析數(shù)學(xué)問題中來討論的。

4 算法的選擇

在求解反演問題時，選擇合適的計算方法是非常重要的。

例如，對于經(jīng)常遇到的待解的線性方程組Aχ＝b，當(dāng)A是滿秩的方陣時，理論上可通過求A的逆矩陣而直接得到解答χ＝A-1b。但對于許多實際反演問題，矩陣A的階數(shù)常常是非常大的，這種情況下，即使A是方陣，也不能用經(jīng)典線性代數(shù)中的克萊姆法則求其逆矩陣A-1，因為這需要耗費大量機(jī)時，甚至不可能算出。但使用高斯消元法則可很快求出解答。這說明解反演問題需要使用合適的算法。

即使能順利求出Aχ＝b解答χ，也還有解答是否穩(wěn)定的問題。例如，對于問題［2］

在實際反演問題中，A通常是行數(shù)大于列數(shù)的矩陣，無法使用消元法之類的解法。在上述地震定位問題中采用了最小二乘法來求解χ，導(dǎo)致求解正規(guī)方程組（14）而得到解答式（15）。這一算法的優(yōu)點是求解的式（14）是一個低階的適定方程組，矩陣ATA是個對稱方陣，其本征值全是非負(fù)實數(shù)，完全可以求逆矩陣。但這一算法也有明顯缺點：①用計算機(jī)解正規(guī)方程時，計算精度需為求解原始方程組的兩倍；② 組成矩陣ATA和矢量ATb時，破壞了原始方程中的某些信息。

5 奇異值分解法

對于一般的線性方程組

式中各量下標(biāo)表示行數(shù)乘列數(shù)，由于n≠m，A不是滿秩方陣，因而無法求逆矩陣A-1。數(shù)學(xué)家蘭克卓斯（C.Lanczos）提出了一種A的奇異值分解方法［1］，由此可計算出A的廣義逆矩陣。奇異值分解法將矩陣A分解為以下3個矩陣的乘積：

式中，p是A的不為零的奇異值的個數(shù)；矩陣S是由p個奇異值（s1，…，sp）構(gòu)成的對角矩陣；矩陣Up由p個正交單位矢量（u1，…，up）構(gòu)成，每個矢量含n個分量；矩陣Vp由p個正交單位矢量（v1，…，vp）構(gòu)成，每個矢量含m個分量。這一分解方法的數(shù)學(xué)推導(dǎo)可參閱專業(yè)書籍，例如《定量地震學(xué)》12.3節(jié)。早有人已編制好奇異值分解法的計算程序，下面給出一個具體A矩陣的奇異值分解實例，可見由右側(cè)3個矩陣的乘積確實可得到左側(cè)矩陣的具體數(shù)值，并可見到矩陣Up和Vp列向量的相互正交性（注意下式 最右端矩陣的行向量即是Vp的列向量）。

數(shù)學(xué)上有人提出了一種關(guān)于廣義逆矩陣的定義如下。

設(shè)A為n×m的實矩陣，若某個n×m的矩陣H滿足以下條件：

則將這樣的矩陣H定義為A的穆爾-彭羅斯（Moore-Penrose）廣義逆矩陣。根據(jù)A矩陣的奇異值分解表示式（17），并注意

Ip是p×p的單位矩陣，很容易驗證矩陣

是滿足式（18）中的諸條件的，因而A-1g就稱為矩陣A的廣義逆矩陣。上式中的矩陣S-1是由矩陣A的p個奇異值（s1，…，sp）的倒數(shù)構(gòu)成的對角方陣。

于是可求得方程組（16）中模型矢量χ的廣義逆解為

在求解線性方程組Aχ＝b時，當(dāng)不存在A矩陣的精確逆矩陣A-1時，就可取式（21）表達(dá)的廣義逆解作為問題的解答。下面看看該解答的一些特性。

（1）分辨率矩陣［1］

根據(jù)式（16）有b＝Aχ，代入式（21）并考慮式（17）對A的分解，可得

（2）解的誤差估計［2］

觀測數(shù)據(jù)的誤差會帶來解答的誤差。待求解的方程組實際是Aχ＝（b±Δb），這里假定A是精確的。對于廣義逆解，已有數(shù)學(xué)家證明［3］，由b的誤差Δb引起的χ的誤差Δχ可由下式給出

式中λ＝smax／smin是矩陣A的最大和最小奇異值的比值，并稱為矩陣A的條件數(shù)；‖·‖表示矢量的范數(shù)，等于各分量平方和的平方根。式（23）表明，A的條件數(shù)λ是觀測量相對誤差放大量的上界。式（23）說明，解的誤差不但與觀測誤差有關(guān)，也與模型矩陣A的結(jié)構(gòu)有關(guān)。

若矩陣A條件數(shù)過大（如＞104），A將有顯著的奇異性，即不能求得穩(wěn)定的解答。對某個實際反演問題，如果矩陣A出現(xiàn)零奇異值（這里的“零”是數(shù)值計算意義上的零，即雖有很小的數(shù)值，但已與計算誤差無法區(qū)別），則不能用上述方法求得廣義逆解，需用另法（例如阻尼最小二乘法）求取反演問題的可能解答。

（作者電子信箱，許忠淮：xuzh＠cea-igp.ac.cn）

［1］［美］安藝敬一，P G理查茲.定量地震學(xué)：理論和方法.北京：地震出版社，1987

［2］［美］李汯鑒，S W斯圖爾特.微震臺網(wǎng)的原理及應(yīng)用.唐美華，譯；葉世元，趙仲和，校.北京：地震出版社，1984

［3］Forsythe G E，Malcolm M A，Moler C B.Computer Methods for Mathmatical Computations.Englewood Cliffs，New Jersey：Prentice-Hall，1977

P315.01；

A；

10.3969／j.issn.0235-4975.2013.04.009

2012-11-07。