孫福玉

摘 要：本文對地球物理學(xué)反演的意義做了簡要說明。以地球物理反問題中的地震波褶積模型為例，通過理論分析和數(shù)學(xué)推導(dǎo)給出超定問題的最小二乘解法、先驗(yàn)信息的應(yīng)用、線性反問題的廣義反演法等多種常用方法，并通過matlab編程實(shí)現(xiàn)做最終成圖對比分析。

關(guān)鍵詞：超定問題;matlab編程;最小二乘法;先驗(yàn)信息;廣義反演

中圖分類號(hào)：P631.4 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1671-2064（2020）06-0186-02

0 引言

在自然界中，存在大量的客觀事物，這些客觀事物在傳統(tǒng)物理學(xué)中，常?？梢钥闯墒且粋€(gè)系統(tǒng)。這些事物包羅萬象，小到原子大到宇宙。一般來說，所有的客觀事物都能用一定的物理量來描述，但不是所有的物理系統(tǒng)所包含的信息都能夠被直接測量，有一部分物理系統(tǒng)的物理量只能通過間接測量。在地球物理學(xué)領(lǐng)域的研究中，研究對象絕大多數(shù)都位于地表以下。例如，研究地下地層的構(gòu)造特征、尋找含油儲(chǔ)集層、尋找礦產(chǎn)資源、尋找含水層。將地下介質(zhì)全部挖開去直接研究這些地質(zhì)信息顯然是不可能的。為了全面的了解一個(gè)物理系統(tǒng)，往往希望獲得盡可能多的描述系統(tǒng)的物理量，往往能夠直接獲取的物理量十分有限。作為一種演繹性的研究手段，地球物理反演就是通過不斷地推演系統(tǒng)中那些可以被直接觀測的物理量與不可直接測量的物理量之間的聯(lián)系。達(dá)到通過已知信息去探索未知信息的一種方法。

反演問題一直以來都是地球物理學(xué)中的核心理論問題。由于人類自身的局限性，獲取深地信息的來源大都來源于地球表面，有極小一部分資料來自于對地球的直接探測，例如鉆井信息。所以人類對于地球的認(rèn)識(shí)和研究終究還是要依靠有限的地球物理信息來實(shí)現(xiàn)，這些信息如何應(yīng)用到實(shí)際研究中，就是反演問題的科研本質(zhì)。本文通過將理論方法用matlab編程實(shí)現(xiàn)，以反演地震子波為例，模擬地球物理反演過程。

1 超定問題的常規(guī)解法[1]

通常，用物理模型去求解物理參數(shù)的過程叫做正演，而給定一系列物理模型的參數(shù)信息去求解物理模型的過程叫做反演。反問題的研究前提是正演問題得到合理的解釋。正演問題大體可分為兩方面的研究，一是通過實(shí)驗(yàn)的方法尋找出經(jīng)驗(yàn)公式，利用此公式建立從模型到數(shù)據(jù)的聯(lián)系。二是通過嚴(yán)密的物理推導(dǎo)，通過給出的模型參數(shù)預(yù)測觀測數(shù)據(jù)。由此可見，只要正演問題不能被正確的解決，反演問題的研究就不能展開。此外，地下介質(zhì)環(huán)境異常復(fù)雜，及時(shí)正演問題即物理理論已經(jīng)被完善的解決，反演問題還是不能得出正確答案，需要考慮多種假設(shè)條件，通常會(huì)添加一些限制條件得出某種意義下的最佳解。

對于線性離散問題，總可以描述為d=GM的形式，可以理解為一系列的線性方程組的通式，其中為已知的觀測數(shù)據(jù)為待求的模型參數(shù)。 為數(shù)據(jù)核矩陣，與模型有關(guān)。

設(shè)觀測數(shù)據(jù)的數(shù)量為M，待定的模型參數(shù)數(shù)量是為N，G為M×N階矩陣，其秩為r，當(dāng)M>N=r時(shí)，觀測數(shù)據(jù)提供了多于模型參數(shù)數(shù)目的信息，此問題稱為超定問題，求解不同類型的線性反演問題，所采用的方法是有不同的。在本文中著重探討一下，超定問題的不同解法。

超定問題是觀測數(shù)據(jù)提供了多于模型參數(shù)數(shù)目的信息，因此求得的每一個(gè)G×M都會(huì)和真實(shí)觀測數(shù)據(jù)存在誤差。要想尋求超定問題真解是不可能的，只能尋求其在不同條件下的最佳解。在這里通常是去構(gòu)造最小誤差解，也就是說尋找到這樣一組模型參數(shù)，它和G核矩陣相乘所得的計(jì)算值和真正的觀測數(shù)據(jù)之間的差異最小。這個(gè)差異通常用誤差的二范數(shù)平方來表征。二范數(shù)平方的最小值解就是常用的最小二乘解。

所以對于超定問題的求解方案可以這樣描述：期望求得一組與觀測數(shù)據(jù)誤差平方和為最小的預(yù)測數(shù)據(jù)所對應(yīng)的模型參數(shù)[3]。

為了不失普適性，我們先假設(shè)觀測數(shù)據(jù)d為M維向量，模型參數(shù)m為N維向量，且有M>N，則上式的求解可轉(zhuǎn)化為一個(gè)線性方程組的求解[3]，即：

通過目標(biāo)函數(shù)對模型參數(shù)求一階導(dǎo)函數(shù)可求得它的極值點(diǎn)，不難驗(yàn)證這個(gè)極值點(diǎn)即為極小值點(diǎn)，也就是最小值點(diǎn)。解得：

下面看一個(gè)地球物理學(xué)上的實(shí)際問題，地震波褶積模型。由于地震子波是延續(xù)時(shí)間非常短的信號(hào)，地震記錄是數(shù)據(jù)相對豐富的信號(hào)。那么，利用反射系數(shù)（測井資料）和地震記錄（井旁道地震資料）求地震子波這一反演問題，就成為經(jīng)典的超定問題。應(yīng)用matlab編程[2]實(shí)現(xiàn)這一實(shí)際問題的步驟如下：

（1）定義模型。設(shè)計(jì)一個(gè)地質(zhì)模型，定義速度、密度、最小偏移距、采樣間隔、道數(shù)等物理量。通過數(shù)學(xué)公式求出雷克子波，一定要注意雷克子波是一個(gè)延續(xù)時(shí)間非常短的信號(hào)，也就是十分短的向量。

（2）模型正演。因?yàn)榉囱輪栴}的求解都是建立在正演問題充分解決的基礎(chǔ)上的。所以要先求出正演數(shù)據(jù)。通過此模型生成反射系數(shù)序列，利用反射系數(shù)序列求出反射系數(shù)矩陣。正演實(shí)質(zhì)上是反射系數(shù)與雷克子波向量的褶積。所以在此問題中的G核矩陣實(shí)際上就是反射系數(shù)矩陣。

（3）模型反演。利用最小二乘法解出地震子波并將其與雷克子波向量做對比。如圖1為引入加權(quán)因子的反演結(jié)果對比。

2 先驗(yàn)信息的應(yīng)用[1]

上文中的解決方案，是在觀測數(shù)據(jù)基礎(chǔ)上去反推模型參數(shù)。實(shí)際工作過程中要想全面確定模型參數(shù)大多數(shù)情況下都是信息不足的，為了使反演問題的解更切合實(shí)際情況，我們需要添加一些在觀測數(shù)據(jù)中未包含的信息，這種附加給反演問題的信息叫先驗(yàn)信息。先驗(yàn)信息實(shí)際上是對觀測數(shù)據(jù)的一種補(bǔ)充。為了使得上文中超定問題的最小二乘解更加真實(shí)，則可以通過添加先驗(yàn)信息來進(jìn)行優(yōu)化。先驗(yàn)信息的類型有許多種，但對于不同的形式，都要使得目標(biāo)函數(shù)以定量的形式出現(xiàn)，并且是不依賴于實(shí)際數(shù)據(jù)的。

在本文中介紹超定問題常用的先驗(yàn)信息，引入加權(quán)因子。通過對觀測數(shù)據(jù)的預(yù)測誤差進(jìn)行加權(quán)來度量。在實(shí)際工作中，經(jīng)常會(huì)碰到某些觀測數(shù)據(jù)的精度比另外一些觀測數(shù)據(jù)的精度高的情況發(fā)生。在這種情況下，通常的解決方式是在總誤差的定量表示中，在比較精確的觀測數(shù)據(jù)所相應(yīng)的預(yù)測誤差上加的權(quán)大于不精確的觀測數(shù)值所加的權(quán)。

為了完成這樣的加權(quán)，通常要定義一個(gè)廣義預(yù)測誤差，為觀測數(shù)據(jù)的加權(quán)因子，確定了每個(gè)單獨(dú)的誤差對總預(yù)測誤差的相對貢獻(xiàn)。加權(quán)因子是一個(gè)對角矩陣，對角線上的元素，就是相應(yīng)觀測數(shù)據(jù)誤差所加權(quán)重的大小度量。因此可以通過廣義預(yù)測誤差最小來估計(jì)此時(shí)的模型參數(shù)。通過和上文類似的導(dǎo)函數(shù)推導(dǎo)即可求得解的形式。

加權(quán)因子的引入是提高某些數(shù)據(jù)的可信度，在地震波褶積模型問題中因?yàn)樗械臄?shù)據(jù)都足夠準(zhǔn)確，在此基礎(chǔ)上反演得到的模型不會(huì)因?yàn)榧訖?quán)因子的改變而改變，不管引入什么樣的加權(quán)因子都可以反演得到雷克子波。因此在本文中為了方便讀者對加權(quán)因子的引入有更加深入的理解，在matlab編程[2]時(shí)更改地震記錄（觀測數(shù)據(jù)）后引入加權(quán)因子。引入如圖2所示的加權(quán)因子，目的是削弱錯(cuò)誤數(shù)據(jù)的權(quán)重。

We=eye（600）;We=（200，200）=0.2;We=（290，290）=0.1;

如圖2所示，引入加權(quán)因子后削弱了錯(cuò)誤數(shù)據(jù)的權(quán)重，可以明顯地看出反演出來的模型更加接近雷克子波。

3 線性反問題的廣義反演[1]

首先本文先介紹什么是廣義逆，設(shè)矩陣存在矩陣滿足一部分的矩陣稱為G的廣義逆。可見，G的廣義逆不唯一共有15種。同時(shí)滿足所有條件的廣義逆，稱為Moore-Penrose逆，記為G+。對于任意矩陣的Moore-Penrose廣義逆G+必然存在，并且是唯一的。這樣的性質(zhì)才使得我們討論問題有意義。

本文著重討論廣義逆的求法之一，奇異值分解方法。任意一個(gè)實(shí)數(shù)矩陣A都可以寫成A×A轉(zhuǎn)置的特征向量構(gòu)成的列向量乘上A×A轉(zhuǎn)置特征值的平方根構(gòu)成的矩陣乘上A轉(zhuǎn)置×A的特征向量構(gòu)成的列向量三個(gè)向量乘積的形式。

圖中深色部分是奇異值已經(jīng)它所對應(yīng)的向量，由奇異值分解可以求得它的廣義M-P逆，它的形式如圖3。

通過公式推導(dǎo)不難得出超定、欠定、混定三種問題解的統(tǒng)一的表達(dá)式，也從另一個(gè)方面證明了解的統(tǒng)一性，實(shí)際上這都是基于2-范數(shù)。通過matlab編程[2]實(shí)現(xiàn)模型構(gòu)建、G核矩陣的求取，以上兩個(gè)步驟在上文已經(jīng)給出，接下來運(yùn)用matlab預(yù)定義的SVD函數(shù)求取A矩陣的奇異值分解，進(jìn)而求出A矩陣的廣義M-P逆矩陣。運(yùn)用廣義M-P逆的核矩陣帶入通用公式即可求解。

如圖4所示，運(yùn)用廣義M-P逆反演出的子波與原雷克子波模型完全一致。

參考文獻(xiàn)

[1] 姚姚.地球物理反演基本理論與應(yīng)用方法[M].武漢：中國地質(zhì)大學(xué)出版社，2002.

[2] 童孝忠，柳建新.MATLAB程序設(shè)計(jì)及在地球物理中的應(yīng)用[M].長沙：中南大學(xué)出版社，2013.

[3] 謝萬學(xué).井地聯(lián)合反演方法研究[D].青島：中國石油大學(xué)，2007.