吳易智， 范宜仁*， 巫振觀， 鄧少貴，張盼， 陳詩宇， 尹中旭

1 中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院， 青島 266580 2 海洋國家實(shí)驗(yàn)室海洋礦產(chǎn)資源評價與探測技術(shù)功能實(shí)驗(yàn)室， 青島 266237

0 引言

憑借著較高的分辨率以及豐富的探測信息，陣列側(cè)向測井被廣泛用于大斜度井/水平井的油氣藏評價(Rasmus, 1982; Itskovich et al., 1998; 鄧少貴等, 2010).然而在致密性儲層電阻率評價過程中，由于地層裂縫和薄互層發(fā)育以及沉積顆粒的大小、磨圓不一等因素，導(dǎo)致儲層各向異性特征明顯，測井響應(yīng)復(fù)雜，為油氣層的識別與評價帶來一定的難度(Salazar and Torres-Verdín, 2009; 高杰等, 2010; Gao et al., 2013).因此，如何從測井響應(yīng)中準(zhǔn)確提取地層真實(shí)電阻率成為電測井領(lǐng)域亟待解決的難題之一.

非線性反演技術(shù)可以充分考慮井眼、侵入等因素對測井響應(yīng)的影響，成為大斜度井/水平井井旁電阻率剖面重構(gòu)的有效手段(Nam et al., 2010; Wang, 2011; 潘克家等, 2013).其算法包括確定性算法和隨機(jī)類算法.其中，確定性算法主要基于最小二乘理論，其核心在于通過不斷迭代模型參數(shù)的方式使得模擬數(shù)據(jù)與觀測數(shù)據(jù)的偏差最小(Osman and Laporte, 1996; Habashy and Abubakar, 2004; Madsen et al., 2004; 姚東華等, 2010; Hu and Fan, 2018)，該類算法包括最速下降法、高斯牛頓法以及列文伯格馬奎特法等(Anderson, 2001).由于能夠有效消除圍巖、層厚和侵入等因素的影響，確定性算法被廣泛用于2D陣列側(cè)向測井資料的快速處理(Wang et al., 2009; 王磊等, 2018).然而對于斜井各向異性等三維地層環(huán)境而言，正演效率過低，導(dǎo)致反演迭代過程中雅可比矩陣或海森矩陣求取耗時，同時該類算法對初值依賴性強(qiáng)，易陷于局部極小值，不利于陣列側(cè)向測井資料的快速精確處理(Pardo et al., 2007; Dubourg et al., 2017; Wu et al., 2020).不同于確定性算法，隨機(jī)類算法由單點(diǎn)或多點(diǎn)啟發(fā)，搜索步長和方向具有隨機(jī)性，可以在求解域內(nèi)全空間搜索，全局尋優(yōu)能力強(qiáng)，進(jìn)而能夠提高資料處理精度，此類算法包括粒子群優(yōu)化算法、模擬退火算法以及差分進(jìn)化算法等(Eiben et al., 1994;Yang and Torres-Verdín, 2011; Li et al., 2020).然而，兩類算法在尋優(yōu)過程中皆要不斷調(diào)用正演，尤其對于三維地層環(huán)境而言，正演速度仍是制約反演效率的關(guān)鍵因素.

近年來，隨著智能算法的興起，有研究學(xué)者將深度學(xué)習(xí)技術(shù)引入地球物理反演當(dāng)中(林年添等, 2018; Shahriari et al., 2020; Wu and Fan, 2021; Zhang et al.,2021)，以期能夠解決反演策略中正演效率低的難題.其中，反向傳播神經(jīng)網(wǎng)絡(luò)(BPNN)作為一種強(qiáng)大的監(jiān)督類學(xué)習(xí)架構(gòu)，自從提出以來被廣泛用于特征識別以及函數(shù)逼近領(lǐng)域(Rumelhart and McClelland, 1986; Xu et al., 2018).然而，面對復(fù)雜輸入輸出問題，BPNN存在超參數(shù)優(yōu)化困難、遷移性差、訓(xùn)練參數(shù)過多以及不易收斂等缺點(diǎn).為此，深層卷積神經(jīng)網(wǎng)絡(luò)(CNN)憑借著局部連接和權(quán)重共享的優(yōu)點(diǎn)逐漸代替BPNN(Fukushima, 1980; LeCun et al., 1998; Krizhevsky et al., 2017).在測井領(lǐng)域，Zhu等(2018)通過建立五層卷積神經(jīng)網(wǎng)絡(luò)實(shí)現(xiàn)了地層巖性特征的快速識別；Li等(2019)基于1D模型和多層卷積神經(jīng)網(wǎng)絡(luò)實(shí)現(xiàn)了大斜度井/水平井隨鉆電磁波測井的快速反演；Jin等(2019)通過優(yōu)化卷積神經(jīng)網(wǎng)絡(luò)的損失函數(shù)實(shí)現(xiàn)了隨鉆測井領(lǐng)域的實(shí)時地質(zhì)導(dǎo)向.然而，上述研究皆是針對簡單輸入輸出問題，而對于陣列側(cè)向測井3D反演問題，國內(nèi)外尚無研究.因此，針對陣列側(cè)向測井資料處理，如何提升高維正演效率，如何保證反演精度是當(dāng)前智能測井領(lǐng)域亟需解決的難題.

為此，本文從三維有限元技術(shù)(3D-FEM)出發(fā)，首先對斜井各向異性地層參數(shù)進(jìn)行敏感性分析，確定侵入、各向異性和地層傾角對陣列側(cè)向測井的敏感性大?。黄浯?，建立三維陣列側(cè)向測井響應(yīng)數(shù)據(jù)庫，引入二維卷積神經(jīng)網(wǎng)絡(luò)(2D-CNN)技術(shù)，并基于二進(jìn)制轉(zhuǎn)換將地層模型轉(zhuǎn)為網(wǎng)絡(luò)輸入層所需的2D圖像，提高網(wǎng)絡(luò)識別局部特征的能力；接著，建立深層2D-CNN結(jié)構(gòu)，并引入“丟棄”(Dropout)和“填充”(Padding)機(jī)制，實(shí)現(xiàn)三維陣列側(cè)向快速正演；最后，基于敏感性分析和快速正演，將多種群遺傳算法與列文伯格馬奎特算法相結(jié)合(MPGA-LM)，對傾斜各向異性地層的電阻率剖面進(jìn)行了快速精確重構(gòu).

1 陣列側(cè)向測井?dāng)?shù)值模擬

1.1 三維有限元數(shù)值模擬

本文以斯倫貝謝研發(fā)的高分辨率陣列側(cè)向測井儀HRLA為例，儀器包括一個主電極A0、六對屏蔽電極A1-A6(A1′-A6′)及兩對監(jiān)督電極M1、M2(M1′、M2′)，屏蔽電極和監(jiān)督電極關(guān)于主電極對稱，通過改變不同電極間的收發(fā)組合方式，可以形成6種具有不同探測深度的視電阻率曲線，探測深度由淺及深分別命名為RLA0，RLA1，RLA2，RLA3，RLA4，RLA5.

不同探測模式視電阻率可以通過測量監(jiān)督電極電位進(jìn)行計(jì)算，其中淺探測模式RLA0需要測量同側(cè)兩個監(jiān)督電極的電勢差值ΔUM1M2(RLA0)，其他探測模式只需測量任意監(jiān)督電極上電勢，其視電阻率計(jì)算公式如下：

(1)

(2)

式中，i為1～5，分別代表上述5種工作方式；Ki為第i種測井模式的電極系數(shù)，可在均勻地層條件下獲得；RRLAi為第i種測井模式的視電阻率；UM1(RLAi)為第i種測井模式的主電極電位；I0(RLAi)為第i種測井模式的主電極電流.

可以看出，求解陣列側(cè)向測井視電阻率的關(guān)鍵在于得到監(jiān)督電極上的電勢，為此必須將整個空間的電位分布函數(shù)求出.由于陣列側(cè)向測井用的是低頻交流電，故可視為穩(wěn)定電流，進(jìn)而該問題轉(zhuǎn)為穩(wěn)流場電勢分布函數(shù)的求解問題.

根據(jù)邊界條件和變分原理，陣列側(cè)向測井響應(yīng)可歸結(jié)為求取泛函極值問題，泛函構(gòu)造式如下：

(3)

式中，Λ為除去電極部分的求解區(qū)域，Ii和Ui分別代表儀器每個電極的電流和電勢大小.

在電阻率不同的子空間交界面上，電流、電位具有連續(xù)性，滿足如下邊界條件：

(4)

式中，“-”表示交界面左側(cè)(上側(cè))，“+”表示交界面右側(cè)(下側(cè))；n為交界面單位法向量.

在無窮遠(yuǎn)邊界和絕緣電極表面還需滿足Dirichlet和Neumann邊界條件，如下所示：

U|Γ1=0,

(5)

式中，Γ1和Γ2分別為無窮遠(yuǎn)邊界和絕緣電極表面.

在主電極和屏蔽電極表面，滿足等值面邊界條件：

(6)

式中，UAi為電極表面電勢，Γ3主電極或屏蔽電極表面，n為電極表面單位法向量，Ii為電極發(fā)射電流，i為0～5.

為求解上述泛函的最小值，一般采用三維有限元方法(3D-FEM)，最后采用前線解法求取含有大型稀疏矩陣的線性方程組，如式(7)，該方法由于可以在求解過程中對剛度矩陣邊消元邊安裝，可有效提升正演計(jì)算速度(張庚驥, 2009).

KΦ=F,

(7)

1.2 影響因素分析

陣列側(cè)向測井受泥漿侵入、地層傾角以及各向異性等多種因素影響，因此，確定各因素對測井響應(yīng)影響強(qiáng)弱對于反演初值選取、約束條件施加具有重要意義.本節(jié)系統(tǒng)分析了侵入深度、各向異性和地層相對傾角等因素對陣列側(cè)向測井響應(yīng)的敏感性，定義敏感因子S如下：

(8)

其中，i為陣列側(cè)向探測模式，此處取1～5，x為地層參數(shù)，分別代表侵入深度Di，各向異性系數(shù)λ和地層傾角θ.

首先，分析侵入對陣列側(cè)向測井響應(yīng)的影響.以無限厚一層模型為例，地層為各向同性，侵入帶電阻率為2 Ωm，侵入深度Di變化范圍為0～3 m，侵入帶與原狀地層電阻率對比度范圍Rat為1～20.敏感性分布圖如圖1所示，橫縱坐標(biāo)分別為侵入深度和電阻率對比度.可以看出，陣列側(cè)向測井對侵入深度敏感因子最大約為2，且隨著探測深度的增大，對侵入深度的敏感區(qū)間也隨之增大：淺、中、深探測模式RLA1、RLA3和RLA5的敏感區(qū)域分別為0.2～0.5 m、0.3～1.0 m和0.6～2.3 m.

圖1 HRLA對侵入深度的敏感性 (a)?RLA1/?Di； (b)?RLA3/?Di； (c)?RLA5/?DiFig.1 The sensitivity of HRLA to invasion depth

接著，分析各向異性和地層傾角對陣列側(cè)向測井的影響.以無限厚模型為例，地層為各向異性，地層水平電阻率為2 Ωm，各向異性系數(shù)變化范圍為1～4，地層相對傾角變化范圍為0°～90°.圖2和圖3分別表示各向異性和地層相對傾角的敏感性分布圖，橫縱坐標(biāo)分別為各向異性系數(shù)和地層相對傾角.可以看出：(1)陣列側(cè)向測井對各向異性系數(shù)和地層相對傾角的敏感因子最大約為1；(2)陣列側(cè)向測井在高角度地層(60°～90°)對各向異性最為敏感，且探測深度較深的測井模式對各向異性的敏感性較大，原因在于當(dāng)?shù)貙酉鄬A角逐漸增大時，儀器電極發(fā)射電流方向逐漸朝著垂直電阻率方向傾斜.因此，垂直電阻率對測井響應(yīng)貢獻(xiàn)變大，各向異性敏感性增強(qiáng)，而深探測曲線因探測范圍最廣，受其影響最大；(3)陣列側(cè)向測井在中高角度(50°～80°)對地層相對傾角最為敏感，且不同測井模式對相對傾角的敏感性與各向異性相似，皆是探測深度較深的測井模式對其敏感性較大，原因在于地層傾角對陣列側(cè)向曲線的影響本質(zhì)是因?yàn)榈貙痈飨虍愋缘挠绊懀驗(yàn)榈貙觾A角增大，導(dǎo)致垂直電阻率的影響增大，進(jìn)而使得陣列側(cè)向五條曲線出現(xiàn)變化.

圖2 HRLA對各向異性的敏感性 (a)?RLA1/? λ； (b)?RLA3/? λ； (c)?RLA5/? λ.Fig.2 The sensitivity of HRLA to anisotropy

圖3 HRLA對地層相對傾角的敏感性 (a)?RLA1/? θ； (b)?RLA3/? θ； (c) ?RLA5/? θ.Fig.3 The sensitivity of HRLA to dipping angle

2 深層卷積神經(jīng)網(wǎng)絡(luò)快速正演模型

2.1 深層卷積神經(jīng)網(wǎng)絡(luò)理論

深層卷積神經(jīng)網(wǎng)絡(luò)憑借著強(qiáng)大的泛化特性以及學(xué)習(xí)能力被廣泛用于圖像分類、語音檢測以及物體識別等復(fù)雜特征識別領(lǐng)域.一般而言，卷積神經(jīng)網(wǎng)絡(luò)包括輸入層、卷積層、匯聚層(池化層)、全連接層和輸出層，本文以經(jīng)典二維卷積神經(jīng)網(wǎng)絡(luò)——LeNet-5為例，其結(jié)構(gòu)如圖4所示.其中，輸入層為三通道RGB圖像或者單通道灰度圖像，卷積層包含多個卷積核，不同的卷積核相當(dāng)于不同的特征提取器，結(jié)構(gòu)為一個二維矩陣，圖像經(jīng)過卷積層可提取多個局部特征并保存；匯聚層的作用在于特征選擇，降低特征數(shù)量，從而減少參數(shù)數(shù)量，其核心思想在于步長間隔的選擇，達(dá)到“降維”的作用；全連接層與全連接前饋神經(jīng)網(wǎng)絡(luò)相同，不再贅述.

圖4 LeNet-5網(wǎng)絡(luò)結(jié)構(gòu)Fig.4 The architecture of LeNet-5

相較于全連接神經(jīng)網(wǎng)絡(luò)，卷積神經(jīng)網(wǎng)絡(luò)最大的特點(diǎn)在于局部連接和權(quán)重共享，如圖5所示，圖5a為傳統(tǒng)的二層全連接神經(jīng)網(wǎng)絡(luò)，包括3個輸入，5個輸出，需訓(xùn)練參數(shù)包括15個權(quán)重和1個偏置共16個參數(shù)；圖5b為二層卷積神經(jīng)網(wǎng)絡(luò)，包括3個輸入，5個輸出，由于輸入層每個神經(jīng)元只需要與下層局部窗口內(nèi)的神經(jīng)元連接，且每個窗口內(nèi)權(quán)重共享，若卷積核大小為3，則需訓(xùn)練參數(shù)包括3個權(quán)重和1個偏置共4個參數(shù)，可大大減少訓(xùn)練參數(shù)的數(shù)量，提高訓(xùn)練效率.

圖5 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu) (a) 全連接層； (b) 卷積層.Fig.5 The connection feature of network (a) Fully connected layer； (b) Convolution layer.

同時，為提高深層神經(jīng)網(wǎng)絡(luò)的訓(xùn)練效率，通常卷積層和全連接層分別引入填充(Padding)和丟棄(Dropout)機(jī)制，如圖6所示，通過在圖像卷積之前填充一層像素，使卷積之后的圖片與原圖尺寸相同，進(jìn)而使得卷積層得以加深，提高訓(xùn)練精度，如圖6a；另外通過隨機(jī)丟棄部分神經(jīng)元，可有效防止網(wǎng)絡(luò)過擬合，加快訓(xùn)練進(jìn)程，如圖6b.

圖6 神經(jīng)元連接機(jī)制 (a) “Padding”機(jī)制； (b) “Dropout”機(jī)制.Fig.6 Neuronal connection mechanism (a) “Padding” mechanism； (b) “Dropout” mechanism.

2.2 數(shù)據(jù)庫構(gòu)建及模型可視化

神經(jīng)網(wǎng)絡(luò)模型的可靠性和穩(wěn)定性很大程度上依賴于訓(xùn)練數(shù)據(jù)的精度，因此，本文基于嚴(yán)格三維有限元算法，構(gòu)建大型數(shù)據(jù)集，并將其載入深層2D-CNN進(jìn)行訓(xùn)練，以此建立快速回歸模型.考慮圍巖和層厚對陣列側(cè)向測井響應(yīng)的影響，數(shù)據(jù)集基于三層地層模型構(gòu)建，如圖7a所示，該模型亦作為2D-CNN的輸入.具體地，令兩側(cè)圍巖為各向同性地層，電阻率相同，中間層為有侵各向異性地層，模型參數(shù)范圍如下：各向異性系數(shù)λ為1～4，分為13組，間隔0.25；侵入深度Di為0～1 m，分為11組，間隔0.1 m；圍巖電阻率Rs、沖洗帶水平電阻率Rxoh和目的層水平電阻率Rth范圍皆為0.1～500 Ωm，將其以對數(shù)形式各分為20組，對數(shù)間隔為0.19；兩側(cè)圍巖為半無限厚地層；目的層厚度H為0～5 m，分為16組，間隔0.33 m；測量點(diǎn)為以目的層中點(diǎn)為中心，上下各取兩個點(diǎn)，共5個測量點(diǎn)，采樣點(diǎn)間隔為0.1 m.需要注意的是，為簡化數(shù)據(jù)庫構(gòu)建流程，將其按地層傾角分為10個子庫，每個子庫地層傾角間隔10°.因此，每個子庫有13×11×203×16，共18,304,000個樣本，每個樣本包括六個模型參數(shù)，共109,824,000個模型數(shù)據(jù)，將其儲存至變量M中，大小為18304000×6；同時每個樣本含有五個測量點(diǎn)，每個測量點(diǎn)包括五個電阻率值(RLA1～RLA5)，因此，每組共包含457,600,000個數(shù)據(jù)點(diǎn)，將數(shù)據(jù)點(diǎn)作為輸出儲存至變量D中，大小為18304000×25.在計(jì)算機(jī)內(nèi)存中，每個數(shù)據(jù)點(diǎn)占8個字節(jié)，因此，每組數(shù)據(jù)大小約為0.944 GB，數(shù)據(jù)庫總大小為9.44 GB.

對于復(fù)雜模型回歸問題，1D-CNN訓(xùn)練精度無法滿足需求，因此，本文采用學(xué)習(xí)和泛化能力較強(qiáng)的2D-CNN結(jié)構(gòu).而對于2D-CNN而言，輸入層普遍為2D圖像，為滿足網(wǎng)絡(luò)需求，本文基于二進(jìn)制轉(zhuǎn)換將一系列一維模型數(shù)組轉(zhuǎn)化為二維字符串，進(jìn)而合成一系列二值圖像(Zhong et al., 2019).具體轉(zhuǎn)換過程如圖7所示，其中圖7a為三層地層模型，具體模型參數(shù)如下：兩側(cè)圍巖為半無限厚各向同性地層，電阻率為2.564 Ωm，中間層層厚為3.156 m，侵入深度為0.428 m，侵入帶和目的層水平電阻率分別為11.278 Ωm和47.433 Ωm，各向異性系數(shù)為1.3，圖7b為模型參數(shù)對應(yīng)的二進(jìn)制矩陣，圖7c為基于二進(jìn)制矩陣轉(zhuǎn)換為的二值圖像，每個模型可唯一對應(yīng)一個二值圖像.由于數(shù)據(jù)庫中最大數(shù)值為8000，且小數(shù)后保留3位，故確定二進(jìn)制字符串長度為32，若字符串長度過長，可自動補(bǔ)零.

圖7 模型可視化過程 (a) 地層模型； (b) 模型參數(shù)對應(yīng)二進(jìn)制矩陣； (c) 二值圖像.Fig.7 The process of model visualization (a) Formation model； (b) The binary matrix correspond to model parameters； (c) Binary image corresponding to formation model.

2.3 2D-CNN損失函數(shù)及數(shù)據(jù)歸一化

由于前饋神經(jīng)網(wǎng)絡(luò)的局限性導(dǎo)致訓(xùn)練精度偏低，為此本文將反向傳播技術(shù)應(yīng)用于2D-CNN網(wǎng)絡(luò)架構(gòu)中，定義損失函數(shù)如下：

(9)

g(w,b,x)=w?x+b,

(10)

σ(z)=max(z,0),

(11)

(12)

式中，L為神經(jīng)網(wǎng)絡(luò)層數(shù)，Kl和Ml分別為第l層卷積核數(shù)量和卷積核尺寸.

接著，通過計(jì)算損失函數(shù)關(guān)于每層參數(shù)的梯度進(jìn)行反向傳播，為不失一般性，本文使用鏈?zhǔn)角髮?dǎo)法則對卷積層第l層的參數(shù)進(jìn)行求導(dǎo):

(13)

由式(10)和式(11)，易得

(14)

(15)

式中，k為第l層第k個特征映射，U為第l+1層的卷積核數(shù)量，X(l,k)=σ(g(l,k)).

同時，由于每組樣本中有5個測量點(diǎn)，即Y= (Y1，Y2，Y3，Y4，Y5)T，每個測量點(diǎn)包含五個視電阻率值，即Yi=(RLAi1,RLAi2,RLAi3,RLAi4,RLAi5)T.為提升訓(xùn)練精度，可將數(shù)據(jù)進(jìn)行歸一化處理，本文采用最小最大值歸一化方法，將輸出數(shù)據(jù)統(tǒng)一歸一化至[0,1]，公式如下：

(16)

需要注意的是，對于卷積神經(jīng)網(wǎng)絡(luò)模型預(yù)測的數(shù)據(jù)而言，需將其反歸一化至正常值.

2.4 深層卷積神經(jīng)網(wǎng)絡(luò)設(shè)計(jì)及測試

本文以9層卷積神經(jīng)網(wǎng)絡(luò)為例，對大型數(shù)據(jù)集進(jìn)行訓(xùn)練測試，網(wǎng)絡(luò)結(jié)構(gòu)及超參數(shù)定義如表1所示.網(wǎng)絡(luò)輸入為侵入深度、各向異性系數(shù)、侵入帶和原狀地層水平電阻率、圍巖層電阻率和目的層厚度所轉(zhuǎn)換為的二值圖像，即Input=BImage6×32×1, 輸出為陣列側(cè)向測井歸一化之后的視電阻率，即Output=RLA1×25.網(wǎng)格架構(gòu)搭建平臺為：Pycharm2020.2.3，Python3.7和Tensorflow2.0(GPU版本).為防止訓(xùn)練模型過擬合以及提升訓(xùn)練精度，引入Dropout和Padding機(jī)制，其中Dropout率為0.5，Padding選擇‘same’模式，即每層輸入輸出維度一致，進(jìn)而使得在每層輸入不丟失信息的前提下將卷積層延續(xù)，達(dá)到深化卷積層、提升訓(xùn)練精度的目的.除此之外，由于匯聚層將卷積層提取的多個特征融合，而對于本文回歸問題，每個模型特征都是唯一的，因此，為保持模型特征完整性，不使用匯聚層.同時反向傳播算法采用自適應(yīng)矩估計(jì)算法(Adam).訓(xùn)練過程中，隨機(jī)選取95%的數(shù)據(jù)集為訓(xùn)練集，剩余5%的數(shù)據(jù)為測試集.

表1 2D-CNN超參數(shù)定義Table 1 The definition of hyperparameters in 2D-CNN

為測試2D-CNN網(wǎng)格訓(xùn)練精度，本文將其訓(xùn)練結(jié)果與BPNN對比分析，其中BPNN網(wǎng)絡(luò)架構(gòu)與其類似，具體為：九層神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，激活函數(shù)為ReLU，隱含層神經(jīng)元個數(shù)分別為32，64，128，256，516，516，516，Dropout率為0.5，損失函數(shù)與反向傳播訓(xùn)練算法與2D-CNN一致.圖8為其訓(xùn)練精度與損失函數(shù)下降量的對比結(jié)果，綠色和黑色散點(diǎn)線分別代表2D-CNN和BPNN訓(xùn)練結(jié)果，可以看出兩者收斂速度大致相同，但訓(xùn)練精度卻相差較大，2D-CNN在100次迭代后訓(xùn)練精度即可達(dá)到99%左右，而BPNN為85%.除此之外，BPNN容易陷入局部極小值，精度無法隨著迭代的進(jìn)行而提升，相較而言，2D-CNN不僅可以快速收斂，而且訓(xùn)練精度隨著迭代的進(jìn)行逐漸接近100%.測試集與之規(guī)律類似，不再贅述.另外，通過計(jì)算2D-CNN和BPNN在測試集上的泛化誤差，計(jì)算公式如式(17).我們發(fā)現(xiàn)，兩者泛化誤差分別為1.28%和16.37%，從而也佐證了2D-CNN算法的高效性.

圖8 神經(jīng)網(wǎng)絡(luò)訓(xùn)練結(jié)果 (a) BPNN和2D-CNN預(yù)測精度； (b) BPNN和2D-CNN損失函數(shù).Fig.8 Neural network training results (a) Prediction accuracy of BPNN and 2D-CNN； (b) Loss function of BPNN and 2D-CNN.

(17)

式中，ntest為測試集樣本數(shù)，此處為211200；RLAij_2D-CNN為2D-CNN計(jì)算結(jié)果，RLAij_3D-FEM為3D-FEM計(jì)算結(jié)果.

除此之外，為驗(yàn)證訓(xùn)練模型的精度，我們隨機(jī)選取20個樣本點(diǎn)，分別測試中間層原狀地層水平電阻率與各向異性系數(shù)對模型預(yù)測精度的影響，如圖9所示，實(shí)線代表3D-FEM計(jì)算結(jié)果，散點(diǎn)表示2D-CNN預(yù)測結(jié)果.樣本點(diǎn)模型皆為三層模型：Rs為2 Ωm，H為4 m，Di為0.5 m，Rxoh為20 Ωm.其中圖9aλ為1.5，Rth變化范圍為1～500 Ωm；圖9b中Rth為100 Ωm，λ變化范圍為1～3.從對比結(jié)果可以看出，兩者吻合度極高，最大相對誤差分別為0.81%和0.24%.因此，基于深層卷積網(wǎng)絡(luò)的快速正演模型滿足精度需求，為斜井各向異性地層的快速反演提供有力支撐.

圖9 2D-CNN精度測試 (a) Rth對預(yù)測精度的影響； (b) λ對預(yù)測精度的影響.Fig.9 Prediction accuracy testing of 2D-CNN (a) The effect of Rth on the prediction accuracy； (b) The effect of λ on the prediction accuracy.

3 混合MPGA-LM算法及快速反演流程

3.1 代價函數(shù)定義

大斜度井/水平井陣列側(cè)向測井是典型的非線性反演問題，確定性算法L-M算法由于結(jié)合高斯牛頓算法和最速下降法的特性，能夠彌補(bǔ)全局尋優(yōu)能力不強(qiáng)或收斂性差的缺陷，反演過程中采用的代價函數(shù)如下：

(18)

式中，x為待反演參數(shù)矢量，x=(Di,Rxoh,Rth,λ,Rs)T，上標(biāo)T表示矩陣的轉(zhuǎn)置，S(x)為陣列側(cè)向儀器正演響應(yīng)，dobs為實(shí)際測量數(shù)據(jù)，‖·‖2表征L2范數(shù)，xp為模型已知參考矢量，通常是上一步迭代參數(shù)矢量，ξ為正則化參數(shù)，最后一項(xiàng)的作用是壓制測量噪聲，同時減小非線性反演出現(xiàn)的病態(tài)問題.反演的目的是得到一個x*使得代價函數(shù)最小，即x*=argmin{C(x)}.

進(jìn)一步，經(jīng)二階泰勒近似，可得第k步迭代步長xk為

xk=-JT(xk){S(xk)-dobs)}[JT(xk)J(xk)+ξkI]-1,

(19)

式中，I為單位矩陣，J為S在xk處的雅克比矩陣，ξk為第k次迭代的正則化參數(shù)大小.

為克服傳統(tǒng)差分方法計(jì)算雅可比矩陣耗時的缺陷，本文采用Broyden近似方法進(jìn)行雅可比矩陣的迭代更新(Broyden, 2000).

·(xk)T,

(20)

式中，J(xk+1)為第k+1步的雅可比矩陣.

3.2 基于MPGA的初值優(yōu)選策略

由于L-M算法對初值依賴性較強(qiáng)，若初值選取不適，算法易困于局部極小值，反演精度低.為此，本文基于卷積神經(jīng)網(wǎng)絡(luò)快速正演模型，同時為克服標(biāo)準(zhǔn)遺傳算法(GA)早熟早收斂的缺陷，利用多種群遺傳算法(MPGA)優(yōu)化初值選取，進(jìn)而提升L-M算法全局尋優(yōu)能力.

如圖10所示，多種群遺傳算法的核心在于將多個種群進(jìn)行初始化，即“二進(jìn)制”編碼，并同時進(jìn)行選擇、交叉和變異操作.種群間通過移民算子進(jìn)行協(xié)同進(jìn)化，最后通過適應(yīng)度函數(shù)分析，借助人工選擇算子選取種群最優(yōu)個體組成精華種群.本文定義適應(yīng)度函數(shù)如下：

圖10 多種群遺傳算法流程Fig.10 The framework of multi-population genetic algorithm

(21)

式中，p= 1,2,…,P，P為種群數(shù)量，n=1,2,…,N，N為每個種群個體數(shù)量，m為每個測量點(diǎn)數(shù)據(jù)量，此處為25.

除此之外，本文所采用的多種群遺傳算法超參數(shù)選擇如下：種群數(shù)目為5，個體數(shù)目為30，最大遺傳代數(shù)為50，變量的二進(jìn)制位數(shù)為32，交叉概率在[0.7，0.9]內(nèi)隨機(jī)產(chǎn)生，變異概率在[0.001，0.05]內(nèi)隨機(jī)產(chǎn)生.

為測試初值選取精度，以30°三層地層模型為例，兩側(cè)為各向同性圍巖，中間層為有侵各向異性地層.表2為測試模型參數(shù)分布以及多種群遺傳算法初值選取結(jié)果，Mode l1-Mode l5層厚分別為1～5 m，間隔為1 m.可以看出，受層厚影響，當(dāng)層厚小于3 m時，圍巖電阻率可準(zhǔn)確選取，平均相對誤差為4.4%；當(dāng)層厚大于3 m時，陣列側(cè)向測井響應(yīng)基本不反映圍巖信息，選取誤差急劇增大，而此時其他參數(shù)選取較為精確，整體而言，侵入深度平均相對誤差為18.6%，各向異性系數(shù)為17.8%，侵入帶和原狀地層水平電阻率分別為18.7%和6.8%，選取精度滿足L-M算法需求.

表2 初值選取結(jié)果Table 2 The selection result of initial value

3.3 混合MPGA-LM快速反演流程及有效性測試

實(shí)際資料處理過程中，針對帶有深度點(diǎn)的大型測井響應(yīng)數(shù)據(jù)集，我們將其由上到下分為若干子區(qū)域，每個子區(qū)域?yàn)槿龑幽Ｐ?，中間層為有限厚，地層厚度由層界面位置決定，而層界面位置可通過視電阻率曲線拐點(diǎn)或伽馬、自然電位等巖性曲線確定，上下地層視為無限厚地層進(jìn)行處理.圖11為混合MPGA-LM反演流程圖，由上節(jié)可知，經(jīng)MPGA優(yōu)選后的參數(shù)分布于反演代價函數(shù)真實(shí)解附近，故可將其作為L-M算法的初值，進(jìn)而提升反演效率；另一方面，由于當(dāng)層厚大于3 m時候，圍巖對陣列側(cè)向響應(yīng)近乎無影響，此時反演參數(shù)可降為4個，即Di、Rxoh、Rth和λ.同時，為提高訓(xùn)練精度，反演可循環(huán)N次(本文N為2)，即處理完整個井段資料后，可進(jìn)行二次處理.

圖11 混合MPGA-LM反演流程Fig.11 The workflow of hybrid MPGA-LM scheme

另外，為測試混合MPGA-LM快速反演算法的有效性，本文以三層模型為例，分別采用簡單初始模型與MPGA提供的初始模型進(jìn)行反演對比，反演算法皆為LM算法.地層傾角選擇30°，模型參數(shù)分布如表3所示.其中，在簡單初始模型中分別將深、淺測井響應(yīng)RLA1和RLA5作為Rxoh和Rth的初始值，Di和λ的初始值分別為0.5 m和1.5，表3為兩者方法反演結(jié)果及相對誤差，可以看出，較之簡單初始模型的LM算法，MPGA-LM可以更好的將地層剖面重構(gòu)，反演精度更高.以中間地層反演結(jié)果為例，基于簡單初始模型的LM算法反演后Di、λ、Rxoh和Rth的相對誤差(RE)分別為36.67%、33%、84.10%和23.08%；MPGA-LM反演后Di、λ、Rxoh和Rth的相對誤差分別為3.33%、0.50%、6.40%和2.86%，精度可以提高數(shù)倍至數(shù)十倍，由此本文提出的混合MPGA-LM快速反演算法的有效性得到了驗(yàn)證，為下文數(shù)值模擬實(shí)例的分析提供了保障.

表3 LM與MPGA-LM算法反演結(jié)果對比Table 3 The inverted results comparison between LM and MPGA-LM algorithm

4 數(shù)值模擬實(shí)例

本節(jié)通過建立21層俄克拉荷馬模型，進(jìn)一步驗(yàn)證上述反演方案的有效性.模型參數(shù)如下：井徑為8 inch，泥漿電阻率為0.1 Ωm，地層相對傾角分別

為0°、30°和60°，侵入深度范圍0～1.0 m，侵入帶水平電阻率分布范圍為0.2～30 Ωm，原狀地層水平電阻率范圍為0.2～100 Ωm，各向異性系數(shù)范圍為1～2.5 Ωm.圖12為該模型參數(shù)分布示意圖，縱坐標(biāo)為真實(shí)垂直深度(TVD).圖13為地層相對傾角分別為0°、30°和60°時陣列側(cè)向測井響應(yīng).本文將地層界面位置作為已知條件，僅反演侵入深度、侵入帶電阻率以及原狀地層電阻率.

圖12 俄克拉荷馬模型參數(shù)分布 (a) 侵入深度； (b) 侵入帶和原狀地層水平電阻率； (c) 各向異性系數(shù).Fig.12 The Oklahoma model parameters (a) Di； (b) Rxoh and Rth；(c) λ.

圖13 陣列側(cè)向測井響應(yīng)結(jié)果 (a) 0°； (b) 30°； (c) 60°.Fig.13 The responses of array laterolog

從測井響應(yīng)結(jié)果可以看出，泥漿濾液侵入對其影響嚴(yán)重，如模型第八層，該層為各向同性地層，厚度為5.2 m，侵入深度為0.3 m，地層真實(shí)電阻率為70 Ωm，而視電阻率值在0°、30°和60°時與地層電阻率的最大偏離度，即max(RLA/Rth)，分別為32.16%、31.81%和31.03%，遠(yuǎn)遠(yuǎn)偏離地層真實(shí)電阻率；除此之外，由于受地層各向異性的影響，導(dǎo)致響應(yīng)結(jié)果介于水平和垂直電阻率之間，如模型第16層，該層無侵入，地層水平和垂直電阻率分別為5 Ωm和20 Ωm，測井響應(yīng)在0°、30°和60°時的視電阻率大小分別在12 Ωm、13 Ωm 和15 Ωm附近，可以看出，60°地層對各向異性更為敏感，由此驗(yàn)證了2.3部分的結(jié)論.由上述可知，單單從測井曲線無法提取地層真實(shí)電阻率，因此需借助反演手段對地層真實(shí)電阻率剖面進(jìn)行重構(gòu).本文基于運(yùn)行內(nèi)存為8 GB 的英特爾i7-9700 處理器對60 m測井?dāng)?shù)據(jù)進(jìn)行反演處理，數(shù)據(jù)采樣間隔為0.1 m，循環(huán)反演次數(shù)為2次，反演總耗時約6.8 min.圖14是基于本文提出的反演策略對測井資料處理后的結(jié)果.

從圖14反演結(jié)果及圖15反演誤差曲線可以看出，就反演精度而言，地層水平電阻率>侵入深度>各向異性系數(shù)，三者的相對誤差范圍分別為0～8.38%、0.11%～6.98%和0～14.35%，該現(xiàn)象可由1.3節(jié)結(jié)論解釋，即測井響應(yīng)對侵入深度的敏感性大于地層各向異性；另一方面，地層傾角對侵入深度和地層水平電阻率的反演精度影響較小，而對地層各向異性影響較大，可以看出，三個角度下各向異性的反演相對誤差最大值分別為14.35%、7.79%和3.85%，該現(xiàn)象亦可由1.3節(jié)結(jié)論解釋，即高角度條件下陣列側(cè)向測井響應(yīng)對各向異性的敏感性較強(qiáng).

圖14 混合MPGA-LM反演結(jié)果 (a) 侵入深度； (b) 侵入帶和原狀地層水平電阻率； (c) 各向異性系數(shù).Fig.14 The inverted results using hybrid MPGA-LM algorithm (a) Di； (b) Rxoh and Rth； (c) λ.

圖15 反演結(jié)果相對誤差 (a) 侵入深度反演誤差； (b) 侵入帶和原狀地層水平電阻率反演誤差； (c) 各向異性系數(shù)反演誤差.Fig.15 The relative error of inverted results (a) The relative error of inverted Di； (b) The relative error of inverted Rxoh and Rth； (c) The relative error of inverted λ.

另外，由于在反演過程中調(diào)用的正演模型是2D-CNN快速計(jì)算模型，而非3D-FEM算法，而2D-CNN模型缺乏對各向異性圍巖層的考慮.因此，為進(jìn)一步分析反演誤差產(chǎn)生原因，我們將2D-CNN與3D-CNN的計(jì)算結(jié)果進(jìn)行誤差分析，平均相對誤差(MRE)計(jì)算公式如式(22)所示，具體結(jié)果如表4所示.

表4 2D-CNN與3D-FEM計(jì)算結(jié)果平均相對誤差Table 4 The mean relative error between the responses calculated by 2D-CNN and 3D-FEM

(22)

式中，RLAi_2D-CNN為2D-CNN計(jì)算結(jié)果，RLAi_3D-FEM為3D-FEM計(jì)算結(jié)果.

可以發(fā)現(xiàn)，整體而言，兩者平均相對誤差小于10%，對于大多數(shù)層而言，相對誤差小于3%，滿足精度要求.然而在薄層處，由于受到圍巖層厚、侵入和各向異性的影響，平均相對誤差普遍大于5%，如第2層、第10層和第18層等.該現(xiàn)象正是造成反演結(jié)果精度較低的原因之一.為此，在后續(xù)研究中，我們將繼續(xù)完善訓(xùn)練模型，提升神經(jīng)網(wǎng)絡(luò)算法對圍巖層特征的拾取能力.

值得注意的是，本文所采用的反演模型均為侵入模型，由于多解性問題的存在，會使得無侵地層的侵入深度反演結(jié)果不為0，但侵入帶電阻率和原狀地層電阻率的反演結(jié)果一致性較高，如模型第五層和第九層.整體而言，該方案可有效提升算法克服反演多解性的能力，進(jìn)而提高傾斜各向異性地層陣列側(cè)向測井反演精度，給陣列側(cè)向測井實(shí)際資料處理提供一種新思路.

5 結(jié)論

針對傾斜各向異性地層，本文基于3D-FEM，分析了侵入深度、各向異性以及地層傾角等因素對陣列側(cè)向測井響應(yīng)的影響；接著，引入2D-CNN網(wǎng)絡(luò)架構(gòu)，針對侵入深度小于1.0 m，各向異性系數(shù)小于4，層厚小于5 m，圍巖和地層水平電阻率位于0.1～500 Ωm的地層，通過建立大型數(shù)據(jù)集以及將模型可視化，使得2D-CNN訓(xùn)練精度和計(jì)算速度能夠滿足正反演需求，訓(xùn)練精度達(dá)到99%左右，計(jì)算一個測井點(diǎn)僅需0.36 ms；在此基礎(chǔ)上，通過引入全局尋優(yōu)能力較強(qiáng)的MPGA算法，可為L-M反演算法提供較為精確的初值，進(jìn)而基于區(qū)域分解技術(shù)使得斜井各向異性地層電阻率反演精度達(dá)到90%以上.

不可否認(rèn)的是，在更為復(fù)雜的三維地層條件下，例如碳酸鹽巖地層，存在裂縫和溶洞發(fā)育、地層各向異性強(qiáng)以及電阻率大等特征，本文提供的數(shù)據(jù)庫及網(wǎng)絡(luò)架構(gòu)將不再適用，讀者僅需將其進(jìn)行擴(kuò)充與完善即可.另外，復(fù)雜三維地層的測井響應(yīng)會使得反演多解性更為強(qiáng)烈，為此，在今后研究中，應(yīng)進(jìn)一步豐富和完善神經(jīng)網(wǎng)絡(luò)架構(gòu)以及反演方案，以期為油氣田開發(fā)提供更為便捷可靠的處理方案.