基于布谷鳥搜索的地震屬性聚類方法及其在塔中地區(qū)碳酸鹽巖礁灘儲層的應用

曹成寅 高 贊

（1. 核工業(yè)北京地質(zhì)研究院，北京 100029；2. 石油工業(yè)出版社有限公司，北京 100011）

0 引 言

目前可以從地震數(shù)據(jù)中提取幾十甚至上百種地震屬性，每種屬性都可以從不同的方面來反映地下儲層的巖性、物性和含油氣性等相關的地質(zhì)特征，它們與儲層特征之間的關系是非常復雜的。在地震解釋時通常需要從中優(yōu)選出對地質(zhì)和儲層特征比較敏感的屬性或者屬性組合來研究地下地質(zhì)特征，這就是常說的地震屬性分析，其中地震屬性聚類是比較重要的一種地震屬性分析方法。目前聚類技術可以分為5 類方法：劃分的方法、層次的方法、基于密度的方法、基于網(wǎng)格的方法和基于模型的方法。J.Dumay 等［1］首先提出將K均值（K‐means）聚類方法應用于地震相分析；陳遵德等［2］將自組織特征映射（SOM）方法應用于地震波形聚類識別油氣；1997 年Paradigm 公司采用神經(jīng)網(wǎng)絡方法對地震波形進行分類，推出的Strati‐magic 軟件，是目前業(yè)界進行地震沉積相劃分應用最為廣泛的軟件；楊培杰等［3］將模糊C均值（FCM）聚類方法應用于地震多屬性聚類中；F.Liu等［4］和劉杏芳［5］將粒子群優(yōu)化算法應用于地震屬性聚類，在碳酸鹽巖儲層地震相劃分中取得了比較好的效果；印興耀等［6］提出了一種基于核子空間的模糊聚類方法，在儲層預測中較好地刻畫了碳酸鹽巖含氣儲層的邊界。

K均值是地震儲層預測中最常用的聚類方法之一，具有很高的計算效率，但是該方法的不足之處是易陷入局部最優(yōu)，計算結(jié)果依賴于初值，計算不穩(wěn)定。群體智能是自然界生物群體通過合作表現(xiàn)出的智能現(xiàn)象，個體行為準則簡單，但是群體合作表現(xiàn)出極強的全局搜索能力，是一類用于解決復雜問題的優(yōu)化方法。本文將群體智能方法中的布谷鳥算法引入到基于K均值的地震屬性聚類中，克服了傳統(tǒng)的K均值聚類方法的不足，在實際應用中取得了一定的效果。

1 基于布谷鳥搜索的聚類算法

聚類分析是一種重要的數(shù)據(jù)挖掘技術，是在劃分的類未知的情況下，根據(jù)“類內(nèi)盡可能相似，類間盡可能相異”的原則將數(shù)據(jù)對象分組成多個類或者簇，K均值方法為其中最常見的方法。

1.1 K均值算法

K均值算法以k為參數(shù)，把n個樣本分為k個簇，以使簇內(nèi)具有較高的相似度，而簇間具有較低的相似度。相似度的計算根據(jù)樣本與所屬簇中樣本的平均值（被看作簇的中心）的距離來進行。

K均值算法的處理流程：首先，隨機選擇k個樣本，每個樣本代表了一個簇的平均值或中心，對剩余的每個樣本，根據(jù)其與各個簇中心的距離，將它賦給最近的簇，然后重新計算每個簇的平均值。這個過程不斷重復，直到準則函數(shù)收斂。通常采用平方誤差準則，其定義為

式中：E(x)——所有樣本平方誤差的總和；xj——樣本空間中的點；—— 簇Ci的 平 均 值，d(xj，)——樣本點xj到中心點xˉi的距離。

這個準則試圖使生成的結(jié)果簇內(nèi)盡可能地緊湊和簇間盡可能的獨立。當簇是密集的，并且簇與簇之間的區(qū)別明顯時，它的效果較好。對處理大數(shù)據(jù)集，該算法是相對可伸縮的和高效的。但是這個算法經(jīng)常收斂于局部最優(yōu)，另外該算法初始中心點的選擇是隨機的，其結(jié)果依賴于初始值的選擇，每次聚類的結(jié)果不一定相同，計算不穩(wěn)定。

1.2 布谷鳥搜索算法

布谷鳥搜索是一種全局尋優(yōu)算法，它是由劍橋大學X.S.Yang 等［7］根據(jù)一些布谷鳥品種如Guria和ani 的“寄生生育”習性提出來的，該算法用Levy‐Flight 過程取代傳統(tǒng)群體智能優(yōu)化算法中的隨機過程，具有更高的搜索效率。與常規(guī)的優(yōu)化算法相比，布谷鳥搜索不要求矩陣求逆或者目標函數(shù)可導，因此具有更廣泛的應用前景。

布谷鳥中Guria 和ani 是一種很獨特的鳥類，它們具有一種獨特的生育方式——寄生生育后代的方式。布谷鳥通常在它生育后代的時候?qū)⒙旬a(chǎn)在食性、卵的顏色、孵化時間與它相差不大的鳥（稱之為宿主鳥）的巢中，為了防止宿主鳥發(fā)現(xiàn)布谷鳥的卵，增加布谷鳥卵孵化的成功率，布谷鳥同時會移走一些宿主鳥的卵。另外在布谷鳥產(chǎn)完卵后，宿主鳥也有可能發(fā)現(xiàn)外來布谷鳥的卵，發(fā)現(xiàn)后它會將布谷鳥的卵給找出來丟掉，或直接將巢廢棄，重新筑巢產(chǎn)卵［7］。X.S.Yang 通過分析布谷鳥尋巢產(chǎn)卵的過程，總結(jié)出3 條規(guī)則［8‐9］：

（1）假定布谷鳥每次只產(chǎn)一個卵，并隨機尋找巢孵化它；

（2）在每次搜索過程中布谷鳥只會保留它認為最優(yōu)的那個巢；

（3）每個巢的宿主鳥發(fā)現(xiàn)巢中布谷鳥卵的概率是pa。

如果發(fā)生這種狀況，那么宿主鳥就會將布谷鳥的卵扔掉或者直接廢棄巢，重新筑巢產(chǎn)卵。簡單的來說，該條規(guī)則可以被近似為賦予每個巢一個被發(fā)現(xiàn)的概率εi，如果εi>pa，則認為產(chǎn)在該巢的卵被發(fā)現(xiàn)，需要重新尋巢，否則保留該巢。

根據(jù)以上3 個規(guī)則，布谷鳥尋巢產(chǎn)卵的過程可以通過公式表示為

式中：xi，xj，xk——3 個不同的巢；t——迭代次數(shù)；λ——步長，它滿足正態(tài)分布的隨機過程；?——哈達瑪?shù)姆e算子；H(u)——Heaviside 函數(shù)；εi——隨機賦予每個巢的被發(fā)現(xiàn)的概率，服從正態(tài)分布。

式中Heaviside 函數(shù)H(u)是δ函數(shù)從負無窮到u的積分，其具體表達式為

利用式（2）可以對實際問題進行優(yōu)化。然而在實際應用中，往往采用Levy‐Flight 過程來取代服從正態(tài)分布隨機過程，使得算法具有更高的收斂速度。

在智能優(yōu)化算法中采用Levy‐Flight 能擴大搜索范圍，增加種群的多樣性，更容易跳出局部最優(yōu)［10‐14］。因此，在實際計算中常常引入Levy‐Flight隨機過程。Levy‐Flight 是一種步長滿足重尾概率分布的隨機搜索，在行走過程中短距離的搜索和長距離的跳躍相間，能擴大搜索范圍、增加種群多樣性，更容易跳出局部最優(yōu)點。其分布滿足積分的表達式為

式中：s——改進后的步長；β——尺度因子。

除了特殊的例子外，該分布不具有解析式。當β= 1 時該積分就成為柯西分布，當β= 2 就變?yōu)檎龖B(tài)分布。在搜索未知的大尺度的空間時，Levy‐Flight 往往比滿足正態(tài)分布的布朗運動更加高效。那么用Levy‐Fight 表示的布谷鳥尋巢產(chǎn)卵過程可以通過公式表示為

式中L(s，β)——Levy‐Flight 隨機過程。

當s比較大時，其表達式為

式中 Γ(x)——伽馬函數(shù)。

Γ(x)的表達式為

在大多數(shù)情況下，步長通常取0~1。

P.Civicioglu 等［15］證明，在目前的群體智能優(yōu)化算法中，布谷鳥搜索算法具有較高的精度和穩(wěn)定性。因此本文將布谷鳥搜索算法引入到地震屬性聚類中，用以解決K均值算法對初值敏感和易陷入局部收斂的問題。

1.3 基于布谷鳥搜索的地震屬性聚類方法

K均值算法對初始值敏感，并且易收斂于局部最優(yōu)，為了避免這些問題，本文就將全局尋優(yōu)的布谷鳥搜索算法和K均值方法結(jié)合起來，進行地震屬性聚類。其中，每個巢代表一個聚類方案，巢位置代表各類別的中心點位置，適應值可以用K均值的目標函數(shù)的倒數(shù)來表示，即

式中ε——為了保證系統(tǒng)的穩(wěn)定性人為賦予的一個非常小的常量。

基于布谷鳥搜索建立地震屬性聚類方法步驟（圖1）。

圖1 基于布谷鳥搜索的地震屬性聚類方法流程Fig.1 Workflow of seismic attribute clustering method based on Cuckoo Search

（1）從地震數(shù)據(jù)中提取多種屬性，并對各屬性進行預處理；

（2）初始化聚類所需的各種參數(shù)，例如，初始的類別數(shù)、巢的個數(shù)（聚類方案的個數(shù)）和巢位置（聚類的中心點），然后計算每個巢的適應值，尋找出當前適應值最大的巢作為最佳巢；

（3）更新每個巢的位置，并計算每個巢新的適應值；

（4）隨機賦予每個巢1 個概率，比較每個巢的概率和預設的被發(fā)現(xiàn)的概率，如果有1 個巢概率大于被發(fā)現(xiàn)的概率（該聚類方案陷入了局部最優(yōu)），則更新該巢位置；

（5）計算新巢的適應值，選擇適應值最大的巢作為當前最優(yōu)的巢；

（6）判斷是否達到終止條件（最佳巢的適應值不變或是迭代次數(shù)達到預設值），如果是，則終止計算，當前的最佳巢位置即為最優(yōu)聚類方案中各類的中心點，否則轉(zhuǎn)至第3 步繼續(xù)下次迭代計算。

最終輸出的巢位置即為最佳聚類的中心點所在的位置，最終的分類方案就是地震屬性聚類結(jié)果。

2 模型試驗

在研究中采用了加州大學提出的UCI 機器學習數(shù)據(jù)庫中的Iris 和Wine 數(shù)據(jù)集對聚類算法的精度進行測試。Iris 數(shù)據(jù)集以鳶尾花的特征作為數(shù)據(jù)來源，數(shù)據(jù)庫中共150 個樣本，4 種特征分別對應著花瓣和花萼的長度和寬度，圖2 為其中的前3 個特征。這些樣本被分為3 個品種，分別為Setosa（第1 種）、Versicolour （第2 種） 和Virginica （第3種），每種含50 個樣本。在3 個品種中，第1 種與其他兩種能比較容易地分開，而第2 種和第3 種比較難區(qū)分。

圖2 Iris數(shù)據(jù)前3個特征Fig.2 Top 3 characteristics of Iris data

本文中分別用K均值方法和基于布谷鳥搜索算法的方法進行聚類，聚類結(jié)果統(tǒng)計見表1。聚類精度的計算公式為

表1 Iris數(shù)據(jù)集的聚類正確率Table 1 Clustering accuracy rates of Iris data set

式中：r——聚類精度；N——樣本數(shù)；er——錯分的樣本數(shù)。

由表1 可見，對于易分的第1 種花，2 種聚類算法都能正確的劃分，沒有分錯的樣本。對于比較難區(qū)分的第2 種花和第3 種花，基于布谷鳥搜索算法的聚類方法正確率明顯高于K均值聚類方法。對于整個Iris 數(shù)據(jù)集來說，基于布谷鳥搜索算法的聚類方法對所有樣本的分類正確率可以達到92.67%，要明顯高于K均值聚類的84.00%。

Wine 數(shù)據(jù)集源自意大利所產(chǎn)的3 種葡萄酒的178 個樣本，有酒精含量、蘋果酸、灰分等13 種特征，圖3 為前3 個特征。全部樣本被分為3 類，每類分別有59、71 和48 個樣本，分別用K均值方法和基于布谷鳥搜索算法的聚類方法對該數(shù)據(jù)集進行聚類，統(tǒng)計結(jié)果見表2。

圖3 Wine數(shù)據(jù)前3個特征Fig.3 Top 3 characteristics of Wine data

表2 Wine數(shù)據(jù)集的聚類正確率Table 2 Clustering accuracy rates of Wine data set

其中，基于布谷鳥搜索算法的聚類方法和K均值聚類方法都能正確地劃分出第1 類，對于第2類，基于布谷鳥搜索算法的聚類方法的正確率要明顯高于K均值聚類的結(jié)果，對于第3 類兩者的正確率相差不大?？傮w來說，基于布谷鳥搜索算法的聚類方法對所有樣本的分類正確率均能達到95.51%，明顯高于K均值聚類的91.57%。

3 應用實例

研究區(qū)位于中國西部塔里木盆地的塔中Ⅰ號坡折帶上，是一個長期繼承性的古隆起，目的層為晚奧陶世在臺地邊緣發(fā)育的一套縱向上多旋回疊置、橫向上多期次加積的較大規(guī)模的礁灘相沉積體系，巖性主要是泥質(zhì)條帶灰?guī)r和顆?；?guī)r，儲層類型為裂縫性和孔洞性儲層，非均質(zhì)性強，埋深普遍超過5 000 m。圖4 為研究區(qū)過井B8 的一條地震剖面，解釋的紅色與綠色地層界面中間部分為該研究區(qū)的目的層，該研究區(qū)地震數(shù)據(jù)分辨率比較低，地震資料主頻約為20 Hz，儲層識別和預測的難度比較大。研究區(qū)內(nèi)共有8 口完鉆井，其中4 口井鉆遇油層（井B5、井B7 和井B8 為高產(chǎn)油流井，井B2 為低產(chǎn)油流井），其余4 口為干井（井B1、井B3、井B4 和井B6）。

圖4 研究區(qū)過井B8的地震剖面Fig.4 Seismic section crossing Well B8 in studied area

碳酸鹽巖沉積相帶的分布與沉積前的古地貌有很大的關聯(lián)性，目的層殘余厚度可以在一定程度上反映目的層沉積前的古地貌，用于分析有利儲層的發(fā)育部位。圖5 為目的層的殘余時間厚度，其中綠色到藍色區(qū)域代表殘余厚度比較大，為古地貌的低洼地，為臺內(nèi)洼地部位；而紅色反映的是殘余厚度比較小，是古地貌的高部位，是礁灘沉積體系中儲層發(fā)育的有利相帶。

圖5 研究區(qū)目的層的殘余時間厚度Fig.5 Residual time thickness of target layer in studied area

利用該工區(qū)的三維地震數(shù)據(jù)提取了目的層瞬時相對振幅、瞬時頻率、瞬時相位、相干體4 種屬性（圖6）。其中瞬時相對振幅屬性能刻畫礁灘體的發(fā)育區(qū)域，但是邊界刻畫得不是很清楚，而瞬時頻率屬性能比較清楚地刻畫出工區(qū)內(nèi)的油氣分布信息，與鉆井已知的油氣分布信息基本一致，瞬時相位與古地貌特征比較相近，較好地刻畫出了礁灘體分布，相干體以及波形差主要描述了工區(qū)內(nèi)斷裂的分布。通過多次試驗，將該研究區(qū)聚類中的類別數(shù)設置為7 時比較合理。隨后分別采用K均值方法、商業(yè)軟件中的SOM 方法和基于布谷鳥搜索算法對這些地震屬性進行聚類，結(jié)果如圖7―圖9 所示，圖中不同顏色反映了不同的地質(zhì)特征。

圖6 研究區(qū)相關地震屬性Fig.6 Seismic attributes from studied area

圖7 為2 次K均值聚類結(jié)果，這2 次K均值聚類結(jié)果都不太一樣。其中在圖7（a）中，類別數(shù)5代表臺內(nèi)洼地，類別數(shù)1、2 和6 代表礁灘體地層發(fā)育位置。但是該結(jié)果刻畫的礁灘體的邊界不是很清楚，無法較好地刻畫油氣分布。在圖7（b）中，類別數(shù)2 代表臺內(nèi)洼地，類別數(shù)4 和6 表示礁灘體地層發(fā)育的位置，其分布與古地貌圖有較好的匹配，能較好地描述礁灘體的分布，但是由于將斷層和油氣聚集區(qū)劃分為同一類別（類別數(shù)5），因而無法較好地區(qū)分斷裂與油氣分布。圖8 是用某商業(yè)軟件中SOM 方法聚類的結(jié)果，圖中類別數(shù)6為臺內(nèi)洼地，類別數(shù)7 為礁灘體的發(fā)育區(qū)域，但是該方法刻畫的礁灘體的輪廓，與古地貌的特征不相匹配，并且將油氣分布與斷層劃分為一類（類別數(shù)2），高產(chǎn)油井和低產(chǎn)油井的分布區(qū)域劃分在同一類別中，無法進行有效區(qū)分。

圖7 K均值聚類結(jié)果Fig.7 K-means clustering results

圖8 某商業(yè)軟件SOM聚類Fig.8 SOM clustering of a business software

圖9 是基于布谷鳥搜索地震屬性聚類結(jié)果，其中類別數(shù)5 代表臺內(nèi)洼地，類別數(shù)6 和7 為礁灘體的發(fā)育位置，與研究區(qū)古地貌特征非常匹配，類別數(shù)3 所代表的區(qū)域應是斷裂分布的位置，類別數(shù)2 應為油氣聚集區(qū)，高產(chǎn)油井B5、B7 和B8 都是分布在這一區(qū)域中，低產(chǎn)油井B2 分布在類別數(shù)3 所在的區(qū)域中，4口不產(chǎn)油井則是分布在其他類別數(shù)代表的區(qū)域中，說明該結(jié)果能較好地刻畫礁灘相儲層的發(fā)育部位，并且能較好地區(qū)分斷層與油氣分布。

圖9 基于布谷鳥搜索的地震屬性聚類Fig.9 Seismic attribute clustering based on Cuckoo Search(CS)

4 結(jié) 論

（1）將群體智能算法中的布谷鳥搜索算法引入到地震屬性聚類中可以較好地解決K均值聚類局部收斂、計算結(jié)果依賴于初值和不穩(wěn)定的問題。

（2）通過理論數(shù)據(jù)集試驗證明，基于布谷鳥搜索的聚類方法聚類精度要明顯優(yōu)于K均值的計算精度。

（3）通過應用實例表明，基于布谷鳥搜索的地震屬性聚類方法能夠更好地刻畫礁灘復合體、儲層、流體和斷裂的分布。