基于混合粒子群的核磁共振反演算法研究

宋 鵬 王國富

（桂林電子科技大學(xué)信息與通信學(xué)院，廣西 桂林 541000）

核磁共振找水反演是地球物理反演的分支，近年來出現(xiàn)了奇異值分解反演、模擬退火反演等方法研究。目前，粒子群優(yōu)化算法廣泛應(yīng)用于調(diào)度、模式識別、模糊系統(tǒng)控制、盲源分離等領(lǐng)域。在地球物理反演領(lǐng)域，主要應(yīng)用于波阻抗反演[1]和磁測資料反演[2]中，而在核磁共振反演領(lǐng)域的應(yīng)用很少。本文在研究核磁共振反演理論的基礎(chǔ)之上，應(yīng)用帶非線性約束優(yōu)化的混合粒子群算法，對核磁共振人工合成數(shù)據(jù)進(jìn)行數(shù)值反演仿真分析，開拓了粒子群優(yōu)化算法的新應(yīng)用領(lǐng)域，同時，為進(jìn)一步研究核磁共振反演方法提供理論依據(jù)。

1 核磁共振找水非線性反演理論

核磁共振找水反演，是將核磁共振數(shù)據(jù)資料解釋成地下水文地質(zhì)參數(shù)，本文主要解釋其中的地下含水量分布。它的主要方法可以分為線性反演方法和非線性反演方法兩大類。傳統(tǒng)的線性反演方法在對地下低電阻介質(zhì)情況下的數(shù)據(jù)進(jìn)行反演時，反演結(jié)果精度較差。非線性反演方法可以解決這一問題，非線性反演方法一般采用最優(yōu)化方法的理論求取地下含水量分布，即使用吉洪諾夫泛函[3]作為最優(yōu)化的目標(biāo)函數(shù)，目標(biāo)函數(shù)表達(dá)式如下所示：

并使目標(biāo)函數(shù)最小化。核磁共振找水非線性反演的

方程中，A為核函數(shù)矩陣，e0表示觀測到的初始振幅值，Nη即為需要求解的含水量值，是使吉洪諾夫泛函最小的解向量。

2 粒子群反演算法

粒子群反演算法[4]是地球物理非線性反演方法中的一種，其中，粒子群算法（Particle Swarm Optimzation，簡稱PSO)是Kennedy和Eberhart提出的基于仿生學(xué)原理發(fā)展起來的智能優(yōu)化算法。為求得目標(biāo)函數(shù)f(x)達(dá)到最優(yōu)解xp（最小值）的相應(yīng)x值，應(yīng)用粒子群算法的具體步驟為：

（1）給定粒子數(shù)目n，自變量x的個數(shù)m。隨機(jī)初始化所有粒子的位置xn×m和速度vn×m。

（2）存儲當(dāng)前各粒子的位置（記為）和相應(yīng)的適應(yīng)度值，找到的最小值，并記錄相應(yīng) f( xp)的位置， f( xp)記為xpg。

（3）按照公式（2）和公式（3）更新粒子的速度和位置。公式中，w為慣性權(quán)重，c1、c2為正的學(xué)習(xí)因子， r1和r2為0到1之間均勻分布的隨機(jī)數(shù)。p為當(dāng)前的個體極小值，pg,ji,j為當(dāng)前的全局極小值。

（4）更新后的適應(yīng)度值f(x)與步驟2中的f(xp)值進(jìn)行比較，較小值的相應(yīng)位置作為當(dāng)前的xp，當(dāng)前的 f(xp)再與f(xpg)進(jìn)行比較，較小值的相應(yīng)位置作為當(dāng)前的xpg。

（5）給定迭代次數(shù)，作為調(diào)節(jié)粒子速度和位置的次數(shù)，若未停止迭代，返回步驟（3）繼續(xù)搜索，若停止迭代，此時的xpg為最優(yōu)位置，f（xpg）為最優(yōu)解。

3 混合粒子群反演算法

為了避免基本粒子群算法陷入局部極小解，引入模擬退火算法[5]，該算法的概率受溫度參數(shù)的控制，具有概率突跳的能力，概率的大小隨溫度下降而減小?；旌狭Ｗ尤悍囱菟惴ǖ木唧w步驟如下：

（1）將種群中的各微粒的位置、速度進(jìn)行隨機(jī)初始化；

（2）對各個微粒的適應(yīng)度（相應(yīng)的函數(shù)值）進(jìn)行評價，存儲當(dāng)前各微粒的位置，設(shè)為pi，存儲全局最優(yōu)個體的位置，設(shè)為pg；

（3）給定初始溫度，根據(jù)公式（4）確定各pi的適配值：

（4）根據(jù)輪盤賭機(jī)制確定全局最優(yōu)值的相應(yīng)粒子位置pgbest，并根據(jù)公式（5）和公式（6）更新微粒的速度和位置：

（5）更新各微粒的局部pi值和全局pg值，并進(jìn)行退溫操作；

（6）若滿足預(yù)設(shè)的運(yùn)算精度或迭代次數(shù)，停止搜索，輸出最優(yōu)位置和最優(yōu)值，否則轉(zhuǎn)到步驟（3）；

為了將粒子群優(yōu)化算法適用于核磁共振找水反演中，設(shè)定非線性約束條件：作為粒子位置的可行域，各粒子在該可行域內(nèi)調(diào)節(jié)相應(yīng)的位置，搜尋最優(yōu)解。

4 測試函數(shù)數(shù)值計算

選擇 3個常用的測試函數(shù)[6]來驗證混合粒子群算法的有效性。

函數(shù)一：Sphere函數(shù)，函數(shù)表達(dá)式為：

函數(shù)一在 xi=0處達(dá)到全局極小點(diǎn)，本文計算結(jié)果為0.000029465。函數(shù)一的迭代過程如圖1所示。

圖 1 sphere函數(shù)迭代過程

函數(shù)二：Rosenbrock函數(shù),函數(shù)表達(dá)式為：

函數(shù)二在 xi=1時達(dá)到全局極小點(diǎn)，本文利用混合粒子群算法的計算結(jié)果為0.9984。函數(shù)二的迭代過程如圖2所示。

圖2 Rosenbrock函數(shù)迭代過程

函數(shù)三：Rastrigin函數(shù)，函數(shù)表達(dá)式為：

函數(shù)三是多峰函數(shù)，在 xi=0處達(dá)到全局極小點(diǎn)，本文計算結(jié)果為0.0000013189。函數(shù)三的迭代過程如圖3所示。

圖3

5 基于混合粒子群反演實驗分析

在均勻半空間條件下，假設(shè)大回線源線圈的半徑為50m，地下電阻率50Ω，地磁傾角60°。對目標(biāo)函數(shù)應(yīng)用粒子群優(yōu)化算法進(jìn)行最優(yōu)化求解，正演擬合初始振幅如圖（a）所示。反演擬合含水量如圖（b)所示。

圖4 粒子群反演擬合結(jié)果

淺水層的反演擬合效果較好，深水層反演擬合效果較差，但是整體的垂直分辨率較高。由此可見，混合粒子群反演算法適合于核磁共振找水反演中。

[1] 易遠(yuǎn)元,王家映.粒子群反演方法[J].工程地球物理學(xué)報,2009,6(4):385-389.

[2] 張大蓮,劉天佑,陳石羨,等.粒子群算法在磁測資料井地聯(lián)合反演中的應(yīng)用[J].物探與化探,2009,33(5):571-576.

[3] 林君,段清明,王應(yīng)吉,等.核磁共振找水儀原理與應(yīng)用[M].北京:科學(xué)出版社,2010.

[4] 方偉,孫俊,須文波.一種多樣性控制的粒子群優(yōu)化算法[J].控制與決策,2008,23(8):863-868.

[5] 王華秋,曹長修.基于模擬退火的并行粒子群優(yōu)化研究[J].控制與決策,2005,20(5):500-504.

[6] 熊勇.粒子群優(yōu)化算法的行為分析與應(yīng)用實例[D].杭州:浙江大學(xué),2005.