基于matlab的地震資料反演模擬

孫福玉

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

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

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

0 引言

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

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

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

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

對于線性離散問題，總可以描述為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ù)據(jù)提供了多于模型參數(shù)數(shù)目的信息，此問題稱為超定問題，求解不同類型的線性反演問題，所采用的方法是有不同的。在本文中著重探討一下，超定問題的不同解法。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

參考文獻

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

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

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