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

孫福玉

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

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

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

0 引言

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

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

1 超定問題的常規解法[1]

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

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

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

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

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

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

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

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

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

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

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

2 先驗信息的應用[1]

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

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

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

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

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

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

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

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

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

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

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

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

參考文獻

[1] 姚姚.地球物理反演基本理論與應用方法[M].武漢：中國地質大學出版社，2002.

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

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