地震波阻抗反演是在勘探與開發期間進行儲層預測的一項關鍵技術。地震波阻抗反演可消除子波影響,僅留下反射係數,再通過反射係數計算出能表徵地層物性變化的物理引數。常用的有道積分、廣義線性反演、稀疏脈衝反演、模擬退火反演等技術。
隨著勘探與開發的深入,研究的地質目標已經從大套厚層砂體轉向薄層砂體,而利用常規波阻抗反演方法刻畫薄層砂體不僅要消耗大量人力、物力,且反演得到的波阻抗精度也難以滿足實際需求。近年來,深度學習在地震反演和解釋等地震領域顯現出了巨大的潛力,其中,折積神經網路(Convolution Neural Networks,CNN)是解決地震反演的一大有力工具。本文將介紹如何搭建CNN網路進行波阻抗反演。
我們將使用Python語言以及Pytorch深度學習框架完成實驗,因此首先需要設定Python環境。
下載地址Anaconda | The World's Most Popular Data Science Platform
Anaconda是一個開源的Python和R語言的發行版本,用於計算科學(資料科學、機器學習、巨量資料處理和預測分析),致力於簡化軟體包管理系統和部署。
安裝過程「下一步、下一步」,在這一步時,勾選新增Anaconda到環境變數(注意:在環境變數中有QT、R語言等時,須在「編輯環境變數」中手動將Anaconda路徑下移,防止變數、路徑被覆蓋)
(tips:Python和Anaconda都有蟒蛇的意思)
選擇Anaconda主要有幾個原因:
隨Anaconda一同下載的有
Anaconda Navigator:是包含在Anaconda中的圖形化使用者介面,使用者可以通過Anaconda Navigator啟動應用,在不使用命令列的情況下管理軟體包、建立虛擬環境和管理路徑。
Jupyter Notebook:是一個基於Web的互動式計算環境,用於建立Jupyter Notebook檔案。這類檔案是一個JSON檔案,包含一個有序的輸入/輸出單元格列表,這些單元格可以包含程式碼、文字(使用Markdown語言)、數學、圖表和富媒體 (Rich media),通常以「.ipynb」結尾擴充套件。Jupyter Notebook的最大的優點就是可以「做一步看一步」,對於初學和開荒Python專案時,可以更簡單、高效地程式設計。
Spyder,是一個開源、免費的Python整合式開發環境(IDE)。其最大優點是模仿MATLAB的「工作空間」的功能,便於觀察變數的值、維度、型別等資訊。
Anaconda Prompt:等同Window PowerShell和cmd。
下載地址Download PyCharm: Python IDE for Professional Developers by JetBrains
由捷克JetBrains公司專為Python打造,PyCharm具備一般IDE的功能,比如,偵錯、語法高亮、Project管理、程式碼跳轉、智慧提示、自動完成、單元測試、版本控制等等。PyCharm是商業軟體,與Visual Stidio類似,有收費的專業版(Professional)和免費的社群版(Community),社群版能夠滿足一般程式設計需要,專業版整合了Jupyter Notebook並且支援遠端開發。
下載Community版後,安裝過程除了調整路徑外,沒有需要操作的地方,非必須的可選項都不打勾,然後「下一步下一步」。
安裝完成後來到新建工程,以及設定Python環境的環節:
第一次建立工程、載入環境需要等上幾分鐘,讀條完成後,這個工程的Python環境及編輯器就設定完成了。
首先我們需要在本地安裝Pytorch,開啟Anaconda Prompt或Window PowerShell或cmd,安裝指令為
pip install torch torchvision torchaudio
或者
conda install pytorch torchvision torchaudio cpuonly -c pytorch
即可安裝CPU版Pytorch(GPU版的安裝稍複雜,此處不作詳細介紹),隨後開始本次專案的實戰部分。
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,:]
資料展示如下:
我們的目標即是通過建立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])))
Step 1 建立網路
採用一維折積神經網路(1D CNN)進行波阻抗的生成與預測,具體步驟與結構如下圖。
建立神經網路通常也就是建立一個繼承自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]之間,保證資料穩定,波動幅度小
缺點:
如果我們在多個層中堆疊sigmoid,則系統學習效率可能低下,並且需要仔細初始化。因此,對於深度神經網路,首選ReLU函數。
ReLU(Rectified Linear Unit,修正線性單元)——nn.ReLU()
\(\text{ReLU}(x) = (x)^{+} = \max(0,x)\)
優點:
缺點:
常用的啟用函數還有許多:
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
...
模型的損失值曲線圖如下:
Step 1 預測測試集
通常測試集也是有標籤的,我們能夠直觀地對比模型的預測效果與精度。
#抽取出待測的地震道
sampleNumber = 25
TestingSetSeismicTrace = Variable(testX[sampleNumber:sampleNumber+1,:,:]) #輸入資料
CNN_ImpedancePrediction = cnn(TestingSetSeismicTrace) #預測結果
圖中展示了4道的真實波阻抗與CNN預測波阻抗:
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)
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()
參照文獻
[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