曹 爽，賀玉珍，安建成

(1.運城學院 數(shù)學與信息技術學院，山西 運城 044000；2.太原理工大學 信息與技術學院，山西 晉中 030600)

0 引 言

一維OTSU閾值法是由日本學者大津于1979年提出的[1]，其求解閾值簡單、分割速度較快、提取效果良好，因此得到了廣泛應用。但一維OTSU法只使用圖像的灰度信息，沒有涉及空間鄰域信息[2]，當圖像有噪聲干擾時無法得到滿意的分割結果。為提高抗噪性和分割準確率，研究者們將OTSU算法拓展到二維和三維[3]。其中三維OTSU閾值法引入了鄰域中值及均值信息，可得到理想的分割效果。

由于中值濾波易濾除椒鹽噪聲，均值濾波易濾除高斯噪聲[1]，因此，三維OTSU法既可對椒鹽噪聲或高斯噪聲有較高的抗噪性，又可對疊加混合噪聲(椒鹽噪聲加高斯噪聲)的圖像有良好的提取結果。然而OTSU法采用窮盡搜索策略，對于三維OTSU法的計算需要三重循環(huán)，導致計算量過多，求解復雜度呈指數(shù)型增長，實時性較差。為此，景曉軍等[4]提出了三維OTSU快速遞推法，使計算效率提高，但分割準確率一般，在此基礎上，范九倫等[5]對遞推公式進行修正，給出新的遞推算法，使求解精度提高。葉志偉等[6]對螢火蟲算法(firely algorithm，F(xiàn)A)進行研究，將萊維飛行引入其中，用改進的螢火蟲算法來搜尋最優(yōu)分割閾值，使分割準確率顯著上升。仇國慶等[7]提出了改進的三維OTSU法，使用一維OTSU法縮小搜索空間以及迭代空間，并采用布谷鳥算法(cuckoo search，CS)求解閾值，既降低了計算量又具有較高的抗噪性。

新的群智能算法—狼群算法(wolf pack algorithm，WPA)是由吳虎勝等[8]于2013年提出的，其參數(shù)設置簡單，已被大量應用在多個領域。對狼群算法進行研究，將花授粉算法(flower pollination algorithm，F(xiàn)PA)[9]求解的最優(yōu)閾值加入到狼群算法的初始閾值，再在圍攻過程中引入高斯變異，將改進的狼群算法與三維OTSU算法結合，可取得較好的分割結果。

1 三維OTSU閾值法

三維OTSU算法不僅考慮了像素點灰度值，還考慮了像素點鄰域均值和中值。設計包含像素灰度及其鄰域均值和中值的三維直方圖，根據(jù)三維直方圖求取分割閾值，從而實現(xiàn)自動閾值選取。

對于m×n的圖像，用h(x,y)表示點(x,y)的像素灰度值。對于某一個選定點k×k區(qū)域內(nèi)的平均灰度值用f(x,y)表示，計算公式如下：

(1)

該選定點k×k區(qū)域內(nèi)的灰度中值是g(x,y)，計算公式如下：

g(x,y)=med{h(x+m,y+n)}

(2)

其中，m=-k/2,…,k/2;n=-k/2,…,k/2。

根據(jù)上面的一系列推導可知，若L是圖像的灰度等級，則f(x,y)和g(x,y)的灰度等級也都為L。將h(x,y)、f(x,y)以及g(x,y)組成的三元組(i,j,k)稱作三維直方圖，如圖1(a)所示。對于直方圖上某一點代表的向量(i,j,k)，其出現(xiàn)的概率是pi,j,k，計算公式為：

(3)

對于一個可行的三維閾值向量(s,t,q)可將三維直方圖分成八個部分，如圖1(b)～(d)所示。其中，將區(qū)域1和0分別劃分成圖像的目標和背景，區(qū)域2～7可看成噪聲和邊緣。由于目標和背景區(qū)域的像素點多于邊緣和噪聲區(qū)域的像素點，從而將區(qū)域2～7的像素點忽略不計。

圖1 三維OTSU直方圖

令C1和C0分別代表目標和背景這兩個區(qū)域，則目標與背景發(fā)生的概率如下：

(4)

(5)

目標與背景均值矢量分別為μ1和μ0：

(7)

直方圖總均值矢量μT為：

μT=(μTi,μTj,μTk)T=

(8)

由于區(qū)域2～7的像素點幾乎為0，則：

p0+p1=1

(9)

p0μ0+p1μ1=μT

(10)

那么目標與背景區(qū)域之間的離散度矩SB為：

SB=p0[(μ0-μt)(μ0-μt)T]+

p1[(μ1-μt)(μ1-μt)T]

(11)

將離散度矩陣跡的值TrSB作為兩類區(qū)域的類間距度量函數(shù)：

TrSB=p0[(μ0i-μTi)2+(μ0j-μTj)2+(μ0k-μTk)2]+p1[(μ1i-μTi)2+(μ1j-μTj)2+(μ1k-μTk)2]

(12)

當Tr(SB(s*,t*,q*))=max{Tr(SB(s,t,q))},0≤s,t,q≤L-1時，對應的(s*,t*,q*)為圖像最佳分割閾值。

2 基于改進狼群算法的三維OTSU法

2.1 狼群算法原理

狼群算法是通過個體狼與整個狼群進行信息交互來搜索全局最優(yōu)解[10-11]。一個狼群體系包括頭狼、猛狼以及探狼這三種成員。頭狼負責指揮其他個體狼搜索獵物。探狼隨機游走去搜尋獵物，根據(jù)獵物氣味濃度判斷獵物源的遠近，濃度越大(離散度矩陣跡的值越大)則狼離獵物越近。探狼發(fā)現(xiàn)獵物會馬上通知頭狼，頭狼向猛狼發(fā)起召喚，猛狼收到指令進行圍攻追捕獵物。為保證種群的生存和發(fā)展，狼群會不斷更新群體，除去瘦弱的個體狼，保存種群實力。

2.2 改進的狼群算法

2.2.1 花授粉算法優(yōu)化狼群算法初值

花授粉算法是由Yang Xinshe于2012年提出的。由于該算法參數(shù)少、結構簡單、全局搜索性能良好，將其引入到狼群算法中，提高狼群算法的全局尋優(yōu)能力。

在花授粉算法初始時，隨機產(chǎn)生M個種群個體p(t)={xit}，其中xit=[xti,1,xti,2,…,xti,j,…,xti,Q]，i=1,2,…,M;Q為優(yōu)化問題的維度，j=1,2,…,Q;t是迭代次數(shù)?；ㄊ诜鬯惴ɡ棉D換概率p∈[0,1]來把控自花授粉和交叉授粉之間的轉換。

自花授粉為局部授粉，其過程為：

x(t+1)i=x(t)i+ε*(x(t)j-x(t)i)

(13)

其中，i,j為[1,M]區(qū)間內(nèi)的隨機整數(shù)，ε為變異因子，是[0,1]上服從均勻分布的隨機數(shù)。

交叉授粉為傳粉者采用萊維飛行的全局授粉，其過程為：

x(t+1)i=x(t)i+γL(x(t)i-g*)

(14)

其中，γ為步長控制因子，g*為當前群體最佳值，參數(shù)L是步長，服從萊維分布。

轉換概率p的計算公式為：

p=0.8+0.2×rand

(15)

其中，rand為[-1,1]之間的隨機數(shù)。

花授粉算法主要步驟如下：

(1)初始化參數(shù)。設置群體數(shù)量M和種群個體p(t)，優(yōu)化問題維度Q，迭代次數(shù)t，變異因子ε，步長控制因子γ。

(2)計算個體適應度值，確定全局最佳適應度值和對應的最佳位置g*。

(3)根據(jù)式(15)計算轉換概率p。

(4)在[0,1]區(qū)間生成一個隨機數(shù)rand1，若p≤rand1，則根據(jù)式(13)進行局部授粉。

(5)若p> rand1，則根據(jù)式(14)進行全局授粉。

(6)判斷算法是否滿足循環(huán)結束的條件，若滿足則輸出結果；否則，返回步驟(2)。

狼群算法初始時位置是隨機的，若位置分布得均勻，則容易尋找全局最佳值；若位置分布得不均勻，則容易陷入局部極值，降低尋優(yōu)準確率。為了提高尋優(yōu)精度，加快全局搜索速度，將花授粉算法求解的較優(yōu)值作為狼群算法的初始值。

具體實現(xiàn)過程如下：令個體狼的總數(shù)為N，待求解的閾值變量數(shù)為D，那么狼群被抽象在N×D的歐氏空間中。設待尋優(yōu)的初始閾值向量為Xi=(xi1,xi2,…,xiD)，則Yi=f(Xi)(Y是目標函數(shù)，即離散度矩陣跡的值)為每條人工狼感知的獵物氣味濃度。將花授粉算法求得的最佳解對應的位置加入到狼群算法的初始位置?；ㄊ诜鬯惴ㄖ校衾侨核惴ǖ娜后w數(shù)量N大于M，可將花授粉算法求解的最佳位置全部隨機加入狼群算法的初始位置；否則，只將花授粉算法求解的個別最佳位置隨機加入狼群算法的初始位置。

待初始位置確定后，選出頭狼(具有最優(yōu)目標函數(shù)的人工狼)記作Ylead。頭狼不參加游走行為、召喚行為和圍攻行為，它直接進入下一輪迭代，直到其被更兇猛的人工狼取代。

2.2.2 游走過程

先選取L_Num只人工狼作為探狼，L_Num取[n/(α+1),n/α]中的整數(shù)(α是探狼比因子)。探狼群體從各自的當前位置隨機向m個方向游走stepa(游走步長)，同時記錄游走后的獵物氣味濃度(目標函數(shù)值)并回到初始位置。向第h(h=1,2,…,m)個方向游走后，探狼i在第d維求解空間中的位置公式如下:

xhid=xid+sin(2π×h/m)*stepha

(16)

Yih(h=1,2,…,m)為探狼i感知的獵物氣味濃度。探狼個體在h個方向中選擇氣味濃度最強的方向游走，對探狼的位置Xi和感知的氣味濃度Yi=max{Yih}進行更新。探狼群體重復上述游走行為，當某只個體狼i感知的氣味濃度Yi>Ylead時，將這只人工狼作為頭狼發(fā)起召喚，否則重復游走行為直到T(游走次數(shù))達到上限。

2.2.3 召喚過程

將余下的(N-1-L_Num)只人工狼作為猛狼。猛狼聽到頭狼的召喚后迅速追捕stepb(奔襲步長)逼近獵物源，那么對于第k+1次迭代第i只猛狼在第d維求解空間中的位置公式如下:

xk+1id=xkid+stepdb*(gkd-xkid)/|gkd-xkid|

(17)

其中，gkd為第k代頭狼在第d維求解空間的位置。追捕過程中，如果猛狼嗅到的氣味濃度Yi>Ylead，那么Ylead=Yi，將該猛狼作為頭狼進行召喚，反之繼續(xù)追捕直至其與頭狼Q之間的距離diq小于dnr時進入圍攻行為。dnr的計算公式如下:

(18)

其中，xmaxd與xmind分別是第d維求解空間的上下限，w是對距離判定的因子。

2.2.4 改進的圍攻過程

令Gkd為第k次迭代，獵物在第d維求解空間里的位置，那么圍攻過程的公式如下:

xk+1id=xkid+λ*stepdc*|Gkd-xkid|

(19)

其中，stepc為圍攻步長，λ為隨機步長因子([-1,1]中均勻分布的隨機數(shù))。如果探狼和猛狼嗅到的獵物氣味濃度超過它們原位置嗅到的氣味濃度，則對人工狼位置進行更新，反之位置不變。stepa、stepb以及stepc之間關系如下:

stepda=stepdb/2=2*stepdc

(20)

stepda=|xmaxd-xmind|/s

(21)

原始圍攻過程中，狼群搜索速度逐漸趨于0，使得算法出現(xiàn)“假死”狀態(tài)，降低算法搜索效率。為此，將高斯變異引入圍攻行為，使得算法后期能夠從局部最佳值附近深度搜索全局最佳值，加快尋優(yōu)速度，并減少狼群陷入局部最佳的次數(shù)。

根據(jù)高斯分布概率密度函數(shù)構造高斯變異函數(shù)：

(22)

其中，σ為高斯分布的標準差，x*為全局最佳值。

通過GMFi(t)的擾動，可使狼群在圍攻行為中跳出局部極值。根據(jù)高斯變異移動的位置公式如下：

x(t)i=GMF(t)i*x(t)i

(23)

設置變異函數(shù)是為了讓所有個體狼都有機會發(fā)生變異移動。變異函數(shù)的公式如下：

(24)

其中，N為種群數(shù)量，n為函數(shù)維度。當0

2.2.5 群體更新

整個狼群根據(jù)“強多弱少”的原則分配獵物。將算法中目標函數(shù)值最差的L匹人工狼除去，并隨機生成L匹狼，L取[n/(2×ξ),n/ξ](ζ為狼群更新因子)內(nèi)的隨機數(shù)。

2.3 改進狼群算法優(yōu)化的三維OTSU法

2.3.1 改進算法原理

利用改進的狼群算法優(yōu)化三維OTSU算法的本質(zhì)即人工狼個體表示一個可行的三維閾值向量[12]。在像素點灰度以及該像素點鄰域灰度均值和中值構成的三維向量空間中[13-14]，狼群通過游走、召喚、圍攻、種群更新這幾種行為的不斷迭代以及狼群間的信息交互來求取一組最優(yōu)解(s*,t*,q*)，使得式(12)的值最大即以目標函數(shù)最大化為基礎來尋找最優(yōu)閾值。

2.3.2 改進算法步驟

(1)利用花授粉算法求取最佳解和對應的最佳位置。

(2)設置狼群算法參數(shù)。確定人工狼位置Xi以及人工狼的數(shù)目N，迭代次數(shù)上限kmax和游走次數(shù)上限Tmax，探狼比因子α，步長因子s，距離判定因子w以及種群更新因子ζ。

(3)種群初始化。將花授粉算法得到的最佳位置加入狼群算法的初始位置。

(4)人工狼的位置為初始化的二維閾值向量，通過式(12)計算人工狼嗅到的獵物氣味濃度。根據(jù)氣味濃度大小選出頭狼、探狼以及猛狼。

(5)若某只探狼嗅到的氣味濃度Yi>Ylead，將這只探狼作為頭狼或達到游走次數(shù)上限Tmax，轉步驟(6)。

(6)猛狼聽到頭狼召喚追捕獵物。如果其嗅到的氣味濃度Yi大于頭狼Ylead，則Ylead=Yi，該猛狼作為頭狼召喚狼群。如果Yi≤Ylead，那么猛狼繼續(xù)追捕直至dis≤dnear，轉步驟(7)。

(7)根據(jù)公式(24)計算變異函數(shù)，若0

(8)計算人工狼個體對應的目標函數(shù)值Yi，同時對狼群群體進行更新。

(9)判斷算法是否達到尋優(yōu)精度或最大迭代次數(shù)，如果達到則將頭狼位置作為圖像最優(yōu)分割閾值，同時二值化分割圖像，否則返回步驟(3)。

3 實驗仿真與分析

實驗運行環(huán)境為Windows 10系統(tǒng)下的MATLAB R2016a。首先，選取Benchmarks函數(shù)集中的Sphere函數(shù)和Griewank函數(shù)對改進的狼群算法進行尋優(yōu)性能測試。測試函數(shù)如表1所示。

表1 測試函數(shù)

由于算法的隨機性會對計算結果產(chǎn)生較大的影響，因此，對于標準狼群算法(WPA)和改進狼群算法分別獨立運行20次。參數(shù)設置：狼群群體個數(shù)N為30，向量空間維數(shù)D為10；游走次數(shù)上限Tmax=10；探狼比因子α=4，更新因子ζ=6，步長因子s=1 000，距離判定因子w=600。表2是兩種算法在不同測試函數(shù)下的最優(yōu)值和平均值。圖2、3分別是兩種算法在Sphere函數(shù)和Griewank函數(shù)上尋優(yōu)的收斂性能對比。

表2 函數(shù)測試結果

通過對比表2中兩個測試函數(shù)的最優(yōu)值及平均值可知，利用改進的狼群算法求解的全局最優(yōu)值更加接近兩個測試函數(shù)的實際最小值，改進算法的求解精度要比標準狼群算法好，尋優(yōu)能力也更強。

從圖2、3中可發(fā)現(xiàn)，改進狼群算法計算全局最優(yōu)解時所需的迭代次數(shù)比標準狼群算法要小，收斂性能更好，能夠快速求解最佳值。其中適應度值的指數(shù)為求取兩個函數(shù)實際最小值的指數(shù)。因此，可以使用改進狼群算法對三維OTSU法進行優(yōu)化。

圖2 Sphere函數(shù)的收斂曲線

圖3 Griewank函數(shù)的收斂曲線

為驗證提出算法的分割效果，將文中算法與傳統(tǒng)三維OTSU算法以及文獻[10]提出的CWPA算法優(yōu)化三維OTSU法進行圖像分割對比。實驗選擇三幅圖片Bubbles、Baboon和Couple，均為512×512像素。狼群群體規(guī)模N設置為50，空間維數(shù)D為3，最大迭代次數(shù)K=100，CWPA算法的信仰空間接受更新的概率門限為0.7。分割效果如圖4所示，運行結果如表3～表5所示。

圖4 三種算法的分割結果

表3 分割閾值對比

從表中可觀察到，對于Bubbles圖片，文中算法的分割閾值與三維OTSU算法相差不大。而對于其他兩張圖片，文中算法的分割閾值均小于另外兩種分割算法。由于算法將高斯變異引入圍攻行為,不易陷入局部最優(yōu)，能精準找到最佳閾值。

表4 分割用時對比 s

從表中可看出，在時間消耗方面，三維OTSU算法用時最長，文中算法和CWPA算法用時較短。這是因為文中算法采用花授粉法優(yōu)化狼群算法后，加快了全局搜索速度。

表5 區(qū)域一致性對比

區(qū)域一致性反映了算法對圖像的分割效果。在同一類別區(qū)域內(nèi)，像素相似度越大，表明一致性越好，分割效果越佳。區(qū)域一致性可定義為：

(25)

其中，Ri是分割完的第i個區(qū)域，Ii是區(qū)域i的面積，C是圖像的像素點總個數(shù)，f(x,y)是圖像(x,y)處的灰度值。區(qū)域內(nèi)部的均勻性越好則V越大(V的范圍在[0,1]內(nèi))。從表5觀察到，對于Bubbles圖片，三維OTSU算法的區(qū)域一致性較好，但分割時間過長。而對于另外兩張圖片，文中算法的區(qū)域一致性均優(yōu)于CWPA算法和三維OTSU法。因此，文中算法不僅減少了搜索閾值的運行時間，而且提高了分割圖像的精度。

4 結束語

三維OTSU法引入了鄰域中值及均值信息，較好地解決了三維OTSU法對噪聲敏感的問題，但該算法計算量龐大，計算復雜度較高。為此，提出應用改進的狼群算法來搜索最優(yōu)三維閾值向量。實驗結果表明，該算法使計算量大大減少，比CWPA算法的求解能力更強，求解精度顯著提高。