李 勤，王 瑋，王瀚林

(西安科技大學(xué) 地質(zhì)與環(huán)境學(xué)院，陜西 西安 710054)

HTI(Transverse Isotropy with a Horizontal axis of symmetry)介質(zhì)是一種具有近似水平對(duì)稱軸的橫向各向同性介質(zhì)，對(duì)于含有垂直裂隙的模型，可將其近似等效為HTI介質(zhì)模型[1]。HTI介質(zhì)是一種典型的各向異性模型，THOMSEN[2]提出通過(guò)各向異性參數(shù)來(lái)描述介質(zhì)的各向異性程度，并推導(dǎo)了各向異性參數(shù)與彈性參數(shù)之間的量化關(guān)系;董守華[3]針對(duì)煤層進(jìn)行了各向異性參數(shù)測(cè)試，討論煤層各向異性參數(shù)與孔隙率的關(guān)聯(lián)，并指出通過(guò)煤層各向異性的大小和方向能夠預(yù)測(cè)煤層的裂隙密度和方位;陳同俊等[4]利用介質(zhì)等效模型理論分析了構(gòu)造煤的各向異性AVO(Amplitude Variation with Offset)特征，對(duì)不同的AVO近似方程應(yīng)用于煤層的效果進(jìn)行比較，對(duì)AVO近似方程進(jìn)行改進(jìn)，使其能適用于各向異性煤層。彭蘇萍和畢銀麗[5]指出煤層在開采過(guò)程中會(huì)產(chǎn)生裂隙。許多專家研究了各向異性參數(shù)與裂隙之間的關(guān)聯(lián)，SCHOENBERG和DOUMA[6]針對(duì)含平行裂隙的介質(zhì)，將其等效為無(wú)限薄無(wú)限柔順的薄層，忽略裂隙的形狀及其空間分布細(xì)節(jié)，假設(shè)裂隙面處位移張量不連續(xù)，但牽引力張量連續(xù)，且2者線性相關(guān)，提出了線性滑動(dòng)理論，建立了裂隙介質(zhì)等效模型;BAKULIN等[7]提出用裂隙弱度來(lái)描述裂隙，分析裂隙流體因子與裂隙充填物的關(guān)聯(lián);QUINTAL[8]研究了飽和流體研究中流體飽和度對(duì)反射系數(shù)的影響;陳懷震等[9]通過(guò)分析流體因子對(duì)不同流體的響應(yīng)關(guān)系，指出HTI介質(zhì)中裂隙飽和水時(shí)，裂隙流體因子值偏小，趨近于0;郎玉泉等[10]利用Gassmann方程進(jìn)行流體替代，通過(guò)AVO正演，估計(jì)頂板砂巖孔隙度和干濕性。AVO反演的基礎(chǔ)是反射系數(shù)近似式，RüGER[11]首先給出不同的橫向各向同性介質(zhì)的縱波反射系數(shù)近似公式，近似式精度較高，是目前AVO理論的基礎(chǔ)[12-13]。利用反射系數(shù)近似公式，可以實(shí)現(xiàn)地震疊前AVO反演[14-16]。MA[17]利用模擬退火法，通過(guò)疊前AVO反演得到了巖石特性;MISRA和SACCHI[18]利用快速模擬退火法解決了基于邊界保護(hù)的AVO反演問(wèn)題，粒子群算法最早由Eberhart提出，源于對(duì)鳥群捕食行為的擬態(tài)[19]，MIKKI[20]將這種方法用于電磁優(yōu)化，使該算法進(jìn)入工程物探領(lǐng)域并得到發(fā)展，SUN等[21]基于粒子群算法進(jìn)行儲(chǔ)層預(yù)測(cè)，收到一定成效;SUN和LIU[22]利用混沌量子粒子群算法求解了疊前AVO彈性參數(shù);MEHRGINI等[23]利用粒子群算法對(duì)具有中子孔隙度、體積密度、水飽和度時(shí)的地層橫波速度VS進(jìn)行了反演并論證;WU等[24]利用粒子群算法對(duì)地層的縱橫波速度、密度進(jìn)行了反演;嚴(yán)哲等[25]利用量子粒子群算法對(duì)Marmousi2模型的縱橫波速度和密度進(jìn)行了反演，反演結(jié)果穩(wěn)定。對(duì)反演算法進(jìn)行改進(jìn)和優(yōu)化是提高反演精度和反演效率的主要手段，對(duì)于疊前AVO反演來(lái)說(shuō)，粒子群算法(PSO,Particle Swarm Optimization)雖然收斂速度快、原理簡(jiǎn)單，但卻容易早熟，陷入局部最優(yōu);模擬退火法(SA,Simulated Annealing)雖然全局尋優(yōu)能力強(qiáng)但其穩(wěn)定性又不能得到保證。因此，在粒子群算法的基礎(chǔ)上結(jié)合模擬退火的思想，對(duì)粒子群算法進(jìn)行改進(jìn)，使新算法有強(qiáng)行跳出局部最優(yōu)解的能力，這種優(yōu)化后的混合算法在反演時(shí)具有明顯的優(yōu)越性和有效性。筆者基于HTI介質(zhì)反射系數(shù)近似公式，建立關(guān)于各向異性參數(shù)的目標(biāo)函數(shù)，采用基于模擬退火的優(yōu)化粒子群算法(SA-PSO)實(shí)現(xiàn)HTI介質(zhì)疊前AVO反演，并通過(guò)裂隙流體因子與各向異性參數(shù)的關(guān)聯(lián)，實(shí)現(xiàn)流體因子預(yù)測(cè)。

1 HTI介質(zhì)各向異性AVO反演理論及步驟

1.1 基于模擬退火的優(yōu)化粒子群算法反演原理

粒子群算法是通過(guò)模擬鳥類覓食的過(guò)程，達(dá)到解決優(yōu)化問(wèn)題的方法，算法收斂速度快。在粒子群算法中，鳥類被抽象為沒(méi)有質(zhì)量和體積的微粒，并延伸到n維空間，粒子i在n維空間的位置表示為向量Xi=(x1,x2,x3,…,xn)，飛行速度表示為向量Vi=(v1,v2,v3,…,vn)，粒子不僅知道自己發(fā)現(xiàn)的最好位置Pb，也知道整個(gè)群體中的其他粒子發(fā)現(xiàn)的最優(yōu)位置Pb,g，粒子們根據(jù)這些信息來(lái)決定接下來(lái)的行動(dòng)方向。根據(jù)目標(biāo)函數(shù)的構(gòu)建，每個(gè)粒子都有其適應(yīng)值范圍。粒子群優(yōu)化算法首先將粒子初始化為一群隨機(jī)粒子，然后通過(guò)迭代找到最優(yōu)解，在每一次迭代中，粒子通過(guò)跟蹤兩個(gè)“最優(yōu)解”(即Pb,Pb,g)來(lái)更新自己，在找到“最優(yōu)解”(Pb,Pb,g)后，粒子通過(guò)下式更新自己的位置:

(1)

式中，Vi(t)為粒子i在t時(shí)刻的速度;ω為慣性因子;c1和c2為學(xué)習(xí)因子，分別為粒子的局部學(xué)習(xí)能力和全局學(xué)習(xí)能力;r1和r2為[0，1]的隨機(jī)數(shù);Pb,i(t),Pb,gi(t)分別為粒子i和群體中其他粒子在t時(shí)刻為止發(fā)現(xiàn)的最好位置;xi(t)為粒子i在t時(shí)刻的位置。

在搜索過(guò)程中，粒子的尋優(yōu)范圍局限在Pb附近，因此不能保證全局收斂得到全局最優(yōu)解，而慣性因子ω表示的就是繼承先前粒子速度的能力，當(dāng)ω較大時(shí)，繼承粒子前一時(shí)刻速度較多，全局尋優(yōu)能力會(huì)有所提高。ω的值如果過(guò)大或過(guò)小就會(huì)造成算法陷入局部尋優(yōu)或者無(wú)法找出最優(yōu)解。因此，對(duì)慣性因子進(jìn)行實(shí)時(shí)調(diào)整很有必要。針對(duì)ω在算法運(yùn)行各個(gè)階段的不同作用，加入隨迭代次數(shù)非線性減小的慣性因子。減小的慣性因子公式為

ω(k)=ωs-(ωs-ωe)(k/Tmax)2

(2)

式中，ωs為初始慣性因子;ωe為迭代至最大次數(shù)時(shí)的慣性因子;k為當(dāng)前迭代次;Tmax為最大迭代次數(shù)。

學(xué)習(xí)因子c1和c2影響算法的全局和局部搜索能力，c1為粒子的個(gè)體認(rèn)知，因此影響局部搜索能力;c2為群體認(rèn)知，影響全局搜索能力。在粒子群優(yōu)化算法迭代過(guò)程中，迭代前期需要粒子群優(yōu)化算法的全局搜索能力強(qiáng)，隨著迭代進(jìn)行至后期，算法需要快速收斂，此時(shí)需要局部搜索能力強(qiáng)。因此，在粒子群優(yōu)化算法加入自適應(yīng)的學(xué)習(xí)因子使得算法找到更優(yōu)解乃至最優(yōu)解。學(xué)習(xí)因子自適應(yīng)公式為

(3)

式中，R1,R2,R3,R4為學(xué)習(xí)因子的自適應(yīng)度參數(shù)。

在粒子群算法中，粒子在發(fā)現(xiàn)的當(dāng)前最優(yōu)位置靠近局部最優(yōu)解的時(shí)候，所有的粒子就會(huì)聚集到最小的區(qū)域，全局搜索能力就會(huì)削弱。而模擬退火法的優(yōu)勢(shì)正在于通過(guò)調(diào)整溫度，控制概率性跳出特性，進(jìn)行新一次的隨機(jī)搜索，將模擬退火的這種特性引入到粒子群算法，可以有效避免陷入局部最優(yōu)解，這種修正主要是針對(duì)全局最優(yōu)位置Pb,g，利用模擬退火的跳出機(jī)制來(lái)選擇每個(gè)粒子的當(dāng)前最優(yōu)位置，即通過(guò)控制模擬退火法對(duì)“變差解”的接受概率，可以提高粒子群算法的靈活性，增加粒子的多樣性。模擬退火法通過(guò)在求解區(qū)間隨機(jī)游走，利用Metropolis抽樣準(zhǔn)則使隨機(jī)搜索最新位置逐漸收斂于最優(yōu)解。Metropolis準(zhǔn)則定義了某一溫度下系統(tǒng)狀態(tài)從狀態(tài)i到狀態(tài)j的能量概率為

(4)

在模擬退火的跳出機(jī)制下，認(rèn)為Pb是和Pb,g懸殊很大的特殊解，可以計(jì)算出溫度為T時(shí)Pb相對(duì)Pb,g的能量概率，即exp(-fPb-fPb,g)T，其中f為目標(biāo)函數(shù)解值。若將所求概率作為Pb的適配值，則用Pb替代Pb,g。

1.2 HTI介質(zhì)反射系數(shù)與裂隙流體因子

HTI介質(zhì)通常是由平行排列的垂直裂隙形成，即對(duì)稱軸接近于水平方向時(shí)的TI介質(zhì)。HTI介質(zhì)彈性矩陣:

(5)

由于HTI介質(zhì)在垂直面彈性參數(shù)相同，根據(jù)矩陣對(duì)稱性，其彈性參數(shù)有如下關(guān)系:

(6)

根據(jù)Thomson各向異性介質(zhì)理論，HTI介質(zhì)的各向異性參數(shù)和彈性參數(shù)可由式(7)進(jìn)行計(jì)算

(7)

其中，VP為縱波速度;ρ為密度;VS為橫波速度;ε(v),δ(v),γ(v)為HTI介質(zhì)的Thomsen各向異性參數(shù)，ε(v)為縱波各向異性程度，δ(v)為縱波在橫向和垂向之間各向異性變化的快慢程度，γ(v)為快、慢橫波速度的差異程度。

基于Thomsen各向異性參數(shù)，RüGER給出了HTI介質(zhì)的PP波反射系數(shù)公式:

(8)

根據(jù)線性滑動(dòng)模型理論中應(yīng)力與裂隙應(yīng)變的關(guān)系，裂隙介質(zhì)彈性系數(shù)矩陣C可以表示成各向同性背景系數(shù)矩陣Ciso加上各向異性擾動(dòng)Cani之和，則描述TI介質(zhì)的裂隙等效彈性矩陣為

C=Ciso+Cani

(9)

其中，令λ+2μ=L，則Ciso,Cani分別為

(10)

(11)

設(shè)g=μ/(λ+2μ)=(VS/VP)2為各向同性巖石骨架橫波速度和縱波速度比值的平方，經(jīng)推導(dǎo)得到裂隙剛度矩陣中的ΔN和ΔT表達(dá)式:

(12)

式中，ΔN和ΔT為Schoenberg線性滑動(dòng)理論中的裂隙法向弱度和切向弱度，其絕對(duì)值在0～1;e為裂隙密度;g為橫縱波速比的平方;Kf和μf為裂隙充填流體體積模量和剪切模量;λ和μ為不含裂隙巖石的拉梅系數(shù);a/c為裂隙縱橫比，a,c分別為裂隙縱向長(zhǎng)度，橫向?qū)挾取?/p>

結(jié)合Thomsen各向異性參數(shù)和線性滑動(dòng)模型等效理論，經(jīng)推導(dǎo)可以得到HTI介質(zhì)各向異性參數(shù)與裂隙柔度參數(shù)之間的表達(dá)式為

(13)

所以，ΔN和ΔT可由各向異性參數(shù)表示為

一直以來(lái)，產(chǎn)業(yè)界人士都在爭(zhēng)議自動(dòng)化技術(shù)能否永久性、大規(guī)模取代人工勞動(dòng)。按照經(jīng)濟(jì)學(xué)中的“比較優(yōu)勢(shì)”理論，即便技術(shù)進(jìn)步，人類也仍然會(huì)在很多領(lǐng)域保持優(yōu)勢(shì)。因此，技術(shù)不會(huì)取代人工勞動(dòng)，而是釋放了工人，讓人類可以從事不具有危險(xiǎn)性，但更具挑戰(zhàn)性的工作。但要承認(rèn)，機(jī)器沒(méi)有人類的弱點(diǎn)和偏見(jiàn)，更加不偏不倚、不主觀臆斷。最重要的是，機(jī)器記憶和處理數(shù)據(jù)的精確度遠(yuǎn)遠(yuǎn)超出人類。而數(shù)據(jù)正以指數(shù)級(jí)增長(zhǎng)。

(14)

假設(shè)裂隙之間互不連通、互不影響,HTI型裂隙介質(zhì)中流體識(shí)別可以用KN/KT來(lái)指示:

(15)

其中，KN為裂隙法向柔度;KT為裂隙切向柔度。對(duì)于液體充填的裂隙介質(zhì)，剪切模量μf=0，體積模量Kf與μ大致相差一個(gè)數(shù)量級(jí)，裂隙縱橫比通常很小,a/c?1，[(Kf+4/3μf)/μ]/(a/c)?1，此時(shí)ΔN≈0，ΔT與裂隙為干燥狀態(tài)時(shí)相同，流體因子KN/KT接近于0。

1.3 HTI介質(zhì)各向異性反演實(shí)現(xiàn)過(guò)程

采用SA-PSO算法對(duì)HTI介質(zhì)進(jìn)行疊前AVO反演，基本思想可概括為:將待反演模型的上、下層的彈性參數(shù)作為種群中尋找最優(yōu)解的粒子，通過(guò)位置公式進(jìn)行迭代計(jì)算，選用不同的慣性因子和學(xué)習(xí)因子，使目標(biāo)函數(shù)求解達(dá)到全局最小。具體流程如下:

(1)將模型確立為HTI介質(zhì)模型;

(2)通過(guò)地震波反射系數(shù)近似公式計(jì)算出模型的PP波反射系數(shù);

(3)確定HTI介質(zhì)模型的輸入?yún)?shù)及待反演參數(shù)，構(gòu)建目標(biāo)問(wèn)題的反演目標(biāo)函數(shù);

(4)將先驗(yàn)信息與待反演參數(shù)輸入反演目標(biāo)函數(shù);

(5)輸入先驗(yàn)地震數(shù)據(jù)信息或?qū)崪y(cè)地震數(shù)據(jù)，利用SA-PSO算法尋找反演目標(biāo)函數(shù)的最佳擬合;

(6)輸出待反演參數(shù)的最優(yōu)解，預(yù)測(cè)儲(chǔ)層的相關(guān)流體因子。

假設(shè)待反演的各向異性參數(shù)符合粒子群分布，將反演的約束條件設(shè)為

(16)

其中，X=(δ(v),ε(v),γ(v))為待反演的各向異性參數(shù);i為地震記錄的道序號(hào);dobs(θi)=[R(Sobs(θi))，I(Sobs(θi))]，dsyn(θi)=[R(Ssyn(θi))，I(Ssyn(θi))]，R和I分別為反射系數(shù)與子波褶積后復(fù)數(shù)的實(shí)部與虛部;θi為該道地震波對(duì)應(yīng)的PP波入射角;Sobs(θi)

和Ssyn(θi)分別為HTI介質(zhì)模型在PP波入射角為θi時(shí)的PP波反射系數(shù)序列與子波序列的褶積。

利用上述方法進(jìn)行反演的技術(shù)流程如圖1所示。

圖1 基于PSO-SA算法的反演流程Fig.1 Schematic diagram of inversion flow based on PSO-SA algorithm

2 模型試算

2.1 多層模型試算

建立一個(gè)4層的模型，從地表往下分別為第1，2，3，4層，橫向?yàn)?00 m。其中，第1層為泥巖，厚度約為300 m;第2層為砂巖，含有垂向飽水裂隙，可近似為HTI型介質(zhì)，層厚約為100 m;第3層為煤層，厚度約為10 m;第4層為泥巖，厚度約為300 m。模型對(duì)應(yīng)的縱、橫波速度體如圖2(a)，3(a)所示，模型彈性參數(shù)和各向異性參數(shù)見(jiàn)表1。根據(jù)模型中第2層為基于HTI型飽水裂隙的假設(shè)，其各向異性參數(shù)ε(v)的值無(wú)限接近于0，δ(v)和γ(v)取值在0.1～0.3。

表1 理論模型的物性參數(shù)和各向異性參數(shù)Table 1 Physical and anisotropy parameters of the theoretical model

圖2 基于SA-PSO算法的模型VP反演剖面對(duì)比Fig.2 Comparison of model VP inversion profiles based on SA-PSO algorithm

圖3 基于SA-PSO算法的多層模型VS反演剖面對(duì)比Fig.3 Model values and inversion results of S-wave velocity based on SA-PSO in the model

對(duì)多層模型設(shè)計(jì)目標(biāo)函數(shù)并進(jìn)行彈性參數(shù)和各向異性參數(shù)反演。首先，為證明本研究所提出的基于SA-PSO算法的科學(xué)性和有效性，對(duì)多層模型的VP進(jìn)行不同迭代次數(shù)的SA-PSO反演和單一的PSO反演進(jìn)行對(duì)比，結(jié)果如圖2所示。圖2(a)為模型的VP，對(duì)比圖2(b)可以看出，利用SA-PSO混合算法迭代次數(shù)較多(100次)時(shí)，尋優(yōu)效果更穩(wěn)定，得到的地層分異明顯，反演速度也與理論速度誤差較小。在圖2(c)中，使用SA-PSO迭代50次時(shí)，此時(shí)全局尋優(yōu)尚不穩(wěn)定，精度不高，所以出現(xiàn)細(xì)小噪點(diǎn);圖2(d)為利用PSO反演，迭代次數(shù)達(dá)到100次時(shí)的結(jié)果，可以看出地層分異相較于圖2(c)效果稍好一些，參數(shù)區(qū)間也更精確一些，但是相比圖2(b)來(lái)說(shuō)，反演結(jié)果不夠理想。綜上所述，利用SA-PSO反演時(shí)，由于SA-PSO反演克服了PSO的早熟現(xiàn)象，在迭代次數(shù)達(dá)到一定要求后反演精度高于單一的PSO反演，本模型中，當(dāng)?shù)螖?shù)達(dá)到100次時(shí)反演結(jié)果最好，得到的結(jié)果更穩(wěn)定。

同時(shí)，對(duì)多層模型中的VS,ρ進(jìn)行迭代次數(shù)為100次的SA-PSO反演，反演結(jié)果如圖3,4所示。分析圖3，4可以看出，從VP和VS的反演剖面來(lái)看，模型的地層能被清晰地分辨，反演結(jié)果較好。對(duì)于ε(v)，δ(v)，γ(v)，在初始模型中，由于除了砂巖層之外的其他地層均等效為各向同性介質(zhì)，其各向異性參數(shù)均設(shè)置為0。在對(duì)砂巖層的各向異性參數(shù)進(jìn)行反演時(shí)，抽取了橫向距離200 m處的數(shù)據(jù)進(jìn)行試算，試算結(jié)果如圖5所示。在圖5中，對(duì)各向異性參數(shù)的反演值與理論各向異性參數(shù)對(duì)比，可以看出，反演值始終圍繞理論值波動(dòng)，反演結(jié)果穩(wěn)定，同時(shí)，ε(v)始終在0附近，這符合對(duì)砂巖層含有飽水裂隙的預(yù)設(shè)。多層模型的彈性參數(shù)和各向異性參數(shù)反演剖面均能指示地層，各向異性參數(shù)ε(v)出現(xiàn)在0附近，對(duì)裂隙流體的指示性強(qiáng)。

圖4 基于SA-PSO算法的多層模型ρ反演剖面對(duì)比Fig.4 Model values and inversion results of ρ based on SA-PSO in the model

針對(duì)抽取的橫向距離200 m處的試算結(jié)果，對(duì)其中縱向深度為350 m、橫向距離為200 m處的反演結(jié)果進(jìn)行誤差定量分析，見(jiàn)表2。由表2可知，對(duì)于多層模型的多參數(shù)反演，ε(v)，VP的反演誤差最小，ρ，δ(v)的反演誤差略大。但整體來(lái)講，反演誤差都較小。

圖5 基于SA-PSO算法的多層模型各向異性參數(shù) 反演結(jié)果對(duì)比Fig.5 Model values and inversion results of anisotropic parameters based on SA-PSO in the model

同時(shí)，為討論該反演方法的抗噪性，依照正演流程計(jì)算模型的反射系數(shù)并采用40 Hz的雷克子波與反射系數(shù)進(jìn)行褶積，生成不同入射角的角道集，并分別添加了信噪比(SNR)為5，10，20的高斯隨機(jī)噪聲，并針對(duì)各向異性層，分析不同信噪比數(shù)據(jù)的反演誤差(表3)。從表3來(lái)看，大致可以認(rèn)為信噪比越高的數(shù)據(jù)反演誤差越小，在反演過(guò)程中也能夠越快收斂，但是在模型試算的結(jié)果中，信噪比與反演誤差并不完全負(fù)相關(guān);即使在SNR=5時(shí)，角道集仍然能夠反映巖石物理模型的層位特征，在信噪比較低的地震記錄中也可以讀取有效信息，也就是說(shuō)，信噪比較低的情況下，即使反演誤差稍大，但是對(duì)VP，ε(v)，δ(v)的反演結(jié)果仍然能滿足反演要求，反演算法具有一定的抗噪性。

表2 理論模型的反演誤差結(jié)果分析Table 2 Analysis of inversion error results of theoretical model

表3 有噪聲的理論模型反演結(jié)果誤差分析Table 3 Analysis of inversion error results of theoretical model with noise

利用模型的物性參數(shù)和各向異性參數(shù)反演值，根據(jù)多層模型縱向深度為301～399 m、橫向距離為100 m處的反演結(jié)果，利用SA-PSO算法對(duì)模型中砂巖層的流體指示因子KN/KT進(jìn)行預(yù)測(cè)，結(jié)果如圖6所示。由圖6可知，利用抽取的部分模型物性參數(shù)和各向異性參數(shù)作為先驗(yàn)值，得到的流體指示因子KN/KT均在0.4以下，最優(yōu)值在0.12左右，符合裂隙飽和水時(shí)流體因子較小，接近于0的規(guī)律。

圖6 基于SA-PSO算法預(yù)測(cè)得到的流體因子Fig.6 Fluid factors predicted based on SA-PSO algorithm

2.2 改進(jìn)的Marmousi2部分模型試算

原始Marmousi2模型如圖7所示,是復(fù)雜二維構(gòu)造的地質(zhì)-地球物理模型，屬于各向同性介質(zhì)，本文選取紅框圈定的部分模型，并根據(jù)不同巖層的特性賦予各向異性參數(shù):將含水砂巖之外的薄層作為各向同性介質(zhì)，假定含水砂巖垂向裂隙發(fā)育，將其近似等效為HTI介質(zhì)，且由于HTI介質(zhì)在飽水裂隙發(fā)育情形下，各向異性參數(shù)ε(v)接近0，用于各向異性參數(shù)反演和裂隙流體因子的預(yù)測(cè)。采用主頻為50 Hz的雷克子波序列與HTI介質(zhì)的PP波反射系數(shù)序列進(jìn)行褶積，并對(duì)橫向距離為2 000 m、深度為 500～800 m的剖面，生成角道集如圖8所示。由圖8可知，反射同相軸大致能夠反映所選模型區(qū)域反射界面物性差異特征，反射振幅在裂隙附近呈現(xiàn)正極性，能量較弱。

圖7 Marmousi2模型及抽取區(qū)域示意[28]Fig.7 Marmousi2 model and the selected area[28]

圖8 Marmousi2改進(jìn)模型抽取的角道集(2 000 m處)Fig.8 Synthetic angle gather in the Marmousi2 model (at the distance of 2 000 m)

采用本文研究的反演方法對(duì)這部分?jǐn)?shù)據(jù)進(jìn)行迭代次數(shù)為100次的SA-PSO反演，反演結(jié)果如圖9～12所示。分析圖9～11可以看出，對(duì)改進(jìn)的Marmousi2部分模型反演得到VP，VS和ρ剖面，雖然反演值與理論值有一定誤差，但仍然能清晰反映所抽取的Marmousi2部分模型的地層變化規(guī)律，分層效果明顯，與實(shí)際模型基本吻合，反演結(jié)果較好。如圖12所示的各向異性參數(shù)反演結(jié)果曲線對(duì)比，各向異性參數(shù)ε(v)和γ(v)反演精度略高于各向異性參數(shù)δ(v)，對(duì)反演目標(biāo)方程的敏感度相對(duì)較好。從多組反演結(jié)果與理論值的對(duì)比剖面可以看出，反演結(jié)果能夠較準(zhǔn)確反映原始剖面形態(tài)特征，且與理論模型的擬合度高，誤差比較小。選取模型中橫向距離為2 000 m、深度為 630～660 m的區(qū)域，對(duì)流體指示因子KN/KT進(jìn)行預(yù)測(cè)，所得結(jié)果如圖13所示，由于抽取的區(qū)域均為飽水裂隙區(qū)域，此時(shí)進(jìn)行間隔采樣得到的流體因子就是PSO-SA算法中的“粒子”，“粒子們”最終迭代出最優(yōu)結(jié)果(即紅色較大圓點(diǎn))，因子粒子尋優(yōu)時(shí)其初始值是隨機(jī)的，粒子的總數(shù)量與KN/KT的總數(shù)量一致,但粒子的取值與KN/KT的值不具有一一對(duì)應(yīng)的關(guān)系，圖中所有粒子的值都小于1，與飽水裂隙中KN/KT的取值范圍一致。

圖9 基于SA-PSO算法的 Marmousi2改進(jìn)模型VP反演剖面對(duì)比Fig.9 Model values and inversion results of P-wave velocity based on SA-PSO in the Marmousi2 model

圖10 基于SA-PSO算法的 Marmousi2改進(jìn)模型VS反演剖面對(duì)比Fig.10 Model values and inversion results of S-wave velocity based on SA-PSO in the Marmousi2 model

圖11 基于SA-PSO算法的 Marmousi2改進(jìn)模型ρ反演剖面對(duì)比Fig.11 Model values and inversion results of ρ based on SA-PSO in the improved Marmousi2 model

圖12 基于SA-PSO算法的Marmousi2改進(jìn)模型各向異性 參數(shù)反演結(jié)果對(duì)比Fig.12 Model values and inversion results of anisotropic para- meters based on SA-PSO in the improved Marmousi2 model

圖13 基于SA-PSO算法的Marmousi2改進(jìn)模型 流體因子預(yù)測(cè)Fig.13 Fluid factor prediction of Marmousi2 improved model based on SA-PSO algorithm

3 結(jié) 論

(1)綜合考慮粒子群算法快速收斂的特點(diǎn)和模擬退火的概率性跳出特性，通過(guò)對(duì)兩種算法的混合使用，實(shí)現(xiàn)了利用基于模擬退火的混合粒子群算法對(duì)HTI介質(zhì)的AVO反演。該算法具有全局尋優(yōu)能力較強(qiáng)，克服早熟現(xiàn)象，不易陷入局部極小值等優(yōu)點(diǎn)。

(2)將基于SA-PSO算法的HTI介質(zhì)疊前AVO反演應(yīng)用于多層模型和復(fù)雜模型的反演計(jì)算中，反演結(jié)果與理論模型擬合較好，計(jì)算結(jié)果穩(wěn)定，抗噪性較好，反演精度較高。

(3)根據(jù)裂隙含水時(shí)流體指示因子KN/KT較小的特點(diǎn)，并借助HTI介質(zhì)各向異性反演，通過(guò)建立裂隙流體因子KN/KT與各向異性參數(shù)的關(guān)系，可實(shí)現(xiàn)流體因子預(yù)測(cè)。