周家香，朱建軍，馬慧云，梅小明

(中南大學 地球科學與信息物理學院，湖南 長沙，410083)

國內(nèi)外在圖像分割方面已取得了很多成果，現(xiàn)有自然圖像分割方法可概括為以下幾個方面：基于全局信息或局部信息；基于區(qū)域或邊緣；針對灰度影像或紋理影像[1]。由于遙感影像所記錄地物的復(fù)雜性和多樣性，數(shù)據(jù)量巨大，將這些分割方法引入到遙感影像分割中面臨很大困難[2?3]。Mean-Shift (MS)算法是一種高效的統(tǒng)計迭代算法，依靠特征空間樣本點的統(tǒng)計特征，不需要任何先驗知識，近年來被廣泛地應(yīng)用于圖像分割[4?7]和跟蹤[8]等計算機視覺領(lǐng)域。傳統(tǒng)MS算法的計算復(fù)雜度是O(kdN2) (其中k為每個數(shù)據(jù)點的平均迭代次數(shù)，d為數(shù)據(jù)維數(shù)，N為1幅圖像的像素)[9?12]，如IKONOS和QuickBird等遙感影像，需要多次迭代才能收斂。同時，若所用的特征向量的維數(shù)多，則收斂速度比較慢，影響遙感圖像的實際應(yīng)用。近年來，一些研究者研究了MS算法的原理，提出減少迭代次數(shù)和降低每次迭代的計算量這2種提高算法效率的方法[10]。為了提高MS算法對遙感影像分割的速度和效果，本文提出以下改進措施：

(1) 對不同空間分辨率遙感圖像給定空間帶寬參考值，不同波段采用不同的灰度值域帶寬?；趥鹘y(tǒng)MS算法編寫的EDISON軟件需要人工輸入空間帶寬和值域帶寬，且所有波段采用相同的值域帶寬。若將遙感影像分割成有意義的區(qū)域且沒有過多的細節(jié)，自適應(yīng)帶寬才能得到好的分割效果。

(2) 面對大數(shù)據(jù)量的遙感影像分割速度慢的問題，本文主要在MS運行效率上采取多項加速措施。

(3) 針對 MS分割存在過分割現(xiàn)象，將檢測出的模點采用基于全局模點融合來得到穩(wěn)定的分割結(jié)果。

(4) 可用于遙感影像分割的特征很多，模點檢測或模點融合時利用全部特征勢必增加計算時間，改進方法將可將特征分組，利用位置和顏色特征進行模點檢測，利用紋理特征進行模點融合。

1 改進MS遙感影像分割算法

1.1 傳統(tǒng)高斯MS算法模點檢測

傳統(tǒng)高斯核MS由高斯核密度估計的梯度推導(dǎo)獲得[4]。給定d維歐拉空間R中n個采樣點{xi，i=1，…，n}，利用高斯核函數(shù)K (x)= ( 2π)?d/2exp(?x2)，則MS迭代向量為：

式中：xs為二維位置特征；xr為多維值域特征；hs為空間帶寬；hr為值域帶寬。則傳統(tǒng)MS計算步驟：

(1) 選用高斯核函數(shù)，帶寬hs和 hr，隨機給定樣本點x；

(2) 根據(jù)式( 1) 計算樣本點的迭代向量m(x)；

開展拓展實驗是培養(yǎng)學生科學探究能力的有效補充途徑。目前高中生物教材中的很多實驗為體驗或驗證性實驗，即使是探究性實驗，也列出了詳細的探究過程與具體內(nèi)容，學生往往在未經(jīng)充分思維后就直接參考了教材。因此，教材實驗在培養(yǎng)學生自主探究能力上還稍顯不足，不能滿足熱愛生物學科的學生的求知需求和探究愿望，而拓展實驗可以作為有效補充。為此，在教學實踐中，筆者結(jié)合校本課進行了拓展實驗教學，積極探索有效開展拓展實驗促進學生科學探究能力提升的實施策略。

(3) 條件判斷，若|m(x)?x|＞ε(ε為閾值，如 0.001)，則移動點x到m(x)，執(zhí)行步驟(2)，否則，結(jié)束循環(huán)；

(4) 將結(jié)束循環(huán)后的 m(x)賦給 x，此時 m(x)稱為樣本點x的模點；

(5) 移動到下一樣本點x，重復(fù)步驟(2)～(4)，直至檢測出影像內(nèi)每個數(shù)據(jù)點的模點為止；

(6) 將收斂至相同模點的數(shù)據(jù)點歸為同一區(qū)域，得到初分割結(jié)果。

這樣，整個影像被分解成m個子區(qū)域Sj(j=1，2，…，m)，每個子區(qū)域像素為nj，中心為cj，則：

1.2 全局模點融合

由于MS分割存在過分割問題，經(jīng)過MS分割后，將以點陣表示的圖像轉(zhuǎn)化為塊狀表示。然后，采用如下過程的基于區(qū)域的MS迭代。

對每一個子區(qū)域Sj(j=1，2，…，m)，進行以下迭代：

式中：np為區(qū)域p的像素；h為帶寬；k為迭代次數(shù)；= c 。用上式進行迭代計算，直至j 為止(εw是閾值，如0.001)，將最終收斂值記錄為。

1.3 改進的MS遙感影像分割算法實現(xiàn)

1.3.1 濾波

圖像噪聲會使檢測出的局部模點數(shù)過多，這樣使后面的模點融合變得困難。濾波能降低噪聲對模點檢測的影響。事實上，MS本身對噪聲不太敏感[11]，模點檢測前的濾波能很好地降低噪聲的影響，這樣能使模點檢測過程更簡單。本文選擇3×3的中值濾波。

1.3.2 帶寬的選擇

在MS算法中最重要的參數(shù)是帶寬矩陣H，它決定迭代速度和分割效果。目前，Plug-in 和Cross-validation是廣泛使用帶寬選擇方法[1]。Plug-in規(guī)則是用估計參數(shù)來代替未知參數(shù)，漸進積分均方誤差是應(yīng)用于Plug-in規(guī)則的通用標準。使?jié)u進積分均方誤差最小便可估計全局值域帶寬：

式中：d為特征向量維數(shù)；j為單個波段；n為數(shù)據(jù)量；σj為特征向量的標準差。

在plug-in準則下獲得的值域帶寬hr=(hR，hG，hB)能盡可能地降低估計密度函數(shù)與原始密度函數(shù)之間的變形。為了將遙感影像分割成有意義的區(qū)域而不含有過多的細節(jié)，一般來說，過小的hr是無意義的。但是，過大的值域帶寬導(dǎo)致忽略更多的影像細節(jié)，同時降低檢測出的模點數(shù)；相反，小的灰度帶寬使更多的地物細節(jié)被提取出來。因此，灰度值域帶寬hr被認為是分割的尺度。

1.3.3 快速高斯模點檢測

應(yīng)用傳統(tǒng)MS算法進行聚類，可得到初始分割結(jié)果。初始分割是對影像中每一個像點進行MS迭代得到其模點，由于遙感影像數(shù)據(jù)量大，迭代計算量非常大，本文采用以下加速策略。

(1) 以待處理點 xs為中心，事先根據(jù)圖像的空間分辨率來選擇空間帶寬hs，以長、寬均為2hs+1的正方形按下式構(gòu)造高斯空間模板；

式中：xs表示待處理點的二維坐標，xsi為以xs為中心，長和寬為2hs+1正方形內(nèi)的采樣點二維坐標?？梢钥醋鳈?quán)值，與xs點近的采樣點權(quán)值大，遠的采樣點權(quán)值小。則式(1)變成：

(2) 為了減少MS迭代次數(shù)，從xn點開始MS迭代，第 τ次迭代后的向量為等于或近似于另一特征向量xm，則認為xn和xm收斂于同一個模點，停止迭代。

(3) 對任意一隨機樣本點x，沿一定的路徑收斂至一固定點(即模點)，將收斂路徑一定寬度內(nèi)的點都認為收斂至相同的模點，這些點不再參加迭代計算。

(4) MS迭代計算量大的1個主要原因是特征向量的維數(shù)，而遙感圖像可用的特征向量的維數(shù)很多，本文將特征向量分為2組，首次迭代時使用多光譜灰度特征，模點融合時使用紋理特征。

1.3.4 基于全局的模點融合

由于傳統(tǒng)MS分割一般存在過分割現(xiàn)象，進一步將分割后的影像進行全局模點融合，得到最終分割結(jié)果。遙感影像的特征有位置、灰度、紋理、形狀等，為了達到理想的分割結(jié)果，應(yīng)盡可能地利用遙感影像的多維特征[9]，這樣勢必會增加計算時間。為了充分利用遙感影像的多維特征又不增加計算時間，在模點融合時使用紋理特征。

此時，用到的紋理特征f1，f2和f3分別為：

式中：i和j為像素點灰度值；p(i,j)為灰度值i和j出現(xiàn)的歸一化概率；L為灰度級。則式(3)中h應(yīng)為紋理帶寬ht。對任一分割子區(qū)域Sj(j=1，2，…，m)，其紋理帶寬由下式計算：

式中：nj為第 j塊的像素，fki(k=1，2，3)為該塊內(nèi)像素i的 fk紋理特征；ckj(k=1，2，3)為第j塊的fk紋理均值。

2 實驗結(jié)果評價與分析

2.1 分割實驗

在改進的分割方法中，空間帶寬hs主要影響計算速度。不同傳感器的影像有不同的空間分辨率。圖1所示為 SPOT5影像采用不同空間帶寬分割結(jié)果的對比。實驗發(fā)現(xiàn)，空間分辨率為10.00 m的SPOT5影像選用hs=5比較合適，空間分辯率為2.44 m的Quickbird影像選用hs=9比較合適。

圖2(a)是1幅Quickbird 影像270×270的子圖像，采用前面提出的參考空間帶寬 hs=9和由式(4)計算的灰度值域帶寬進行快速高斯MS迭代得到圖2(b)所示的分割結(jié)果(290個模點)。圖2(c)所示是用EDSION軟件(1個基于傳統(tǒng)MS算法編寫的軟件)分割結(jié)果。將本文的分割結(jié)果與EDISON軟件分割結(jié)果進行對比，由圖2(c)可見：本文分割結(jié)果要優(yōu)于EDISON軟件分割結(jié)果。

圖1 SPOT5影像不同空間帶寬對比實驗Fig.1 Contrast experiment of SPOT5 image with various space bandwidths

圖2 改進分割方法與EDISON分割結(jié)果對比Fig.2 Comparative segmentation between presented approach and EDISON of Quickbird image

從圖2(b)可以看出：本文的分割結(jié)果比較理想，但與分割參照圖(圖3(b))比較還存在過分割現(xiàn)象。對初分割結(jié)果進行全局模點融合，保持原有空間帶寬不變，采用式(11)計算每塊的紋理帶寬，得到圖3(a)所示的實驗結(jié)果。由圖3(c)可見：本文分割結(jié)果比較理想。

2.2 本文改進方法與傳統(tǒng)MS算法的時間比較

通過實驗比較傳統(tǒng)MS方法與本文MS改進方法的算法時間。采用圖2(a)所示Quickbird圖像，計算時間見表1，從表1可以看出：本文所采用加速策略是有效的。

2.3 分割效果討論

為了評價分割的效果，本文對圖2(a)所示的Quickbird影像進行定量對比實驗。Polak等[13]提出了一種基于對象一致性誤差(EOC)的客觀有效的圖像分割評判標準。EOC(Ig，Is)∈[0，1](其中Ig為地面參考對象，Is為實際分割對象)。當EOC=0時，表示分割完全與地面參考對象一致；EOC越小，表示分割效果越好。實驗所用圖2(a)的研究范圍較小，圖像空間分辨率較高，且各類地物相對集中，通過目視判讀得到實驗參考圖(圖3(b))。將地物分為房屋、林地、水面、農(nóng)田、綠地以及陰影6類。表2所示為本文方法與EDISON軟件分割的EOC，可以看出：采用本文方法所得房屋、林地的分割效果明顯優(yōu)于EDISON的分割效果，水面、農(nóng)田、綠地和陰影4類地物的分割效果相當。

圖3 模點融合后分割結(jié)果Fig.3 Segmentation results after mode merging

表1 傳統(tǒng)MS與改進MS計算時間比較Table 1 Comparative time of traditional MS and improved MS s

表2 改進分割方法與EDISON分割方法的Eoc對比Table 2 Comparison of Eoc between presented improved approach and EDISON

3 結(jié)論

(1) 對常用的高分辨率遙感影像如 SPOT5和Quickbird的MS模點檢測的空間帶寬的選擇給出了參考值，其他的遙感影像根據(jù)其空間分辨率參照這2類影像來確定。

(2) 遙感影像不同波段分別用公式計算值域帶寬，以實現(xiàn)自適應(yīng)的MS遙感影像分割。

(3) 模點檢測時采用 4個有效的加速措施，減少了MS迭代時間；同時，將遙感影像的多維特征分組，模點檢測時采用空間特征和多光譜特征，模點全局融合時采用紋理特征，這樣加快了MS迭代速度。

(4) 采用基于全局的模點融合提高遙感影像分割精度。用 MATLAB實現(xiàn)了本文改進的遙感圖像分割方法，實驗結(jié)果與由傳統(tǒng)算法所得結(jié)果相比，本文方法具有分割效果好等優(yōu)點。

[1]HONG Yi-ping, YU Jian-qiang, ZHAO Dong-bin. Improved mean shift segmentation approach for natural images[J]. Applied Mathematics and Computation, 2007, 185(2): 940?952.

[2]van der Sande C J, de Jong S M, de Roo A P J. A segmentation and classification approach of IKONOS-2 imagery for land cover mapping to assist flood risk and flood damage assessment[J]. International Journal of Applied Earth Observation and Geoinformation, 2003, 4: 217?229.

[3]楊朝云, 陳光儒, 呂嫦艷, 等. 一種改進的高分辨率遙感影像分割方法及應(yīng)用[J]. 電腦信息技術(shù), 2010, 18(4): 45?48.YANG Zhao-yun, CHEN Guang-ru, Lü Chang-yan, et al. An improved high resolution RS images partition way and its applying[J]. Computer and Information Technology, 2010, 18(4):45?48.

[4]Comaniciu D, Meer P. Mean-Shift: A robust approach towards feature space analysis[J]. IEEE Transaction on Pattern Analysis and Machine Intelligence, 2002, 24(5): 603?619.

[5]Comaniciu D, Ramesh V, Meer P. The variable bandwidth mean shift and data-driven scale selection[C]//Proceedings of Eighth IEEE International Conference on Computer Vision. Vancouver,United States, 2001: 438?445.

[6]王兆虎, 劉芳, 焦李成. 一種基于視覺特性的遙感圖像分割[J]. 計算機學報, 2005, 28(10): 1687?1691.WANG Zhao-hu，LIU Fang，JIAO Li-cheng. A remote sensing image segmentation algorithm based on vision information mean shift algorithm accelerated mean shift algorithm[J]. Chinese Journal of Computers, 2005, 28(10): 1687?1691.

[7]Comaniciu D. An algorithm for data-driven bandwidth selection[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2003, 25(2): 281?288.

[8]Comaniciu D, Ramesh V, Meer P. Real-time tracking of non-rigid objects using mean shift[C]//Proceedings of IEEE Conference on Computer Vision and Pattern Recognition. Hilton Head Island, United States, 2000: 142?149.

[9]Georgescu B, Shimshoni I, Meer P. Mean shift based clustering in high dimensions: a texture classification example[C]//Proceedings of IEEE International Conference on Computer Vision. Nice, France, 2003: 456?463.

[10]Carreira-Perpinan M A. Acceleration strategies for Gaussian mean—shift image segmentation[C]//Proceedings of the 2006 IEEE Computer Society Conference on Computer Vision and Pattern Recognition. New York: IEECS Press, 2006: 543?549.

[11]Kaftan J N. Mean shift segmentation evaluation of optimization techniques[C]//International Conference on Computer Vision Theory and Application. Madeira, Portugal, 2008: 365?374.

[12]Carreira-Perpinan M A. Gaussian mean shift is an EM algorithm[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2007, 29(5): 767?776.

[13]Polak M, ZHANG Hong, PI Ming-hong. An evaluation metric for image segmentation of multiple objects[J]. Image and Vision Computing, 2009, 27: 1223?1227.