粒子群優(yōu)化交叉偶極子聲波測(cè)井地層各向異性的反演方法

李鵬舉， 魏雙寶， 田 甜

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

0 引 言

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

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

1 地層各向異性反演原理

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

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

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

(1)

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

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

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

(2)

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

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

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

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

(3)

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

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

T——處理的時(shí)間寬度；

s2——慢橫波的慢度。

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

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

(4)

式中：Anis——地層各向異性大?。?/p>

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

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

2 地層各向異性反演方法

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

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

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

設(shè)在一個(gè)D維的搜索空間中由m個(gè)粒子組成的種群，其中第i個(gè)粒子位置為xi=(xi1,xi2,…,xiD)，其速度為vi=(vi1,vi2,…,viD)。第i個(gè)粒子搜索到的最優(yōu)位置為pi=(pi1,pi2,…,piD)。整個(gè)種群搜索到的最優(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。此外，為防止粒子速度過(guò)大，可設(shè)定速度上限vmax，即當(dāng)vid>vmax時(shí)，取vid=vmax；當(dāng)vid<-vmax時(shí)，取vid=-vmax。粒子群的初始位置和初始速度都是隨機(jī)產(chǎn)生的，按式(5)和(6)進(jìn)行迭代，直至兩代全局最優(yōu)位置之差小于給定精度。

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

(7)

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

2.2 反演方法的具體步驟

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

(2)根據(jù)實(shí)際地層給定慢橫波慢度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í)將每個(gè)個(gè)體的初始位置設(shè)為當(dāng)前每個(gè)粒子的最優(yōu)位置pi，把所有粒子的最優(yōu)位置作為當(dāng)前整個(gè)種群的最優(yōu)位置pg。

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

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

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

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

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

3 實(shí)際測(cè)井資料的處理

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

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

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

4 結(jié)束語(yǔ)

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