李鵬舉， 魏雙寶， 田 甜

(1.東北石油大學(xué) 地球科學(xué)學(xué)院， 黑龍江 大慶 163318； 2.東北石油大學(xué) 非常規(guī)油氣成藏與開發(fā)省部共建國家重點實驗室培育基地， 黑龍江 大慶 163318)

0 引 言

交叉偶極子聲波測井作為測量儲層不同方向上聲波特性的測井方法，自20世紀(jì)90年代初出現(xiàn)以來對聲波測井技術(shù)的發(fā)展與應(yīng)用具有十分重要的意義，利用交叉偶極子測井?dāng)?shù)據(jù)計算各向異性參數(shù)也成為了一種重要手段。橫波在各向異性地層中傳播時會分裂成與質(zhì)點偏振方向相互垂直的快橫波和慢橫波，由此可以分析出地層的各向異性。1986年R.M.Alford[1]提出的旋轉(zhuǎn)分析法成功得到了地層的各向異性方位角，但是該方法在處理各向異性較小的地層時會出現(xiàn)多解的情況，多解性是由于兩主波陣列分別處理時的速度誤差之和可能會大于對應(yīng)的快橫波與慢橫波的速度之差。1999年X.M.Tang等[2]提出了波形反演的方法來反演地層的各向異性，該方法能夠大大壓制或消除波形數(shù)據(jù)中的噪音影響，但是由于波形反演算法數(shù)據(jù)量大，反演函數(shù)又具有多個極值點，所以快速、準(zhǔn)確地求取全局最優(yōu)解成為了關(guān)鍵。國內(nèi)學(xué)者以波形反演法公式為目標(biāo)函數(shù)提出了一些相關(guān)的方法。2001年陶果等[3]利用模擬退火法求解偶極橫波時差。2005年何峰江等[4]將改進(jìn)的模擬退火算法應(yīng)用于正交偶極子各向異性反演中。2007年王才志等[5]利用極快速模擬退火算法反演地層各向異性。2010年劉宇[6]綜合極快速模擬退火法和變區(qū)間局部模擬退火法進(jìn)行各向異性計算，均取得了一些應(yīng)用效果。

粒子群優(yōu)化(particle swarm optimizer，PSO)算法是當(dāng)前研究與應(yīng)用最廣泛的群智能算法之一，由J.Kennedy和C.Eberhart在1995年提出[7]。PSO算法在求解最優(yōu)化問題，特別是非線性、多峰、不可導(dǎo)等復(fù)雜優(yōu)化問題中具有較強(qiáng)的全局搜索能力，現(xiàn)已成功應(yīng)用于函數(shù)優(yōu)化、信號處理、神經(jīng)網(wǎng)絡(luò)訓(xùn)練、模式分類、數(shù)據(jù)挖掘等多學(xué)科領(lǐng)域。目前，在地球物理測井方面的應(yīng)用還很少，2014年江沸菠[8]通過將粒子群優(yōu)化算法與BP神經(jīng)網(wǎng)絡(luò)相結(jié)合反演非線性電阻率成像，2016年周超[9]對改進(jìn)的粒子群優(yōu)化算法識別裂縫屬性進(jìn)行了研究，2018年蔡偉等[10]提出應(yīng)用粒子群優(yōu)化算法反演瑞雷波頻散曲線。文中提出了將粒子群優(yōu)化算法應(yīng)用于提取交叉偶極子聲波測井地層各向異性參數(shù)的方法，由于充分利用了PSO算法全局尋優(yōu)、收斂速度快、求解精度高的特點，因此該方法能夠得到全局最優(yōu)解，在各向異性比較弱的地層仍然適用，計算的結(jié)果準(zhǔn)確可靠。

1 地層各向異性反演原理

1.1 Alford旋轉(zhuǎn)法提取快、慢橫波原理

交叉偶極子聲波測井采集4個分量的聲波陣列數(shù)據(jù)，如圖1所示。單個源距的接收器所測量的4個波形信號分別為

圖1 各向異性地層中四分量交叉偶極子測井示意

(1)

式中，F(xiàn)p(t)和Sp(t)分別代表快、慢橫波。把式(1)寫成矩陣形，即

將帶有Fp(t)和Sp(t)的矩陣移到等式左邊，可以得到

因此，可以得到快、慢橫波波形的表達(dá)式為

(2)

1.2 反演目標(biāo)函數(shù)的構(gòu)建

通過上述分析得到了快、慢橫波波形表達(dá)式。為了增加反演時目標(biāo)函數(shù)對θ角變化的敏感性，對式(2)中的角度θ求導(dǎo)數(shù)，得到兩個輔助主波表達(dá)式為

由于分裂的快、慢橫波的極性相同和波形相似，所以可以通過相同和不同接收器間的快、慢橫波移動來完成匹配，如圖2所示，先把接收器n上的慢橫波傳播到接收器m的位置，再在時間上向前移動一定時間，來匹配該位置上的快橫波。即有

其中，m和n表示第m個和第n個接收器;Δz是相鄰兩個接收器的間距;Δs為快、慢橫波的慢度差；zm是第m個接收器到聲源的距離。然后對匹配誤差進(jìn)行累加可得到如下的反演目標(biāo)函數(shù)：

(3)

式中：E——目標(biāo)函數(shù)；

N——接收器的總數(shù)；

T——處理的時間寬度；

s2——慢橫波的慢度。

當(dāng)目標(biāo)函數(shù)值達(dá)到全局最小時，在處理時窗T內(nèi)快、慢橫波才能得到最佳匹配，此時的θ角為快橫波方位角。

反演出Δs、s2和θ以后，就可以得到地層各向異性，其中θ為各向異性方向，各向異性大小按式(4)計算：

(4)

式中：Anis——地層各向異性大??；

s1、s2——快、慢橫波的慢度。

圖2 接收器陣列快、慢橫波匹配處理示意

2 地層各向異性反演方法

如上所述，地層各向異性反演就是采用一定的優(yōu)化方法尋找目標(biāo)函數(shù)(3)的最小值。文中提出采用粒子群優(yōu)化算法對目標(biāo)函數(shù)(3)進(jìn)行優(yōu)化求解，同時反演出各向異性的方位與大小。

2.1 粒子群優(yōu)化算法

粒子群優(yōu)化算法是一種穩(wěn)定、適用性強(qiáng)的全局優(yōu)化算法，它是通過個體間的協(xié)作與競爭實現(xiàn)復(fù)雜空間中最優(yōu)解的搜索。PSO算法是在可行解空間中初始化一群粒子，在每一代中粒子通過追隨兩個極值而動，一個是粒子本身迄今找到的最優(yōu)解，另一個是種群迄今找到的最優(yōu)解。其數(shù)學(xué)描述如下：

設(shè)在一個D維的搜索空間中由m個粒子組成的種群，其中第i個粒子位置為xi=(xi1,xi2,…,xiD)，其速度為vi=(vi1,vi2,…,viD)。第i個粒子搜索到的最優(yōu)位置為pi=(pi1,pi2,…,piD)。整個種群搜索到的最優(yōu)位置為pg=(pg1,pg2,…,pgD)。按追隨當(dāng)前最優(yōu)粒子的原理，粒子每一維速度和位置都按照式(5)、(6)進(jìn)行變化：

(5)

(6)

其中，d=1,2,…，D，D為維數(shù)；i=1,2,…，m，m為種群規(guī)模；t為當(dāng)前進(jìn)化代數(shù)；r1和r2為分布于[0,1]之間的隨機(jī)數(shù)；c1和c2為加速常數(shù)，是粒子向自身極值和全局極值推進(jìn)的加速權(quán)數(shù)，文中設(shè)置為2。此外，為防止粒子速度過大，可設(shè)定速度上限vmax，即當(dāng)vid>vmax時，取vid=vmax；當(dāng)vid<-vmax時，取vid=-vmax。粒子群的初始位置和初始速度都是隨機(jī)產(chǎn)生的，按式(5)和(6)進(jìn)行迭代，直至兩代全局最優(yōu)位置之差小于給定精度。

為了更好的控制算法的探測和開發(fā)能力克服PSO算法在后期搜索速度慢等缺點，Y.Shi等[11]在式(7)的基礎(chǔ)上引入了慣性權(quán)重w，即速度更新公式變?yōu)?/p>

(7)

由式(6)和(7)組成的迭代算法為標(biāo)準(zhǔn)粒子群優(yōu)化算法。慣性權(quán)重w表示粒子上一代速度對當(dāng)前代速度的影響，加強(qiáng)了PSO算法的全局搜索能力?？刂苭取值大小可調(diào)節(jié)PSO算法的全局與局部尋優(yōu)能力。w值較大，全局尋優(yōu)能力強(qiáng)，局部尋優(yōu)能力弱；反之，則局部尋優(yōu)能力增強(qiáng)，而全局尋優(yōu)能力減弱[12]。針對文中問題，通過試驗得出設(shè)置w為0.9時，處理效果較好。

2.2 反演方法的具體步驟

(1)數(shù)據(jù)預(yù)處理：對波形數(shù)據(jù)進(jìn)行去增益、濾波及均衡化處理。

(2)根據(jù)實際地層給定慢橫波慢度s2和快慢橫波慢度差Δs范圍，例如，慢橫波慢度s2為131.2～393.6 μs/m、快慢橫波慢度差Δs為0～65.6 μs/m，快橫波方位角θ為0～180°。

(3)初始化：在給定范圍內(nèi)隨機(jī)生成初始位置xi=(Δsi,θi,s2i)及其速度vi=(vΔsi,vθi,vs2i)，其中i=1,2，…，m，m=50為種群規(guī)模。同時將每個個體的初始位置設(shè)為當(dāng)前每個粒子的最優(yōu)位置pi，把所有粒子的最優(yōu)位置作為當(dāng)前整個種群的最優(yōu)位置pg。

(4)讀取某一深度點的8道四分量波形數(shù)據(jù)，根據(jù)每個粒子的位置按照式(3)計算各向異性目標(biāo)函數(shù)E。

(5)評價每一個粒子：如果粒子函數(shù)值小于其歷史最小值，則將pi設(shè)置為該粒子的位置，并更新個體極值。如果所有粒子的個體極值中最小值小于當(dāng)前的全局極值，則將pg設(shè)置為該粒子的位置，并更新全局極值。

(6)粒子的更新：對每一個粒子的速度和位置按照式(6)和(7)進(jìn)行更新。

(7)檢驗是否結(jié)束：如果當(dāng)前迭代次數(shù)達(dá)到了給定的最大迭代次數(shù)或兩代全局最優(yōu)位置之差小于給定精度ε，則停止迭代，結(jié)束該深度點的計算，輸出快橫波方位角度、快慢橫波慢度及各向異性大??；否則迭代次數(shù)t=t+1，轉(zhuǎn)到步驟(5)。

(8)重復(fù)步驟(3)到步驟(7)，直到整個深度段全部處理完畢。

3 實際測井資料的處理

為進(jìn)一步驗證文中方法的有效性，處理了多口交叉偶極子陣列聲波測井(XMAC)資料，其中圖3為某井的處理成果，選取2 900～2 950 m深度段。圖中第1道為伽馬測井曲線和井徑曲線；第2道為深度道；第3道是同向分量(XX)的變密度波形；第4道是各向異性大小，其中紅色為文中方法計算的各向異性大小，藍(lán)色為CIFLog軟件快速模擬退火法計算的各向異性大小；第5、6道分別為兩種方法計算的快、慢橫波慢度曲線；第7道為兩種方法計算的快橫波方位角曲線。

圖3 解釋成果與各向異性參數(shù)對比

從圖3中可以看出，在各向異性大小Anis較小，即地層各向異性較弱的位置，相比于快速模擬退火法，文中計算的快橫波方位角連續(xù)穩(wěn)定，因此文中方法在弱各向異性地層的計算結(jié)果更穩(wěn)定可靠。在多口井的地層各向異性處理上，文中方法和快速模擬退火法的計算結(jié)果基本一致，且文中方法對弱各向異性地層處理的穩(wěn)定性和計算速度都有所提高，驗證了文中方法求解各向異性參數(shù)的有效性。

4 結(jié)束語

筆者提出并實現(xiàn)了一種利用交叉偶極子聲波測井資料反演地層各向異性參數(shù)的方法。該方法通過交叉偶極子聲波測井四分量數(shù)據(jù)推導(dǎo)出快、慢橫波表達(dá)式，將所有接收器的波形進(jìn)行匹配構(gòu)建反演目標(biāo)函數(shù)，利用粒子群優(yōu)化算法求解反演目標(biāo)函數(shù)的最小值，得到各向異性大小及方向。實際測井?dāng)?shù)據(jù)處理結(jié)果表明，文中提出的粒子群優(yōu)化算法反演交叉偶極子聲波測井地層各向異性參數(shù)的方法是準(zhǔn)確、有效、可靠的，并且提高了計算速度及處理結(jié)果的穩(wěn)定性。因此，該方法具有較高的實用價值，為交叉偶極子聲波測井?dāng)?shù)據(jù)各向異性反演提供了新的思路和方法。