基于改進(jìn)的量子粒子群算法的變差函數(shù)擬合方法及應(yīng)用

謝 源, 李 瓊, 陳 杰, 何榮勝

(成都理工大學(xué) 地球物理學(xué)院，成都 610059)

在進(jìn)行儲層平面非均質(zhì)性研究中，變差函數(shù)是最常見的研究方法。通過計(jì)算實(shí)驗(yàn)變差函數(shù)，然后對其進(jìn)行擬合得到相關(guān)的參數(shù)，再對其進(jìn)行地質(zhì)解釋。然而在實(shí)驗(yàn)變差函數(shù)擬合的過程中，由于理論變差函數(shù)參數(shù)眾多且通常為非連續(xù)可導(dǎo)的，所以對其擬合一直是一個難點(diǎn)。最初對于變差函數(shù)的擬合都采用人工的方法[1]，但是這種方法得到的變差函數(shù)參數(shù)沒有一個統(tǒng)一的評價(jià)標(biāo)準(zhǔn)，擬合的結(jié)果往往因人而異，很難得到最佳的擬合效果。當(dāng)前，變差函數(shù)擬合的方法主要有最小二乘法[2]、極大似然法[3]、線性規(guī)劃法[4]、加權(quán)多項(xiàng)式回歸法[5]和遺傳算法[6]等，但是各種方法都有其相應(yīng)的局限性。最小二乘方法雖然簡單易行，但是由于無法確定待擬合參數(shù)的正負(fù)號，故無法得到最優(yōu)解；而極大似然法則計(jì)算復(fù)雜，而且有時不能得到最優(yōu)解；加權(quán)多項(xiàng)式以及線性規(guī)劃法依然不能解決參數(shù)正負(fù)號的問題；遺傳算法雖然在對于求解非線性優(yōu)化時具有全局尋優(yōu)的特點(diǎn)，但是由于其可調(diào)參數(shù)過多且結(jié)構(gòu)復(fù)雜從而影響了效率。

J.Kennedy教授和R.C.Eberhar教授在1995年提出了粒子群算法(Particle Swarm Optimization, 簡稱為PSO)[7]。它是一種基于群智能的全局尋優(yōu)算法。其基本原理是模擬鳥群在覓食的過程中其個體間共享所了解的食物源信息，從而在準(zhǔn)確的信息下使整個群體向著食物源的方向飛行。由于PSO算法收斂速度快、設(shè)置參數(shù)少，而且不受函數(shù)是否連續(xù)可導(dǎo)的條件的限制，使其在非線性擬合方面得到了廣泛的應(yīng)用。但是傳統(tǒng)的PSO算法具有易陷入局部極值的缺陷，有時不能很好地解決問題。因此，本文提出了一種改進(jìn)的PSO算法，并利用改進(jìn)的算法實(shí)現(xiàn)了對實(shí)驗(yàn)變差函數(shù)進(jìn)行自動擬合。

1 變差函數(shù)原理

1.1 變差函數(shù)的定義

變差函數(shù)是J.Motheron等[8]提出的一種矩估計(jì)方法。其表示為區(qū)域化變量的增量平方的數(shù)學(xué)期望，也就是區(qū)域化變量的增量的方差。根據(jù)地質(zhì)統(tǒng)計(jì)學(xué)理論,通過半變差函數(shù)，利用適當(dāng)?shù)臄?shù)理模型(如球狀模型、高斯模型、指數(shù)模型等)，可對由自然地質(zhì)規(guī)律形成的儲層參數(shù)的空間相關(guān)性進(jìn)行定量描述。下式為實(shí)驗(yàn)變差函數(shù)的計(jì)算公式

(1)

式中：N(h)表示距離為h的點(diǎn)個數(shù)；Xi為第i個點(diǎn)坐標(biāo)；Z(Xi)、Z(Xi+h)分別為Xi和Xi+h兩點(diǎn)處的觀測值；h為兩觀測點(diǎn)間的距離；r*(h)為實(shí)驗(yàn)變差函數(shù)值。

1.2 模型的選擇

變差函數(shù)的模型有多種，常用的模型主要有球狀模型、指數(shù)模型、高斯模型、冪函數(shù)模型、對數(shù)模型等。各種模型的儲層參數(shù)在同一方向上的變化速度不同，變程的表達(dá)形式也不同，變程的大小能反映儲層參數(shù)沿某一方向的變化速度的快慢。基臺值是反映儲層參數(shù)在某一方向變化的幅度的大小，基臺值越大，參數(shù)變化幅度越大，非均質(zhì)性越強(qiáng)；反之，基臺值越小，參數(shù)變化幅度也越小，非均質(zhì)性越弱[9]。球狀模型能夠更快速地達(dá)到基臺值，當(dāng)達(dá)到基臺值后，變程的增加不再會影響函數(shù)值，因此實(shí)際應(yīng)用中多為球狀模型。球狀模型的表達(dá)式

(2)

式中：C0為塊金常數(shù)；C為拱高；C0+C為基臺值；a為變程。

2 PSO算法及其改進(jìn)

2.1 標(biāo)準(zhǔn)粒子群算法

在PSO算法中，首先假設(shè)在D維搜索空間中，初始化M個粒子，將Xi=(Xi1,Xi2, …,XiD)記為第i個粒子的所在位置，可以根據(jù)目標(biāo)函數(shù)計(jì)算得到Xi的適應(yīng)度值，從而通過適應(yīng)度值來判斷其位置的優(yōu)劣；將Vi=(Vi1,Vi2, …,ViD)記為第i個粒子的飛行速度；將pbi=(pbi1,pbi2, …,pbiD)記為第i個粒子所經(jīng)過的位置中最優(yōu)的位置；將gbi=(gbi1,gbi2, …,gbiD)定義當(dāng)前群體經(jīng)過的最優(yōu)位置。然后粒子通過下式進(jìn)行速度和位置的迭代更新，其更新公式如下

(3)

(4)

其中:d(1, …,D)是空間的維度；i(1, …,m)為種群的粒子數(shù)；R1、R2是介于(0,1)區(qū)間的隨機(jī)數(shù)；k(1, …,n)是迭代次數(shù)；c1和c2為學(xué)習(xí)因子，其可以控制種群中個體向全體中其他優(yōu)秀個體學(xué)習(xí)的能力，通常取c1與c2的值為2。式(3)中，第一項(xiàng)是粒子當(dāng)前速度；第二項(xiàng)是粒子個體所經(jīng)歷的最好位置，這可以看作是粒子個體對自身經(jīng)驗(yàn)的總結(jié)；第三項(xiàng)代表了種群的當(dāng)前最優(yōu)位置，可以看作是粒子間信息共享的體現(xiàn)。

2.2 改進(jìn)粒子群算法

PSO算法在進(jìn)行函數(shù)擬合時雖然簡單易行，但是標(biāo)準(zhǔn)的PSO算法不能很好地提高算法的收斂性。其主要原因在于，PSO算法中粒子的收斂形式是以軌道的方式實(shí)現(xiàn)的；而且它在搜索過程中受最大速度的限制，使得PSO算法不是以概率1進(jìn)行全局收斂。基于此，Sun等[10]提出了量子粒子群算法(Quantum-behaved Particle Swarm Optimization)。QPSO在一維空間中引入δ勢阱，假設(shè)粒子在以P點(diǎn)為中心的δ勢阱中，在粒子收斂過程中其速度和位置在量子空間中不能一起確定，其狀態(tài)可以用波函數(shù)φ(Y)來表示，由式(5)確定,其中L是δ勢阱的特征長度。

(5)

使用蒙特卡洛方法對粒子進(jìn)行模擬，得到粒子的位置方程如下

(6)

其中:X為粒子的位置；P為吸引勢的中心點(diǎn)；u是(0, 1)上服從均勻分布的隨機(jī)數(shù)。由于L為時間t的函數(shù)，則量子粒子群算法中粒子的基本進(jìn)化方程如下

(7)

然后我們將其運(yùn)用于n維空間中,則吸引子Pi可以表示為Pi=(Pi1, …,Pin)，在每一維中用Pn為坐標(biāo)，以Pi,j為中心建立一位δ勢阱，則粒子i的波函數(shù)表示為

φ[xi,j(t+1)]=

(8)

位置方程

(9)

定義平均最好位置，其計(jì)算公式

(10)

則有

Li,j(t)=2α|mbestj(t)-xij(t)|

(11)

于是得到粒子的位置更新方程

(12)

其中:ui,j(t)～u(0,1);α為收縮擴(kuò)張因子。

綜上可以得出量子粒子群的更新方程如式(13)～(15)。

(13)

(14)

(15)

雖然由均勻分布產(chǎn)生的粒子具有較好的全局搜索能力，但是其局部搜索能力較弱。而柯西分布很容易得到距離原點(diǎn)較遠(yuǎn)的隨機(jī)數(shù)，因此，采用柯西分布產(chǎn)生隨機(jī)數(shù)可以使得算法很快跳出局部極值。標(biāo)準(zhǔn)QPSO算法中式(14)～(15)中參數(shù)φ1、φ2、u的值是由均勻分布產(chǎn)生的隨機(jī)數(shù)來得到的。在改進(jìn)算法中，我們采用柯西分布產(chǎn)生的隨機(jī)數(shù)來確定φ1、φ2、u的值，并且采用線性變化的收縮擴(kuò)張因子來對QPSO進(jìn)行改進(jìn)。收縮擴(kuò)張因子α的取值如下

(16)

其中:tmax為最大迭代次數(shù);t為當(dāng)前迭代次數(shù)。

2.3 算法的實(shí)現(xiàn)流程

在變差函數(shù)擬合時，將所求參數(shù)一組解看作是一個粒子，將其個數(shù)看作是粒子的維度。對于一階球狀模型可以得到算法中粒子i的適應(yīng)度函數(shù)如式(17)。

一階球狀模型

F(i)=

(17)

其中:hj為第j個滯后距；r*(h)為計(jì)算得到的實(shí)驗(yàn)變差函數(shù)值。

改進(jìn)QPSO算法具體的擬合流程如圖1。

3 仿真測試與實(shí)例計(jì)算

3.1 仿真測試

為了測試所改進(jìn)的量子粒子群算法的性能，我們進(jìn)行了一個簡單的仿真測試。選擇的測試函數(shù)及其性質(zhì)如表1所示。

采用上述3個函數(shù)來對所改進(jìn)的量子粒子群算法進(jìn)行測試：模型一中設(shè)定粒子的數(shù)量為20、維度為10、迭代次數(shù)為 2 000次，分別用標(biāo)準(zhǔn)QPSO算法和改進(jìn)量子粒子群算法運(yùn)行30次并計(jì)算最優(yōu)值的平均值；模型二中設(shè)定粒子數(shù)量為60、維度為20，迭代次數(shù)不變，運(yùn)行30次并計(jì)算最優(yōu)值的平均值。2種模型得到的最優(yōu)值的平均值見表2。

圖1 改進(jìn)量子粒子群擬合流程Fig.1 Fitting process of the improved quantum particle swarm

函數(shù)名稱表達(dá)式自變量范圍最小值Sphere函數(shù)f(x)=∑Di=1x2i(-100, 100)0Rosenbrock函數(shù)f(x)=∑Di=1(100(xi+1-x2i)2+(xi-1)2)(-100, 100)0Rastrigrin函數(shù)f(x)=∑Di=1(x2i-10cos(2πxi)+10(-5.12, 5.12)0

從測試的結(jié)果我們可以看出，改進(jìn)的量子粒子群算法相比較于標(biāo)準(zhǔn)的QPSO算法所求得的函數(shù)極值更加接近于0。其在求解函數(shù)的極值的精度上有明顯的提高，特別是對高維函數(shù)其優(yōu)越性更加明顯。

表2 模型測試結(jié)果Table 2 Model test results

3.2 實(shí)例計(jì)算

對塔河油田某井的聲波時差數(shù)據(jù)按照式(2)計(jì)算實(shí)驗(yàn)變差函數(shù)值。由于所選層段的非均質(zhì)性較強(qiáng)，在實(shí)際計(jì)算過程中發(fā)現(xiàn)，當(dāng)滯后距為30 m時計(jì)算得到的實(shí)驗(yàn)變差函數(shù)的穩(wěn)健性稍好。其計(jì)算結(jié)果見表3。

首先使用傳統(tǒng)的加權(quán)最小二乘擬合對得到的實(shí)驗(yàn)變差函數(shù)進(jìn)行擬合，擬合模型選擇球狀模型。得到擬合參數(shù)值為：C0=4.810 9,C=7.923 3,a=294.533 9。

在使用改進(jìn)的QPSO算法擬合時，我們選擇一階球狀模型，粒子的群體規(guī)模為30，維度設(shè)置為3，迭代次數(shù)設(shè)置為 1 500次，其中參數(shù)的取值范圍分別設(shè)置為：C0=[0,15],C=[0,15],a=[0,500]。獨(dú)立運(yùn)行程序10次，得到的待擬合的參數(shù)取值為：C0=4.325 4,C=7.342 4,a=233.416 5，此時得到F(i) =37.75。

然后將2種方法得到的參數(shù)代入球狀模型得到理論變差函數(shù)值(表3)。

表3 塔河油田某井的實(shí)驗(yàn)變差函數(shù)及理論變差函數(shù)值Table 3 Experimental variogram and theoretical variogram value from a drilling well in Tahe area

h為滯后距；r(h)為計(jì)算得到的實(shí)驗(yàn)變差函數(shù)值；r*(h)1為改進(jìn)的量子粒子群算法所得到的理論變差函數(shù)值；r*(h)2為最小二乘算法得到的理論變差函數(shù)值

2種擬合方法得到的變差函數(shù)如圖2，從圖中可以看出，藍(lán)色的折線為計(jì)算所得的實(shí)驗(yàn)變差函數(shù)值，紅線為利用改進(jìn)的量子粒子群算法擬合得到的理論變差函數(shù)曲線，綠線為用加權(quán)最小二乘擬合得到的理論變差函數(shù)曲線。從中我們可以看出利用改進(jìn)的量子粒子群算法擬合得到的理論變差函數(shù)曲線比加權(quán)最小二乘得到的曲線具有更小的變程，即其非均質(zhì)性更強(qiáng)，擬合結(jié)果也更加貼合計(jì)算得到的實(shí)驗(yàn)變差函數(shù)曲線，所以本文提出的改進(jìn)量子粒子群算法對變差函數(shù)的擬合取得了比較理想的效果。

圖2 實(shí)驗(yàn)變差函數(shù)擬合結(jié)果Fig.2 Results of the experimental variogram fitting

我們引入擬合度和絕對偏差來對2種擬合方法的結(jié)果進(jìn)行評價(jià)，其計(jì)算公式

(18)

(19)

分別計(jì)算改進(jìn)的量子粒子群算法和加權(quán)最小二乘擬合得到的理論變差函數(shù)值的擬合度和絕對偏差見圖3。從圖中可以看出，采用改進(jìn)的量子粒子群算法擬合得到的理論變差函數(shù)的擬合度和絕對偏差都是小于由加權(quán)最小二乘方法所得到的理論變差函數(shù)的。因此，本文的方法可以用于實(shí)驗(yàn)變差函數(shù)的擬合。

3.3 儲層平面非均質(zhì)性預(yù)測

對某儲層的孔隙度數(shù)據(jù)在4個方向(0°, 45°, 90°, 135°)上分別計(jì)算實(shí)驗(yàn)變差函數(shù)，然后采用改進(jìn)QPSO算法對實(shí)驗(yàn)變差函數(shù)進(jìn)行擬合。在擬合時，設(shè)定粒子數(shù)量為50，維度為3，迭代次數(shù)為2 000次，得到的擬合效果如圖4。

由圖4可以看出，在4個方向上該儲層的孔隙度變差函數(shù)具有較明顯的差異，表明儲層存在非均質(zhì)性。圖中可以看出，在0°和90°方向，變程a的值明顯小于45°和135°方向的變程值，這表明在這2個方向上該儲層具有較強(qiáng)的非均質(zhì)性。在45°和135°方向，變程a值較大，表明該方向上非均質(zhì)性較弱。

4 結(jié) 論

a.本文采用柯西分布產(chǎn)生隨機(jī)數(shù)并采用線性變化的收縮擴(kuò)張因子的方法來改進(jìn)QPSO算法，對改進(jìn)算法進(jìn)行了仿真，結(jié)果表明，改進(jìn)的QPSO算法比標(biāo)準(zhǔn)的QPSO算法在求解函數(shù)的極值上具有更高的精度。

b.運(yùn)用改進(jìn)的量子粒子群算法對實(shí)驗(yàn)變差函數(shù)進(jìn)行擬合，相比較于傳統(tǒng)的加權(quán)最小二乘擬合變差函數(shù)的方法，改進(jìn)算法得到的結(jié)果的擬合度和絕對偏差更小，擬合的結(jié)果更接近于實(shí)驗(yàn)計(jì)算值，是一種更優(yōu)的變差函數(shù)擬合方法。

圖3 兩種方法的擬合度和絕對偏差對比Fig.3 Contrast of the fitting and absolute deviations of the two methods

圖4 某儲層4個方向孔隙度變差函數(shù)Fig.4 Porosity variation function in four directions of a reservoir

c.運(yùn)用改進(jìn)QPSO方法對塔河油田某儲層的平面非均質(zhì)性進(jìn)行研究，對儲層孔隙度實(shí)驗(yàn)變差函數(shù)進(jìn)行擬合，從擬合結(jié)果看出該儲層在0°和90°方向具有較強(qiáng)的非均質(zhì)性，在45°和135°方向非均質(zhì)性較弱。

[參考文獻(xiàn)]

[1] 王仁鐸，胡光道.線性地質(zhì)統(tǒng)計(jì)學(xué)[M].北京：地質(zhì)出版社，1998.

Wang R D, Hu G D. Linear Geostatistics[M]. Beijing: Geological Publishing House, 1998. (in Chinese)

[2] 李正文，李瓊，吳朝容.沉積盆地有效儲集層綜合識別技術(shù)[M]．成都：四川科學(xué)技術(shù)出版社，2002.

Li Z W, Li Q, Wu C R. Comprehensive Identification Technology of Effective Reservoirs in Sedimentary Basins[M]. Chengdu: Sichuan Science and Technology Press, 2002. (in Chinese)

[3] Cressie N. Fitting variogram models by weighted least squares[J]. Journal of the International Association for Mathematical Geology, 1985, 17(5): 563-586.

[4] Mardia K V, Marshall R J. Maximum likelihood parameter estimation of models for residual covariance in spatial regression[J]. Biometrika, 1984, 71(1): 135-138.

[5] 矯希國，劉超.變差函數(shù)的參數(shù)模擬[J].物探化探計(jì)算技術(shù)，1996，16(2)：157-161．

Jiao X G, Liu C. Estimation of variation parameter[J]. Computing for Geophysical and Geochemical Exploration, 1996, 16(2): 157-161. (in Chinese)

[6] 王仁鐸.用加權(quán)多項(xiàng)式回歸進(jìn)行球狀模型變差圖的最優(yōu)擬合[J].地球科學(xué)，1986，11(2)：143-151.

Wang R D. On the best fitting of the spherical model variogram with the weighted polynomial regression[J]. Earth Science, 1986, 11(2): 143-151. (in Chinese)

[7] 馮國慶，周涌沂，李允，等.利用遺傳算法擬合實(shí)驗(yàn)變差函數(shù)[J].西南石油學(xué)院學(xué)報(bào)，2005，27(5)：23-25．

Feng G Q, Zhou Y Y, Li Y,etal. Application of genentic algorithm to match theoretical variogram[J]. Journal of Southwest Petroleum Institute, 2005, 27(5): 23-25. (in Chinese)

[8] Kennedy J, Eberhart R C. Particle swarm optimization[C]//Proceeding of IEEE International Conference on Neural Networks. Piscataway: IEEE Press, 1995: 1942-1948．

[9] Matheron G. Principles of geostatistics[J]. Economic Geology, 1963, 58(8): 1246-1266.

[10]Sun J, Xu W, Feng B. Adaptive parameter control for quantum-behaved particle swarm optimization on individual level[C]//2005 IEEE International Conference on Systems, Man and Cybernetics. IEEE, 2005, 4: 3049-3054.

敬告作者

為適應(yīng)我國科技信息化建設(shè)需要，擴(kuò)大作者學(xué)術(shù)交流渠道，本刊已加入《中國學(xué)術(shù)期刊(光盤版)》和《中國知網(wǎng)》(http://www.cnki.net)、萬方數(shù)據(jù)電子出版社的《萬方數(shù)據(jù)——數(shù)字化期刊群》(http://www.wanfangdata.com.cn)、教育部科技發(fā)展中心的《中國科技論文在線》、重慶維普資訊有限公司的《中文科技期刊數(shù)據(jù)庫》、華藝數(shù)位藝術(shù)股份有限公司的《CEPS中文電子期刊》、北京書生網(wǎng)絡(luò)技術(shù)有限公司的《書生數(shù)字期刊》、北京世紀(jì)超星公司的“域出版”平臺、中教數(shù)據(jù)庫等。作者著作權(quán)使用費(fèi)與稿酬由本刊一次性給付。如果作者不同意將文章編入上述數(shù)據(jù)庫，請?jiān)趤砀鍟r聲明，本刊將做適當(dāng)處理。

《成都理工大學(xué)學(xué)報(bào)(自然科學(xué)版)》編輯部