Pytorch實現波阻抗反演

GeoFXR發表於2022-06-22

Pytorch實現波阻抗反演

1 引言

地震波阻抗反演是在勘探與開發期間進行儲層預測的一項關鍵技術。地震波阻抗反演可消除子波影響,僅留下反射係數,再通過反射係數計算出能表徵地層物性變化的物理引數。常用的有道積分、廣義線性反演、稀疏脈衝反演、模擬退火反演等技術。

Pytorch實現波阻抗反演

隨著勘探與開發的深入,研究的地質目標已經從大套厚層砂體轉向薄層砂體,而利用常規波阻抗反演方法刻畫薄層砂體不僅要消耗大量人力、物力,且反演得到的波阻抗精度也難以滿足實際需求。近年來,深度學習在地震反演和解釋等地震領域顯現出了巨大的潛力,其中,卷積神經網路(Convolution Neural Networks,CNN)是解決地震反演的一大有力工具。本文將介紹如何搭建CNN網路進行波阻抗反演。

2 Python環境配置

我們將使用Python語言以及Pytorch深度學習框架完成實驗,因此首先需要配置Python環境。

(1)安裝Anaconda

下載地址Anaconda | The World's Most Popular Data Science Platform

Anaconda是一個開源的Python和R語言的發行版本,用於計算科學(資料科學、機器學習、大資料處理和預測分析),致力於簡化軟體包管理系統和部署。

安裝過程“下一步、下一步”,在這一步時,勾選新增Anaconda到環境變數(注意:在環境變數中有QT、R語言等時,須在“編輯環境變數”中手動將Anaconda路徑下移,防止變數、路徑被覆蓋)

Pytorch實現波阻抗反演

(tips:Python和Anaconda都有蟒蛇的意思)

選擇Anaconda主要有幾個原因:

  • Anaconda的包使用軟體包管理系統Conda進行管理,提供了conda install xxx命令,與pip install xxx功能相同。當pip安裝庫出現一些未知錯誤時,Conda可以作為補充,省去非必要的報錯糾結;
  • Anaconda整合了本機中的Python環境與Python包路徑,並且提供虛擬環境管理,使得不同Python版本(Python2.x,Python3.x)、CPU/GPU版第三方庫可以並列共存。避免了在今後使用中出現版本、路徑、環境混亂的問題;
  • 下載安裝Anaconda可以一鍵獲取新手套裝,無需從Python官網從零開始。

隨Anaconda一同下載的有

  • Anaconda Navigator:是包含在Anaconda中的圖形使用者介面,使用者可以通過Anaconda Navigator啟動應用,在不使用命令列的情況下管理軟體包、建立虛擬環境和管理路徑。

  • Jupyter Notebook:是一個基於Web的互動式計算環境,用於建立Jupyter Notebook文件。這類文件是一個JSON文件,包含一個有序的輸入/輸出單元格列表,這些單元格可以包含程式碼、文字(使用Markdown語言)、數學、圖表和富媒體 (Rich media),通常以“.ipynb”結尾擴充套件。Jupyter Notebook的最大的優點就是可以“做一步看一步”,對於初學和開荒Python專案時,可以更簡單、高效地程式設計。

Pytorch實現波阻抗反演
  • Spyder,是一個開源、免費的Python整合開發環境(IDE)。其最大優點是模仿MATLAB的“工作空間”的功能,便於觀察變數的值、維度、型別等資訊。

  • Anaconda Prompt:等同Window PowerShell和cmd。

(2)安裝Pycharm

下載地址Download PyCharm: Python IDE for Professional Developers by JetBrains

由捷克JetBrains公司專為Python打造,PyCharm具備一般IDE的功能,比如,除錯、語法高亮、Project管理、程式碼跳轉、智慧提示、自動完成、單元測試、版本控制等等。PyCharm是商業軟體,與Visual Stidio類似,有收費的專業版(Professional)和免費的社群版(Community),社群版能夠滿足一般程式設計需要,專業版整合了Jupyter Notebook並且支援遠端開發。

下載Community版後,安裝過程除了調整路徑外,沒有需要操作的地方,非必須的可選項都不打勾,然後“下一步下一步”。

安裝完成後來到新建工程,以及配置Python環境的環節:

Pytorch實現波阻抗反演

第一次建立工程、載入環境需要等上幾分鐘,讀條完成後,這個工程的Python環境及編輯器就配置完成了。

3 Pytorch實戰

首先我們需要在本地安裝Pytorch,開啟Anaconda Prompt或Window PowerShell或cmd,安裝指令為

pip install torch torchvision torchaudio 
或者
conda install pytorch torchvision torchaudio cpuonly -c pytorch 

即可安裝CPU版Pytorch(GPU版的安裝稍複雜,此處不作詳細介紹),隨後開始本次專案的實戰部分。

(1)資料準備

Step 1 匯入資料

匯入的資料包括地震記錄,此處的資料已經被整理儲存為.mat格式,通過scipy.io庫將其讀入。地震道與波阻抗資料都是檔案中的鍵值對,從原始資料中抽取出部分作為實驗資料集。

dataframe = sio.loadmat('Train_DataSyn_Ricker30.mat')

#從檔案中分別提取提取地震道與波阻抗資料
Seismic_data = dataframe['Seismic']				#地震道
Impedance_data = dataframe['Imp']/1e6		#波阻抗

#隨機抽取部分道,作為訓練集
howMany = 2020 
np.random.seed(9)	#隨機種子,便於復現
indxRand = [randint(0,dataframe['Seismic'].shape[1]-1) for p in range(0,howMany)]	#隨機索引

#地震道
Seismic_data = Seismic_data.transpose()	#轉置
Seismic_data = Seismic_data[indxRand,:]	#通過索引抽取
#波阻抗
Impedance_data = Impedance_data.transpose()	
Impedance_data = Impedance_data[indxRand,:]

資料展示如下:

Pytorch實現波阻抗反演

我們的目標即是通過建立CNN模型挖掘規律,建立地震振幅屬性\(\Longrightarrow\)波阻抗的對映,在未知波阻抗的地方可用地震記錄進行預測,實現波阻抗反演。

Step 2 分割資料集

通常一個機器學習專案會需要我們對資料集進行分割,劃分為訓練集(建模)、驗證集(調整超引數與初步評估)和測試集(評估模型)。

#其中驗證集500個,測試集1000個,剩餘(520)為訓練集
howManyToValidate = 500
howManyToTest = 1000

#對輸入Seismic與標籤Imp進行相同的索引與處理
#用numpy索引切片的方式進行劃分
valX = (Seismic_data[:howManyToValidate,:])
testX = (Seismic_data[howManyToValidate:howManyToValidate+howManyToTest,:])
trainX = (Seismic_data[howManyToValidate+howManyToTest:,:])

valImp = (Impedance_data[:howManyToValidate,:])
testImp = (Impedance_data[howManyToValidate:howManyToValidate+howManyToTest,:])
trainImp = (Impedance_data[howManyToValidate+howManyToTest:,:])

#轉為torch中的Tensor格式
#此時資料為(道數,取樣點數)的二維陣列,按照torch的輸入格式整理為(道數,資料高度,資料長度),便於後續輸入道CNN網路中
valX = torch.FloatTensor(np.reshape(valX, (valX.shape[0], 1, valX.shape[1])))
testX = torch.FloatTensor(np.reshape(testX, (testX.shape[0], 1, testX.shape[1])))
trainX = torch.FloatTensor(np.reshape(trainX, (trainX.shape[0], 1, trainX.shape[1])))

valImp = torch.FloatTensor(np.reshape(valImp, (valImp.shape[0], 1, valImp.shape[1])))
testImp = torch.FloatTensor(np.reshape(testImp, (testImp.shape[0], 1, testImp.shape[1])))
trainImp = torch.FloatTensor(np.reshape(trainImp, (trainImp.shape[0], 1, trainImp.shape[1])))

(2) 建模

Step 1 建立網路

採用一維卷積神經網路(1D CNN)進行波阻抗的生成與預測,具體步驟與結構如下圖。

Pytorch實現波阻抗反演
建立神經網路通常也就是建立一個繼承自***torch.nn.Module***的類,而將網路層及其連線定義在類中的方法中
noOfNeurons = 60		#定義卷積核個數
dilation = 1			
kernel_size = 300		#卷積核尺寸
stride = 1				#卷積核滑動步長
padding = int(((dilation*(kernel_size-1)-1)/stride-1)/2)	#0填充個數

class CNN(nn.Module):					
    def __init__(self):					#建構函式
        super(CNN, self).__init__()		#前面三行為固定格式
        self.layer1 = nn.Sequential(	#nn.Sequential為一個順序的容器
            #Conv1d的前兩個引數分別表示輸入通道數與輸出通道數,又有卷積核個數=卷積層輸出通道數
            nn.Conv1d(1, noOfNeurons, kernel_size=kernel_size, stride=1, padding = padding+1),#卷積層
            nn.ReLU())	#ReLu啟用函式

        self.layer2 = nn.Sequential(
            nn.Conv1d(noOfNeurons, 1, kernel_size=kernel_size, stride=1, padding = padding+2),
            nn.ReLU())

    def forward(self, x):	#在forward中將網路像搭積木一樣連線起來
        out = self.layer1(x)	
        out = self.layer2(out)

        return out

cnn = CNN()		#例項物件

注:①在PyTorch中,可以把神經網路類中forward函式看作一個專用函式,它專門用於編寫前向傳播的計算方法,並且已經寫在了nn.Module的“出廠設定”中,在傳入資料時便會開始執行,如上例中,使用cnn(x)本質上等於cnn.forward(x),顯式使用後者反而報錯。

②關於啟用函式

啟用函式給神經元引入了非線性因素,使得神經網路可以任意逼近任何非線性函式,能夠完成極其複雜的分類或迴歸任務。若沒有啟用函式,則每層就相當於矩陣乘法。

Sigmoid——nn.Sigmoid()

\(\text{Sigmoid}(x) = \sigma(x) = \frac{1}{1 + \exp(-x)}\)

優點:

能夠將函式壓縮至區間[0,1]之間,保證資料穩定,波動幅度小

缺點:

  • 函式在兩端的飽和區梯度趨近於0,當反向傳播時容易出現梯度消失或梯度爆炸
  • 輸出不是0均值(zero-centered),導致收斂緩慢
  • 運算量較大

如果我們在多個層中堆疊sigmoid,則系統學習效率可能低下,並且需要仔細初始化。因此,對於深度神經網路,首選ReLU函式。

ReLU(Rectified Linear Unit,修正線性單元)——nn.ReLU()

\(\text{ReLU}(x) = (x)^{+} = \max(0,x)\)

優點:

  • 梯度不飽和,收斂速度快
  • 減輕反向傳播時梯度彌散的問題
  • 無需進行指數運算,運算速度快、複雜度低

缺點:

  • 輸出不是0均值(zero-centered)
  • 對引數初始化和學習率非常敏感,容易出現神經元死亡。

常用的啟用函式還有許多:

Pytorch實現波阻抗反演
**Step 2 定義訓練引數**
#定義與訓練有關的超引數
num_epochs = 500		#迭代輪數
batch_size = 1			#批次大小
learning_rate = 0.001	#學習率
batch_size_tot = trainX.shape[0]	
no_of_batches = int((batch_size_tot - batch_size_tot%batch_size)/batch_size)	#總批數

loss_fn = nn.MSELoss()	#損失函式,此處為一個迴歸任務,採用均方根誤差作為損失值
optimizer = torch.optim.Adam(cnn.parameters(), lr=learning_rate)	#優化器選擇為Adam

注:①pytorch中的nn模組提供了很多可以直接使用的loss函式,常用的有:

損失函式 名稱 適用場景
torch.nn.MSELoss() 均方誤差損失 迴歸
torch.nn.L1Loss() 平均絕對值誤差損失 迴歸
torch.nn.CrossEntropyLoss() 交叉熵損失 多分類
torch.nn.BCELoss() 二分類交叉熵損失 二分類
torch.nn.MultiLabelMarginLoss() 多標籤分類的損失 多標籤分類

②使用Loss函式會對每個樣本計算其損失,然後開始梯度下降去降低損失\(w+=-\alpha dx\),最基礎的是一次優化輸入一個樣本的隨機梯度下降(Stochastic gradient descent,SGD),或是一個小批次(Mini-batch)。除此之外還有很多種不同的優化器,torch.optim是一個實現了各種優化演算法的庫。它們採取不同的策略更新梯度,比初始的梯度下降更加高效:

函式 描述 公式
torch.optim.SGD SGD演算法 \(w=w-{\frac {\eta }{n}}\sum _{i=1}^{n}\nabla J_{i}(w)\)
torch.optim.SGD
(Momentum=0.9...)
Momentum演算法 \(m=b_1*m-\alpha dx\)
\(w=w+m\)
torch.optim.AdaGrad AdaGrad演算法 \(v=v+dx^2\)
\(w=w-\alpha *\frac{dx}{\sqrt v }\)
torch.optim.RMSProp RMSprop演算法
(Root Mean Square Prop)
Momentum+AdaGrad
\(v=b_1*v+(1-b_1)*dx^2\)
\(w=w-\alpha *\frac{dx}{\sqrt v }\)
torch.optim.Adam Adam演算法
(Adaptive momentum)
\(m=b_1*m+(1-b_1)*dx\)
\(v=b_2*v+(1-b_2)*dx^2\)
\(w=w-\alpha *\frac{m}{\sqrt v }\)

Step3 訓練

for epoch in range(num_epochs):		#開始迭代
    for i in range(no_of_batches):
     	#地震道資料
        #通過手動索引的方式建立batch
        #使用Variable對Tensor進行封裝,便於改變Tensor的.data、.grad、.grad_fn屬性
        traces = Variable(trainX[i*batch_size:(i+1)*batch_size,:,:])
        imp_label = Variable(trainImp[i*batch_size:(i+1)*batch_size,:,:])

        '''以下5行程式碼為固定格式,幾乎所有pytorch網路都是同樣的'''
        outputs = cnn(traces)	#將訓練資料輸入道cnn模型中,前向傳播

        loss = loss_fn(outputs, imp_label)	#計算損失
        optimizer.zero_grad()			#每一批次計算完成後梯度清零
        loss.backward()					#反向傳播
        optimizer.step()				#梯度更新


    #然後在每一批次訓練完成後,用交叉驗證集檢驗
    traces_val = Variable(valX)
    outputs_val = cnn(traces_val)
    imp_val = Variable(valImp)
    loss_val = loss_fn(outputs_val, imp_val)

    #列印
    print ('Epoch [%d/%d], Iter [%d], Train Loss: %.6f, Validation Loss: %.6f'
           %(epoch+1, num_epochs, i+1, loss.data.cpu().numpy(), loss_val.data.cpu().numpy()))

輸出如下:

Epoch [1/500], Iter [1400], Train Loss: 0.357334, Train Loss: 0.099698, Validation Loss: 0.489520, Validation Loss: 0.121744
Epoch [2/500], Iter [1400], Train Loss: 0.279390, Train Loss: 0.088157, Validation Loss: 0.453510, Validation Loss: 0.117180
Epoch [3/500], Iter [1400], Train Loss: 0.249917, Train Loss: 0.083377, Validation Loss: 0.435483, Validation Loss: 0.114828
Epoch [4/500], Iter [1400], Train Loss: 0.277595, Train Loss: 0.087873, Validation Loss: 0.431772, Validation Loss: 0.114337
...

模型的損失值曲線圖如下:

Pytorch實現波阻抗反演
### (3)預測

Step 1 預測測試集

通常測試集也是有標籤的,我們能夠直觀地對比模型的預測效果與精度。

#抽取出待測的地震道
sampleNumber = 25
TestingSetSeismicTrace = Variable(testX[sampleNumber:sampleNumber+1,:,:])	#輸入資料

CNN_ImpedancePrediction = cnn(TestingSetSeismicTrace)		#預測結果

圖中展示了4道的真實波阻抗與CNN預測波阻抗:

Pytorch實現波阻抗反演
**Step 2 儲存&載入模型**

訓練一次模型往往需要花費很長時間,將模型儲存下來,在需要的時候載入,避免關閉程式後再重新訓練。

#儲存模型
with open('cnn.pkl', 'wb') as pickle_file:
     torch.save(cnn.state_dict(), pickle_file)

#載入模型
#載入時需要先例項化物件
cnn_new = CNN()
with open('cnn.pkl', 'rb') as pickle_file:
     cnn_new.load_state_dict(torch.load(pickle_file))

Step 3 應用

應用時往往沒有標籤(波阻抗),需要模型由已知資料預測未知。

dataframe = sio.loadmat('HardTest_DataSyn_Ricker30.mat')

TestingSeismic = dataframe['wz_with_multiples']
TestingSeismic = TestingSeismic.transpose()[:, 0:333]

TestingImpedance = dataframe['IpTimeVec']/1e6
TestingImpedance = TestingImpedance.transpose()[:, 0:333]

print(TestingImpedance.shape[1])

sampleNumber = 25

#抽取應用資料的地震道
AppSeismicTrace = AppSeismic[sampleNumber,:]

#Numpy.array-->Tensor-->Variable
AppSeismicTraceTorch = torch.FloatTensor(np.reshape(AppSeismicTrace, (1,1,  AppSeismicTrace.shape[0])))	#輸入尺寸與訓練集保持一致,第一個1表示應用當前單個樣本

AppSeismicTraceTorch = Variable(AppSeismicTraceTorch)	#封裝成Variable
#用載入的cnn_new模型去預測
CNN_ImpedancePrediction = cnn_new(AppSeismicTraceTorch)
Pytorch實現波阻抗反演
於是,我們便完成了,通過建立神經網路模型,從地震記錄反演波阻抗的整個過程。

完整程式碼

import torch.nn as nn
import torch.autograd
from torch.autograd import Variable
import scipy.io as sio
import matplotlib.pyplot as plt
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import numpy as np
from random import randint

#******************匯入資料*********************
dataframe = sio.loadmat('Train_DataSyn_Ricker30.mat')

#從檔案中分別提取提取地震道與波阻抗資料
Seismic_data = dataframe['Seismic']				#地震道
Impedance_data = dataframe['Impedance']/1e6		#波阻抗

#隨機抽取部分道,作為訓練集
howMany = 2020
np.random.seed(9)	#隨機種子,便於復現
indxRand = [randint(0,dataframe['Seismic'].shape[1]-1) for p in range(0,howMany)]	#隨機索引

#地震道
Seismic_data = Seismic_data.transpose()	#轉置
Seismic_data = Seismic_data[indxRand,:]	#通過索引抽取
#波阻抗
Impedance_data = Impedance_data.transpose()
Impedance_data = Impedance_data[indxRand,:]

dt = 4.3875e-4
time = np.linspace(0,(Seismic_data.shape[1]-1)*dt,Seismic_data.shape[1])

#*********************分割資料集****************************

#其中驗證集500個,測試集1000個,剩餘(520)為訓練集
howManyToValidate = 500
howManyToTest = 1000

#對輸入Seismic與標籤Imp進行相同的索引與處理
#用numpy索引切片的方式進行劃分
valX = (Seismic_data[:howManyToValidate,:])
testX = (Seismic_data[howManyToValidate:howManyToValidate+howManyToTest,:])
trainX = (Seismic_data[howManyToValidate+howManyToTest:,:])

valImp = (Impedance_data[:howManyToValidate,:])
testImp = (Impedance_data[howManyToValidate:howManyToValidate+howManyToTest,:])
trainImp = (Impedance_data[howManyToValidate+howManyToTest:,:])

#轉為torch中的Tensor格式
#此時資料為(道數,取樣點數)的二維陣列,按照torch的輸入格式整理為(道數,資料高度,資料長度),便於後續輸入道CNN網路中
valX = torch.FloatTensor(np.reshape(valX, (valX.shape[0], 1, valX.shape[1])))
testX = torch.FloatTensor(np.reshape(testX, (testX.shape[0], 1, testX.shape[1])))
trainX = torch.FloatTensor(np.reshape(trainX, (trainX.shape[0], 1, trainX.shape[1])))

valImp = torch.FloatTensor(np.reshape(valImp, (valImp.shape[0], 1, valImp.shape[1])))
testImp = torch.FloatTensor(np.reshape(testImp, (testImp.shape[0], 1, testImp.shape[1])))
trainImp = torch.FloatTensor(np.reshape(trainImp, (trainImp.shape[0], 1, trainImp.shape[1])))


#********************建模*********************
noOfNeurons = 60		#定義卷積核個數
dilation = 1
kernel_size = 300		#卷積核尺寸
stride = 1				#卷積核滑動步長
padding = int(((dilation*(kernel_size-1)-1)/stride-1)/2)	#0填充個數

class CNN(nn.Module):
    def __init__(self):					#建構函式
        super(CNN, self).__init__()		#前面三行為固定格式
        self.layer1 = nn.Sequential(	#nn.Sequential為一個順序的容器
            #Conv1d的前兩個引數分別表示輸入通道數與輸出通道數,又有卷積核個數=卷積層輸出通道數
            nn.Conv1d(1, noOfNeurons, kernel_size=kernel_size, stride=1, padding = padding+1),#卷積層
            nn.ReLU())	#ReLu啟用函式

        self.layer2 = nn.Sequential(
            nn.Conv1d(noOfNeurons, 1, kernel_size=kernel_size, stride=1, padding = padding+2),
            nn.ReLU())

    def forward(self, x):	#在forward中將網路像搭積木一樣連線起來
        out = self.layer1(x)
        out = self.layer2(out)

        return out

cnn = CNN()		#例項物件

#*************************訓練************************
#定義與訓練有關的超引數
num_epochs = 500		#迭代輪數
batch_size = 1			#批次大小
learning_rate = 0.001	#學習率
batch_size_tot = trainX.shape[0]
no_of_batches = int((batch_size_tot - batch_size_tot%batch_size)/batch_size)	#總批數

loss_fn = nn.MSELoss()	#損失函式,此處為一個迴歸任務,採用均方根誤差作為損失值
optimizer = torch.optim.Adam(cnn.parameters(), lr=learning_rate)	#優化器選擇為Adam

for epoch in range(num_epochs):		#開始迭代
    for ii in range(no_of_batches):
     	#地震道資料
        #通過手動索引的方式建立batch
        #使用Variable對Tensor進行封裝,便於改變Tensor的.data、.grad、.grad_fn屬性
        traces = Variable(trainX[ii*batch_size:(ii+1)*batch_size,:,:])
        imp_label = Variable(trainImp[ii*batch_size:(ii+1)*batch_size,:,:])

        '''以下5行程式碼為固定格式,幾乎所有pytorch網路都是同樣的'''
        outputs = cnn(traces)	#將訓練資料輸入道cnn模型中,前向傳播

        loss = loss_fn(outputs, imp_label)	#計算損失
        optimizer.zero_grad()			#每一批次計算完成後梯度清零
        loss.backward()					#反向傳播
        optimizer.step()				#梯度更新


    #然後在每一批次訓練完成後,用交叉驗證集檢驗
    traces_val = Variable(valX)
    outputs_val = cnn(traces_val)
    imp_val = Variable(valImp)
    loss_val = loss_fn(outputs_val, imp_val)

    #列印
    print ('Epoch [%d/%d], Iter [%d], Train Loss: %.6f, Validation Loss: %.6f'
           %(epoch+1, num_epochs, ii+1, loss.data.cpu().numpy(), loss_val.data.cpu().numpy()))

#**********************測試集視覺化對比**************************
#取出4道
sampleNumbers = np.array([21,50,162,206])

fig, axs = plt.subplots(1, 4, sharey=True)
axs[0].invert_yaxis()
fig.suptitle('Samples of Testing Data Predictions')

for i in range(4):
    sampleNumber = sampleNumbers[i];

    TestingSetSeismicTrace = Variable(testX[sampleNumber:sampleNumber+1,:,:])

    CNN_ImpedancePrediction = cnn(TestingSetSeismicTrace)

    #Numpy資料作圖用
    TestingSetSeismicTrace = testX[sampleNumber,:].numpy().flatten()
    TestingSetImpedanceTrace = testImp[sampleNumber,:].numpy().flatten()
    CNN_ImpedancePrediction = CNN_ImpedancePrediction.data.cpu().numpy().flatten()

    line1, = axs[i].plot(TestingSetImpedanceTrace, time, 'r-')
    axs[i].set_xlabel('Impedance')
    if i==0:
        axs[i].set_ylabel('Time(s)')
    line2, = axs[i].plot(CNN_ImpedancePrediction, time, 'k--')
lgd = plt.legend((line1, line2), ('True Impedance', 'CNN Impedance'), bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)

#********************模型儲存&載入*********************

#儲存模型
with open('cnn.pkl', 'wb') as pickle_file:
     torch.save(cnn.state_dict(), pickle_file)

#載入模型
#載入時需要先例項化物件
cnn_new = CNN()
with open('cnn.pkl', 'rb') as pickle_file:
     cnn_new.load_state_dict(torch.load(pickle_file))

#**********************應用************************
Appdataframe = sio.loadmat('HardTest_DataSyn_Ricker30.mat')

AppSeismic = Appdataframe['Seismic']
AppSeismic = AppSeismic.transpose()

sampleNumber = 25

#抽取應用資料的地震道
AppSeismicTrace = AppSeismic[sampleNumber,:]

#Numpy.array-->Tensor-->Variable
AppSeismicTraceTorch = torch.FloatTensor(np.reshape(AppSeismicTrace, (1,1,  AppSeismicTrace.shape[0])))	#輸入尺寸與訓練集保持一致,第一個1表示應用當前單個樣本

AppSeismicTraceTorch = Variable(AppSeismicTraceTorch)	#封裝成Variable
#用載入的cnn_new模型去預測
CNN_ImpedancePrediction = cnn_new(AppSeismicTraceTorch)

#作圖
fig, ax1 = plt.subplots()
ax1.plot(time, AppSeismicTrace, 'b-')
ax1.set_xlabel('time (s)')

ax1.set_ylabel('Seismic Amplitude', color='b')
ax1.tick_params('y', colors='b')

ax2 = ax1.twinx()
ax2.plot(time, CNN_ImpedancePrediction.detach().numpy().flatten(), 'r-')
ax2.set_ylabel('Impedance', color='r')
ax2.tick_params('y', colors='r')

fig.tight_layout()
plt.show()

資料與原始碼CNN_based_impedance_inversion/Challenge_testing_Facies_variograms at master · vishaldas/CNN_based_impedance_inversion (github.com)

引用文獻

[1] Vishal Das;Ahinoam Pollack;Uri Wollner;Tapan Mukerji.Convolutional neural network for seismic impedance inversion[J].Geophysics,2019,Vol.84(6): R869-R880

[2]王治強. 稀疏脈衝反褶積及其波阻抗反演研究[D].中國石油大學(北京),2018.DOI:10.27643/d.cnki.gsybu.2018.001091.

[3] Activation and loss functions (part 1) · Deep Learning (atcold.github.io)

[4] https://www.w3cschool.cn/article/78828381.html

相關文章