一種改進(jìn)的Otsu多閾值SAR圖像分割方法

(遼寧工程技術(shù)大學(xué) 測繪與地理科學(xué)學(xué)院，遼寧 阜新 123000)

0 引言

主動遙感技術(shù)中的合成孔徑雷達(dá)(synthetic aperture radar，SAR)具有全天時、全天候等光學(xué)傳感器所不具有的圖像采集能力，因此作為其信息載體的SAR圖像已在環(huán)境監(jiān)測，地球資源探測和災(zāi)害評估中得到了廣泛的應(yīng)用[1]。圖像分割將圖像細(xì)分為非重疊的均勻區(qū)域，是SAR圖像處理和應(yīng)用的基本步驟[2]。閾值法由于其簡單直觀的特點(diǎn)，在圖像分割中占有重要地位[3]。在實(shí)際使用中，閾值的選取有多種方法，如直方圖法[4]、最大熵法[5]、Otsu法[6]、模糊聚類法[7]等。其中，Otsu法具有原理簡單、穩(wěn)定有效的優(yōu)點(diǎn)，是使用最廣泛的閾值方法之一。Otsu法由OTSU[8]提出，利用圖像的光譜測度直方圖，使用窮舉法選取使圖像類間方差最大的閾值將圖像劃分為前景和背景2個部分。然而，由于圖像中的對象具有各自不同的光譜分布，因此需要較多的閾值才能從其背景中分割出感興趣的對象。與單閾值Otsu法相比，多閾值Otsu法可以更好地提取圖像目標(biāo)，然而，由于光譜測度直方圖的零階和一階累積矩的大量重復(fù)計算，計算量顯著增加[9]，特別是在需較多閾值的圖像分割中更為突出。

Otsu多閾值圖像分割可視為一個組合優(yōu)化問題，其基本原理為在圖像的光譜測度范圍內(nèi)，根據(jù)特定的準(zhǔn)則和類屬數(shù)搜索系列單閾值的最佳組合。目前，研究人員開始嘗試將優(yōu)化方法和Otsu多閾值分割方法相結(jié)合，以減少計算量，縮短運(yùn)算時間，提高算法的魯棒性。現(xiàn)代優(yōu)化算法主要包括蟻群算法[10]、粒子群算法[11]、遺傳算法(genetic algorithm，GA)[12]和禁忌搜索算法[13]等。其中，GA具有過程簡單、尋優(yōu)能力強(qiáng)、易與其他算法相結(jié)合等優(yōu)點(diǎn)，是一種常用的優(yōu)化算法。GA是模擬自然界的進(jìn)化機(jī)制，對個體由確定的適應(yīng)度函數(shù)逐代演化，以最終演化的個體作為問題的最優(yōu)解[14]。近幾十年來，研究人員已提出了幾種基于GA的圖像分割方法，如Lee等[15]提出的一種基于GA的圖像分割自動閾值選取，以類間方差為進(jìn)化方向，利用GA自動演化能量函數(shù)參數(shù)，使輸出分割結(jié)果滿足定義的感知能量。趙夫群等[16]提出的一種基于GA-Otsu法的圖像閾值分割，其以圖像類間方差為適應(yīng)度函數(shù)；通過適應(yīng)度值確定進(jìn)一步的搜索方向和范圍，最終確定最優(yōu)閾值實(shí)現(xiàn)圖像分割。但由于GA易陷入局部最優(yōu)，出現(xiàn)早熟現(xiàn)象，且迭代過程中計算量大，搜索速度較慢，降低了該類方法的實(shí)際應(yīng)用價值。

量子計算[17]可以解決傳統(tǒng)計算領(lǐng)域的許多難題。Narayanan和Moore[18]首次將量子力學(xué)的概念和原理引入到GA中，提出了量子啟發(fā)遺傳算法，并成功應(yīng)用于求解TSP(traveling salesman problem)問題。由于這種類型的算法不使用量子比特來編碼染色體，因此它仍然屬于傳統(tǒng)意義上的GA。Han和Kim[19]將量子比特和量子門的概念引入到GA，并提出了量子遺傳算法(quantum genetic algorithm，QGA)，其采用量子比特編碼染色體增加種群多樣性和全局尋優(yōu)能力；采用量子門演化量子染色體加快算法收斂，有效地避免早熟現(xiàn)象。目前，QGA已成功用于許多領(lǐng)域，如信號處理[20]、機(jī)組組合問題[21]等。

然而，由于SAR的特殊成像機(jī)制，使SAR圖像具有強(qiáng)相干斑噪聲，降低了圖像質(zhì)量[1]。使得直接對SAR圖像進(jìn)行Otsu多閾值分割的結(jié)果中含有大量的雜散點(diǎn)[22]，這無疑會對后續(xù)的信息提取造成不可逆的錯誤。而提前對圖像進(jìn)行去噪則是一種消除這種現(xiàn)象的有效操作，但傳統(tǒng)的去噪方法都存在不同程度的局限性，如Lee濾波[23]和曲波濾波法[24]在有限的迭代次數(shù)內(nèi)很難去除原圖中無用的紋理信息，無法起到去除大量雜散點(diǎn)的效果，而增加迭代次數(shù)則會導(dǎo)致對象的精細(xì)邊緣全部被平滑掉，降低分割精度。SRAD濾波則可在精確地保持邊緣的前提下以很少的迭代次數(shù)平滑掉紋理信息，最大限度地去除雜散點(diǎn)，特別適用于具有復(fù)雜精細(xì)結(jié)構(gòu)的目標(biāo)。但傳統(tǒng)的SRAD由于無法自動控制迭代次數(shù)使得其會過度平滑掉部分精細(xì)結(jié)構(gòu)。

綜上，為了解決Otsu多閾值法對SAR圖像分割具有的噪聲敏感且計算量大的問題，本文采用圖像間平均結(jié)構(gòu)相似性指標(biāo)來自適應(yīng)地控制SRAD的迭代次數(shù)，以最小迭代次數(shù)達(dá)到既保持精細(xì)邊緣又不過度平滑掉紋理結(jié)構(gòu)的目的。后經(jīng)自適應(yīng)QGA優(yōu)化閾值組合，獲得最佳的圖像分割效果。

1 算法描述

1.1 改進(jìn)的SRAD濾波

設(shè)I0={I0(x，y)，(x，y)∈Ω}為定義在圖像域Ω、大小為M×N的原始SAR圖像，其中，(x，y)表示像素位置，I0(x，y)∈{0，…，L}為像素(x，y)強(qiáng)度的光譜測度，L為圖像的光譜測度量化級。將I0(x，y)作為圖像初始輸入I(x，y；0)，通過偏微分方程建立SRAD模型，可表示為：

(1)

(2)

式中：q(x，y；t)為濾波的瞬時擴(kuò)散系數(shù)，定義為：

(3)

(4)

q0(t)≈q0exp(-ρt)

(5)

式中：ρ為指數(shù)衰減常數(shù)；q0是為初始擴(kuò)散閾值。

對SRAD模型，使用雅可比迭代法進(jìn)行求解。假設(shè)較小的時間步長Δt和在坐標(biāo)x和y方向上較小的空間步長h，則時間和空間坐標(biāo)可被離散為：t=nΔt，n=0，1，2，…；x=ih，i=0，1，2，…，M′-1；y=jh，j=0，1，2，…，N′-1；M′h×N′h=M×N。

(6)

(7)

(8)

(9)

(10)

其中，j∈{0，…，N′-1}，i∈{0，…，M′-1}。

(11)

(12)

其滿足對稱邊界條件，

(13)

其中，j∈{0，1，…，N′-1}，i∈{0，1，…，M′-1}。

圖1 擴(kuò)散示意圖

使用前向差分來近似時間導(dǎo)數(shù)，則SRAD模型的數(shù)值近似為：

(14)

為了在平滑圖像的同時保持圖像的結(jié)構(gòu)特征，在迭代過程中引入平均結(jié)構(gòu)相似性指標(biāo)(mean structural similarity index measurement，MSSIM)[25]作為定量的標(biāo)準(zhǔn)來控制迭代過程。表達(dá)式為：

(15)

(16)

|MSSIM(n)-MSSIM(1)|/MSSIM(1)≤ε

(17)

其中，ε為迭代終止閾值，一般取ε=0.5。

當(dāng)滿足式(17)時，得到平滑圖像z={z(x，y)，(x，y)∈Ω }，其中，z(x，y)為像素(x，y)的光譜測度值，z(x，y)∈{1，…，L}。

1.2 基于自適應(yīng)QGA的Otsu多閾值圖像分割

利用QGA在函數(shù)尋優(yōu)方面具有計算量小，收斂速度快等特點(diǎn)，將其在大小為L的搜索空間中，以圖像類間方差最大為準(zhǔn)則，優(yōu)化得到C個最優(yōu)閾值，并據(jù)此閾值將圖像域Ω分割為C+1個同質(zhì)區(qū)域，即Ω={Ωk，k=1，…，C+1}，其中，k為同質(zhì)區(qū)域索引，Ωk表示第k個同質(zhì)區(qū)域，滿足Ωk∩ Ωk′=?，?k，k′∈{1，…， C+1}，k≠ k′， Ω1∪… ∪ ΩC+1=Ω。在尋優(yōu)過程中，包含以下步驟：量子染色體編碼、量子染色體測量和解碼、適應(yīng)度評價、量子染色體進(jìn)化。

(1)量子染色體編碼。QGA中的染色體用量子比特表示，一個量子比特的狀態(tài)可以表示為：

|ψ〉=α|0〉+β|1〉

(18)

式中：α和β為概率幅，滿足以下歸一化條件：

|α|2+|β|2=1

(19)

式中：|α|2和|β|2分別表示量子比特坍縮到態(tài)|0〉和|1〉的概率。

(20)

(21)

對二進(jìn)制串進(jìn)行解碼，得到閾值的十進(jìn)制值，其原理為在第K個二進(jìn)制串中，第d個二進(jìn)制值1表示十進(jìn)制值2d-1。則閾值與二進(jìn)制串的關(guān)系可表示為：

(22)

(23)

式中：Tk是第k個閾值，其滿足T1

(24)

整個圖像光譜測度級的平均值μ為：

(25)

則基于自適應(yīng)QGA的多閾值Otsu的適應(yīng)度函數(shù)F(·)為：

(26)

圖2 量子旋轉(zhuǎn)門調(diào)整量子比特

(27)

式中：φK，v為第K個染色體中第v個量子比特的旋轉(zhuǎn)角，其轉(zhuǎn)換過程如下：

(28)

因此，第t代種群的更新過程可描述為：

Q(t+1)=G(t)×Q(t)

(29)

式中：G(t)為第t代的量子門；Q(t)和Q(t+1)分別為第t代及其進(jìn)化后第t+1代的種群。

表1 旋轉(zhuǎn)角設(shè)置方案

為了使量子染色體更快地向適應(yīng)度高的個體方向進(jìn)化和避免局部最優(yōu)，當(dāng)量子染色體與歷史最優(yōu)染色體的適應(yīng)度差異較大時，旋轉(zhuǎn)角將適當(dāng)?shù)卦龃螅粗瑢⑦m當(dāng)減小。同時，根據(jù)相鄰兩代量子染色體的差異，逐代地調(diào)整量子旋轉(zhuǎn)角的大小。具體的調(diào)整策略為：

(30)

1.3 算法流程

算法流程圖如圖3所示。

圖3 流程圖

綜上，算法的基本步驟可總結(jié)為：

①對原始SAR圖像進(jìn)行改進(jìn)的SRAD濾波；

②定義Otsu多閾值分割時，圖像的類間方差；

③利用量子比特將閾值組合編碼為種群Q(t)；

④對種群Q(t)進(jìn)行測量和解碼，得到閾值的十進(jìn)制值；

⑤將②的類間方差作為適應(yīng)度函數(shù)，根據(jù)適應(yīng)度大小評價各個個體的優(yōu)劣，保留適應(yīng)度最高的個體作為進(jìn)化目標(biāo)；

⑥利用量子旋轉(zhuǎn)門G(t)作用于Q(t)，得到新種群Q(t+1)；

⑦對新種群Q(t+1)執(zhí)行④、⑤和⑥；

⑧將算法執(zhí)行到最大進(jìn)化代數(shù)或相鄰兩代的最大類間方差的差值小于0.05作為終止條件，如滿足終止條件，則輸出分割結(jié)果，否則繼續(xù)計算。

2 實(shí)驗(yàn)與分析

為了驗(yàn)證該方法的可行性和有效性，對模擬和真實(shí)SAR強(qiáng)度圖像進(jìn)行了實(shí)驗(yàn)。

2.1 模擬SAR圖像

圖4(a)示出了尺度為128像素×128像素的模擬圖像模板，其中，數(shù)字Ⅰ-Ⅴ分別表示不同的均勻區(qū)域。表2示出了對應(yīng)于模擬SAR圖像的均勻區(qū)域的Gamma分布形狀參數(shù)γ和尺度參數(shù)λ，得到的模擬SAR圖像如圖4(b)所示。

圖4 模擬圖像

表2 模擬圖像各同質(zhì)區(qū)域的Gamma分布參數(shù)

利用該方法對圖4(b)進(jìn)行實(shí)驗(yàn)，分割結(jié)果如圖5所示，其中，圖5(a)和圖5(b)分別為SRAD濾波結(jié)果和分割結(jié)果。為了驗(yàn)證提出方法的優(yōu)越性，利用文獻(xiàn)[16]和文獻(xiàn)[22]的方法對圖4(b)進(jìn)行實(shí)驗(yàn)，分割結(jié)果如圖5(c)和圖5(d)所示。從圖5(b)～圖5(d)的比較可以看出，文獻(xiàn)[16]的分割結(jié)果含有較多的分割噪聲，而文獻(xiàn)[22]的方法在區(qū)域Ⅱ-V內(nèi)的分割精度有所提高，但在區(qū)域Ⅰ內(nèi)仍較低。另外，3種方法對圖4(b)的運(yùn)算時間分別為6.35 s、22.52 s、35.81 s，因此，本文方法能較快得到結(jié)果。

圖5 分割結(jié)果

對于定性評估，提取結(jié)果的輪廓線(圖6(a))，將其疊加在圖4(b)的模擬圖像上，如圖6(b)所示。由疊加結(jié)果可以看出，所提輪廓線與實(shí)際輪廓線非常一致，表明本文方法能較好地處理各均勻區(qū)域的邊界。

圖6 視覺評價

對于定量評估，根據(jù)模板圖像分別計算分割結(jié)果圖5(b)～圖5(d)的混淆矩陣[26]，如表3所示。在此基礎(chǔ)上，計算產(chǎn)品精度、用戶精度、總體精度和Kappa值(表4)，其中各項指標(biāo)越高，表明其分割精度越高。

從表4可以看出，本文方法得到的產(chǎn)品精度，用戶精度和總精度均最高，都達(dá)到91.1%以上，Kappa值為0.966，而一般分類器Kappa值達(dá)到0.8以上就被視為高質(zhì)量分類器[26]。上述定量分析結(jié)果表明了該算法的有效性和準(zhǔn)確性。

表3 模擬圖像分割結(jié)果的混淆矩陣

表4 產(chǎn)品精度、用戶精度、總精度和Kappa值

2.2 真實(shí)SAR圖像

選取4幅真實(shí)SAR強(qiáng)度圖像如圖7(a)～圖7(d)所示。圖7(a)和圖7(b)是分辨率為25 m，HV極化、尺度為128像素×128像素的Radarsat-Ⅱ圖像，其中，圖7(a)是河口圖像，人為判讀含有2個同質(zhì)區(qū)域，圖7(b)是城市圖像，人為判讀含有3個同質(zhì)區(qū)域。圖7(c)和圖7(d)是分辨率為30 m、VV極化、尺度分別為128像素×128像素和256像素×256像素的Radarsat-I海冰圖像，人為判讀含有3個同質(zhì)區(qū)域。

從圖7可以看出，原始圖像含有大量斑點(diǎn)噪聲，為此，利用改進(jìn)的SRAD濾波原始SAR圖像，其設(shè)置參數(shù)Δt=0.1；ρ=0.3；h=1，由圖8濾波結(jié)果可發(fā)現(xiàn)，圖像斑點(diǎn)噪聲被很好的濾除且邊緣得到了較好的保留。

圖7 真實(shí)SAR圖像

圖9顯示了4幅真實(shí)SAR強(qiáng)度圖像的分割結(jié)果。其中圖9(a1)～圖9(e1)為文獻(xiàn)[16]的分割結(jié)果；圖9(a2)～圖9(e2)為文獻(xiàn)[22]的分割結(jié)果；圖9(a3)～圖9(e3)為本文方法的分割結(jié)果。從圖9可以看出，文獻(xiàn)[16]的分割結(jié)果中含有較多的分割噪聲，文獻(xiàn)[22]的分割結(jié)果的分割噪聲有所降低，但在邊界處存在誤分割現(xiàn)象，而本文方法在分割噪聲和邊界方面均能得到較好地實(shí)現(xiàn)。次外，對圖7中4幅圖像分割的3種方法的平均運(yùn)行時間分別為23.22 s、35.42 s及6.03 s，該方法的速度明顯高于其他2種方法。

對于定性評估，提取分割結(jié)果的同質(zhì)區(qū)域的邊界線并將其疊加在原始圖像上，如圖10所示。由圖10中可以看出，文獻(xiàn)[16]和文獻(xiàn)[22]中獲得的輪廓線與實(shí)際輪廓線不一致，而本文方法獲得的輪廓線可精確地匹配實(shí)際輪廓。因此，驗(yàn)證了該方法的可行性和有效性。

圖9 分割結(jié)果

圖10 視覺評價

3 結(jié)束語

Otsu法是一個簡單而有效的閾值法，具有良好的性能，被廣泛應(yīng)用于圖像分割。然而，它對噪聲敏感，并且當(dāng)它從單閾值推廣到多閾值時，因?yàn)樾栌嬎忝恳活惖睦鄯e概率和平均值，存在計算量大，耗時過長的缺點(diǎn)，這使其在實(shí)際應(yīng)用中受到限制。為了提高分割精度，同時實(shí)現(xiàn)更好的時間性能，本文提出一種結(jié)合降斑各向異性擴(kuò)散和自適應(yīng)量子遺傳算法的Otsu多閾值SAR圖像分割方法。該方法克服了傳統(tǒng)方法對噪聲敏感以及邊緣細(xì)節(jié)丟失等問題，通過在邊緣兩側(cè)以異性擴(kuò)散來最大限度地去除同質(zhì)區(qū)域的斑點(diǎn)噪聲，并采用平均相似度有效的控制迭代過程，同時兼顧了去除斑點(diǎn)噪聲與保持精度。將改進(jìn)的自適應(yīng)量子遺傳算法作為Otsu多閾值的優(yōu)化工具，快速精確地確定最佳閾值組合，使Otsu多閾值圖像分割更加實(shí)用。但本文需人為確定類別數(shù)，使其分割自動化降低，因此，在未來工作中，擬對本文方法進(jìn)行改進(jìn)，快速地實(shí)現(xiàn)SAR圖像的可變類分割。