基于區(qū)域統(tǒng)計的遙感地物半自動提取

張 煜，陳文龍，魏思奇，葉 松，曹 波

(長江科學(xué)院 空間信息技術(shù)應(yīng)用研究所，武漢 430010)

?

基于區(qū)域統(tǒng)計的遙感地物半自動提取

張 煜，陳文龍，魏思奇，葉 松，曹 波

(長江科學(xué)院 空間信息技術(shù)應(yīng)用研究所，武漢 430010)

由于遙感影像的復(fù)雜性，面狀地物的提取存在著一定的難度。區(qū)域統(tǒng)計活動輪廓法有著較強(qiáng)的穩(wěn)定性和適應(yīng)性，且在一定的概率密度分布下存在著快速的算法。針對交互式操作的特點，提出了一種基于多邊形區(qū)域統(tǒng)計活動輪廓的面狀地物半自動提取方法，并使用最小描述長度法簡化面狀地物的輪廓。實驗結(jié)果顯示了該方法在面狀地物提取問題上的有效性和實用性。

統(tǒng)計活動輪廓；面狀地物；區(qū)域提取；隨機(jī)搜索；最小描述長度

1 研究背景

圖像分割與目標(biāo)輪廓提取是遙感信息提取的重要任務(wù)之一，常常作為其較低層次的前處理部分，有著廣泛的應(yīng)用價值。湯育紅[1]在地理國情普查中，基于遙感影像的地表覆蓋分類提取，需要針對不同地物類型對地表影像進(jìn)行分割，獲取不同地物和人文要素的分布；黃俊等[2-3]基于遙感的水土保持監(jiān)測需要通過遙感分類對不同地物覆蓋的逐年變化進(jìn)行分析，實現(xiàn)測區(qū)內(nèi)水土侵蝕類型、等級、流失面積、覆蓋變化等信息的提取。然而，盡管存在著眾多的分割提取方法，卻沒有一個通用的方法來解決所有的圖像分割任務(wù)，只能針對不同的實際任務(wù)設(shè)計不同的算法。基于交互式的半自動提取通常被認(rèn)為是一種充分結(jié)合人工目視識別方式實現(xiàn)快速可靠的分類提取方法，是實際遙感處理工作中常采用的有效的操作方法。

本文基于區(qū)域統(tǒng)計的活動輪廓法[4]，增加了可變面積項，并針對半自動交互式操作的特點，提出了基于多邊形統(tǒng)計區(qū)域的實現(xiàn)方法，同時使用最小描述長度法描述和簡化影像目標(biāo)的輪廓。實驗結(jié)果顯示了本文方法對遙感影像中面狀地物提取和分割有著良好的效果。

2 區(qū)域統(tǒng)計活動輪廓及其快速算法

圖1 輪廓區(qū)域的定義Fig.1Definition of contour

對于給定的影像I={I(x，y)|(x，y)∈[1，Nx]×[1，Ny]}，如圖1所示，假設(shè)只有一個感興趣的目標(biāo)(圖中陰影部分)，實線多邊形為輪廓的初始位置。輪廓將影像分為輪廓內(nèi)區(qū)域a和輪廓外區(qū)域b，在影像灰度為高斯分布條件下，區(qū)域統(tǒng)計活動輪廓的能量函數(shù)[4]為

(1)

式中：Na和Nb分別為a，b區(qū)域的像素數(shù)量；σa和σb為a，b區(qū)域內(nèi)像素的標(biāo)準(zhǔn)差。而T被定義如下(以區(qū)域a為例)：

(2)

(3)

(4)

式中：針對前述的3個累加量，k分別為0，1，2；圖2(a)顯示了不同參數(shù)的解釋，ymax,ymin分別為輪廓在y方向上的最大、最小坐標(biāo)；而xmax，xmin分別為當(dāng)前y值下x方向上的最大值和最小值；δ為輪廓；c(x,y)為系數(shù)，其值取0，-1，1。

圖2 輪廓示意、節(jié)點Freeman碼和c(x,y)值Fig.2 Illustrations of contour, the direction of Freeman code and the value of c(x,y)

式(4)是實現(xiàn)快速算法的關(guān)鍵，可以在算法迭代前預(yù)先處理而求出所有的F[x,y]，而只在迭代過程中通過不同的x,y查找出F[x,y]值即可。于是，對于二維的雙重累加被輪廓線上一維累加代替。而系數(shù)c(x,y)取1還是-1實際上分別對應(yīng)著在y取固定值下輪廓線上x方向上最大值和最小值。輪廓線的編碼采用如圖2(b)所示的Freeman碼，通過對于某一點的鏈入和鏈出(即上一節(jié)點到當(dāng)前節(jié)點和當(dāng)前節(jié)點到下一節(jié)點的鏈碼值)，可以建立一個二維表格，針對不同的鏈入鏈出值取不同的c(x,y)值，于是跟蹤輪廓可以獲得所有的c(x,y)值，如圖2(a)中1所指示的矩形區(qū)域，其鏈表結(jié)構(gòu)及c(x,y)值如圖2(c)所示。于是式(3)可以通過跟蹤輪廓快速地求得，并最終能夠得到區(qū)域a中的3個累加量。因此，實際的二維掃描只在算法迭代前進(jìn)行了1次，以后的每次迭代計算累加量都是沿輪廓進(jìn)行的，這樣大大提高了計算效率。該算法的具體流程和詳細(xì)說明可以參見文獻(xiàn)[4]。

3 基于多邊形區(qū)域統(tǒng)計活動輪廓的目標(biāo)提取

前述方法中區(qū)域b為區(qū)域a以外的整個影像，然而在實際的影像分割任務(wù)，目標(biāo)地物與背景是很難自動區(qū)分的，在多數(shù)情況下，目標(biāo)地物附近可能存在多個其它地物，它們之間的相互位置也不是固定的，因此需要人工目視給出種子點。

本文針對此特點設(shè)計了對于任意單連通多邊形區(qū)域的統(tǒng)計活動輪廓目標(biāo)提取算法，人工給定的簡單初始輪廓將包含完整的感興趣面狀地物以及少量的背景，這樣做的好處為：①人為避開干擾物；②減少數(shù)據(jù)處理量，提高計算速度。其整體的算法流程如圖3所示。

圖3 算法流程Fig.3 Flow chart of the algorithm

通過種子點給出背景多邊形區(qū)域，該區(qū)域包含完整的希望提取出的某個地物，因此必須將該多邊形標(biāo)記出來，如圖4所示，矩形區(qū)域為多邊形的包圍盒，式(4)的F[x,y]在給出種子點后無需在整個影像區(qū)域內(nèi)計算，而只用在被標(biāo)記的多邊形區(qū)域內(nèi)進(jìn)行。

圖4 多邊形區(qū)域示意圖Fig.4 Illustration of the polygon

多邊形的直線段繪制采用數(shù)值微分法，在繪制的同時建立Freeman碼，而后基于二維查表法建立c(x,y)值鏈表，多邊形區(qū)域的標(biāo)記直接借用了c(x,y)值鏈表結(jié)構(gòu)。如圖4所示，建立與圖中包圍盒大小相同的標(biāo)記區(qū)，掃描包圍盒標(biāo)記區(qū)，開始時標(biāo)記區(qū)內(nèi)每點值為0，若掃描線經(jīng)過c(x,y)鏈表值-1，則開始將標(biāo)記區(qū)內(nèi)設(shè)置為1，若再掃到c(x,y)鏈表值為1則停止設(shè)置。如此反復(fù)，最終建立起多邊形區(qū)域的標(biāo)記。

此外，種子點給出的多邊形區(qū)域不僅要包含目標(biāo)地物，也要容納一定數(shù)量的背景像素，這就導(dǎo)致了種子點遠(yuǎn)離了感興趣的目標(biāo)地物，盡管活動輪廓演變時利用隨機(jī)搜索法，但在遠(yuǎn)離前景的地方仍然容易陷入局部極小值。本文使用了可變面積項使活動輪廓在演變的初始階段具備較強(qiáng)的向內(nèi)壓力而收縮[5]，而在演化后期該壓力逐漸減少，即

(5)

式中：r與隨機(jī)搜索范圍相關(guān)，而該范圍隨后面介紹的多級循環(huán)次數(shù)的增加而逐漸減??；Na+b為多邊形區(qū)域內(nèi)的總像素數(shù)量，迭代過程中保持不變；Na為當(dāng)前輪廓內(nèi)像素數(shù)量。

活動輪廓演化的過程使用了隨機(jī)最小化方法，具體步驟如下：

(2) 進(jìn)行一次某個輪廓點的移動，等概率隨機(jī)選擇某個輪廓點，將該點由開始位置(x,y)移動到(x+dx，y+dy)，其中dx，dy的變化范圍為(-r，r)，而dx或dy的取值也是等概率的在范圍(-r，r)中選取。

(3) 進(jìn)行移動點的交叉測試和有效區(qū)域測試，交叉測試確保輪廓在演化過程中不會出現(xiàn)自相交，而有效區(qū)測試通過使用多邊形標(biāo)記區(qū)確保輪廓節(jié)點的移動不會超出多邊形區(qū)域的邊界。如果測試不通過則放棄當(dāng)前的移動。

(4) 依式(5)，如果當(dāng)前計算出的E′小于以前的則移動有效，否則放棄移動。

將以上步驟(2)—步驟(4)迭代指定次數(shù)后為一次循環(huán)。由于人工給定的初始點通常較少，直接使用這些控制點作為輪廓點無法滿足復(fù)雜地物的邊緣擬合，算法采用多尺度方法，具體的實現(xiàn)過程是：第1次循環(huán)的多邊形為人工給定的少量種子點，此時由于種子點離前景較遠(yuǎn)，可以使用較大的r值范圍，即較大的隨機(jī)搜索步伐，于是一個近似的較少點的邊緣輪廓產(chǎn)生且盡可能地與目標(biāo)地物吻合;第2次循環(huán)，活動輪廓的節(jié)點將根據(jù)特定的閾值th而增加，凡是多邊形的某線段的距離大于th的就在該線段中間增加一個節(jié)點，如此反復(fù)地增加直到所有的邊長度都小于該閾值th，然后再運行隨機(jī)搜索算法，而活動輪廓隨機(jī)搜索循環(huán)中使用的搜索范圍r也依循環(huán)次數(shù)相應(yīng)地減小。后續(xù)的循環(huán)過程與第2次循環(huán)類似，th與r都相應(yīng)減小。一般而言，經(jīng)過3—4次循環(huán)，活動輪廓可以收斂到較好的結(jié)果。

圖5顯示的是4次循環(huán)迭代結(jié)果，其中，圖5(a)為人工給定初始多邊形輪廓，共6個種子點，經(jīng)過300次隨機(jī)動態(tài)演化的循環(huán)，6點形成的多邊形緊靠目標(biāo)地物的邊緣。注意圖5(b)中輪廓頂端出現(xiàn)的尖點，由于考慮活動輪廓隨機(jī)演變效率而未采用平滑項，尖點出現(xiàn)是難免的，因此每次迭代完成后都需要檢測和刪除這些尖點，也可以使用最小描述長度法在簡化輪廓的時候自動去除。

圖5 多尺度方法不同循環(huán)階段結(jié)果Fig.5 Results of different stages of multiple scale method

4 基于最小描述長度準(zhǔn)則的輪廓簡化

為了達(dá)到與目標(biāo)地物邊緣一致，活動輪廓的節(jié)點數(shù)量越多越好，但是節(jié)點數(shù)量過多對于輪廓的存儲和描述都不利。本文采用了最小描述長度法來簡化輪廓。最小描述長度法是由J.Rissanen[6]提出，用來描述數(shù)據(jù)的平均復(fù)雜度，該準(zhǔn)則認(rèn)為，評估一組數(shù)據(jù)復(fù)雜度的最好方法是同時考慮對數(shù)據(jù)模型的編碼碼長以及模型參數(shù)本身的編碼碼長。本文中，數(shù)據(jù)模型編碼包括對區(qū)域a和b的編碼，而參數(shù)編碼即是對活動輪廓節(jié)點位置的編碼，其總的碼長可以表達(dá)為

(6)

式中：Pa與Pb分別為a，b區(qū)域像素灰度值分布；k為活動輪廓節(jié)點數(shù)量；L()為總碼長，在影像灰度為高斯分布條件下，由式(1)的推導(dǎo)過程，基于最小描述長度準(zhǔn)則的能量函數(shù)可以寫為

Nblgσb+klgNa+b=E+klgNa+b。

(7)

圖6 簡化前后輪廓疊加顯示Fig.6 MDL contour superimposed on contour without MDL

可見，使用基于最小描述長度準(zhǔn)則只需在原來的能量函數(shù)基礎(chǔ)上增加一項活動輪廓的編碼長度即可。實際應(yīng)用時，該準(zhǔn)則并不適合于一開始就使用，這樣會導(dǎo)致在節(jié)點數(shù)量不多時無法準(zhǔn)確匹配邊緣，有效可行的方案是先使用式(5)過度增加節(jié)點，然后基于式(7)逐個刪除多余節(jié)點，即每次刪除能使EMDL下降最快的輪廓點。圖6顯示了采用EMDL對輪廓進(jìn)行簡化后與使用E’的簡化前的輪廓疊加，可見二者的差別非常小，但是簡化前輪廓節(jié)點數(shù)為173，簡化后為51。

5 實驗與驗證

本文實驗數(shù)據(jù)為航空影像，在假設(shè)影像灰度為高斯分布的基礎(chǔ)上，使用式(5)和式(7)并基于文中的方法對圖像中的面狀地物進(jìn)行提取。圖7顯示的是人工給定的初始多邊形(即提取前的圖像)和對應(yīng)的提取結(jié)果。對于相對孤立的面狀地物如圖7(a)和圖7(b)，提取可直接使用矩形作為初始多邊形，而圖7(c)和圖7(d)的圖像前后景關(guān)系復(fù)雜，更適合采用多邊形以獲得正確的結(jié)果。圖中的面狀地物普遍沒有明顯的邊緣，前后景的差別也不是特別清晰，直接使用基于邊緣或基于簡單的均勻區(qū)域分割的方法都很難將感興趣的區(qū)域提取出來。

圖7 影像提取前后的實驗結(jié)果Fig.7 Experiment results before and after image extraction

6 結(jié) 論

本文所用的多邊形區(qū)域統(tǒng)計活動輪廓對這類物體的提取有著很好的適應(yīng)性，算法的抗干擾能力較強(qiáng)，對種子點的位置也沒有嚴(yán)格的要求，算法的收斂經(jīng)過了4次循環(huán)，分別經(jīng)過了300，400，500，600次迭代，節(jié)點隨機(jī)移動范圍以像素為單位分別是20，10，4，2，分級閾值th除了第1次由于人工給定而無th外，其它3次分別是30，10，5。實驗結(jié)果顯示了本文方法的有效性和實用性。

[1] 湯育紅. 地理國情普查地表覆蓋與國情要素信息的提取方法探討[J]. 測繪與空間地理信息，2013，36(12):89-91.

[2] 黃 俊，申邵洪. 結(jié)合高分辨遙感影像和GIS數(shù)據(jù)的土地利用變化監(jiān)測[J]. 長江科學(xué)院院報，2012，29(1):49-52.

[3] 孟廣濤, 方向京, 和麗萍，等. 3S技術(shù)在水土保持動態(tài)監(jiān)測中的應(yīng)用[J]. 水土保持研究，2007，14(2)：8-10.

[4] CHESNAUD C, REFREGIER P, BOULET V. Statistical Region Snake-based Segmentation Adapted to Different Physical Noise Models[J]. IEEE Transactions on Pattern Analysis and Machine Intelligence，1999，21(11)：1145-1157.

[5] RONFARD R. Region-based Strategies for Active Contour Models[J]. International Journal of Computer Vision, 1994, 13(2):229-251.

[6]RISSANEN J.Hypothesis Selection and Testing by the MDL Principle[J]. Computer Journal, 1999, 42(4)：260-269.

(編輯：黃 玲)

Target Extraction Based on Polygon StatisticalRegion-based Active Contour

ZHANG Yu，CHEN Wen-long，WEI Si-qi，YE Song，CAO Bo

(Spatial Information Technology Application Department, Yangtze River Scientific Research Institute, Wuhan 430010, China)

Active contour is an important method to segment images. Its segmentation strategy is usually edge-based or region-based and based on the contour energy minimization. But the classical active contour is difficult to deal with the region which has weak edges or is in the complicated background. In this paper, we implement an improved statistical region-based active contour and propose a strategy to extract region of interest by polygon. Within the family of exponential probability density function, statistical region-based active contour can run by fast algorithm. We give details of the realization of the algorithm and simplify the contour by the minimum description length. Experiment results illustrate the efficiency of our proposed method in contour extraction.

statistical active contour; planar features; region extraction; random search; minimum description length

2016-08-20

云南省水利重大科技項目(CKSK2015852/KJ)；云南省省級水資源費項目(CKSK2015720/KJ)

張 煜(1971-)，男，山東陽谷人，高級工程師，博士，研究方向為攝影測量與遙感，(電話)18986061273(電子信箱)zhangyu_1999@126.com。

10.11988/ckyyb.20160855

2016,33(11):12-16

P237

A

1001-5485(2016)11-0012-05