利用掩膜主動輪廓模型提取水庫水面面積

苗立新,唐守正,李 霞,紀中奎,文 強,熊志明

（1.中國林業(yè)科學研究院資源信息研究所,北京100091;2.二十一世紀空間技術應用股份有限公司,北京 100096）

密云水庫作為北京惟一的地表水飲用水源,承擔著北京城區(qū)60%的自來水供應任務,對保障北京穩(wěn)定與發(fā)展起著至關重要的作用。近年來,密云水庫入庫水量連年遞減、水面逐年萎縮,對城市的用水安全造成影響,成為人們關注的焦點。水面面積是衡量水庫儲水量的一個重要指標[1]。遙感技術具備有宏觀、實時及成本低等特點,從衛(wèi)星遙感影像中快速、準確地提取水體信息已成為水資源調查、水資源宏觀監(jiān)測及濕地保護的重要手段。北京一號小衛(wèi)星多光譜數據具有時間分辨率高的優(yōu)點,可以滿足對密云水庫進行密集連續(xù)動態(tài)監(jiān)測的需求。但是由于監(jiān)測周期比較短,水庫水面變化不明顯,同時水面提取的精度受系統(tǒng)誤差和人為誤差的影響,水面的實際變化趨勢不能很好地反映出來。為了最大限度地減少系統(tǒng)誤差和人為誤差,提高水庫邊界提取的精度,有必要尋找一種穩(wěn)定的水體提取方法。

針對幾種常見分辨率遙感影像的光譜特點,學者們提出了單波段閾值法[1-2]、譜間關系法[3-4]、色彩空間變換法等多種水體提取方法。例如,劉建波等[1]利用密度分割法從TM影像中提取水體的分布范圍。陸家駒等[2]分別利用閥值法、色度判別法、比率測算法從TM影像中識別水體。楊存建等[3]和鐘春棋等[4]分別以TM影像中TM2+TM3和TM4+TM5的差值及二者的比值輔以適合的閾值來提取水體。王剛等[5]提出兩種水體提取方法:一種是將譜間關系法與IHS彩色空間構建模型相結合的方法,另一種是基于衛(wèi)星圖像數據的LBV變換與歸一化植被指數NDVI的提取方法。都金康等[6]及鄧勁松等[7]利用決策樹方法分別從SPOT 4和SPOT 5影像中提取水體信息;趙書河等人[8]利用迭代混合分析法從中巴資源一號衛(wèi)星數據中提取水體。韓棟等人[9]針對北京一號小衛(wèi)星多光譜數據的特點,引入特征波段PRWI=(B3+B1)/(B3-B1),通過設置不同的PRWI閾值提取水體信息。在以上各種水體提取方法中,一般都存在一個必不可少的步驟,即閾值的選取。由于閾值的選取一般受人的主觀影響比較大,不僅需要反復試驗來確定,而且當閾值選擇不合適時,提取的水體邊界準確度會受到很大影響。同以上方法不同,本文提出的利用掩膜主動輪廓模型提取水體邊界的方法,在水體提取過程中不需要人為選取任何閾值,而是充分利用影像自身內在的信息(梯度信息以及對象間的異質性信息),利用水平集演化的方法提取水庫水面邊界,為水庫水面面積提取方法的探索提供了一種新的思路。

1 掩膜主動輪廓模型

1.1 主動輪廓模型

主動輪廓模型(Active Contour Model)最早由Kass[10]等人提出,用來解決圖像邊緣檢測問題。主動輪廓模型的基本原理是在滿足給定圖像約束的情況下,為了檢測出該圖像中的對象,而讓一條曲線沿著一定的規(guī)則進行變形和移動。這種為了恢復對象的形狀而在數字圖像內進行變形的曲線稱為主動輪廓(Active Contour)。在Kass等人提出主動輪廓模型之后,國內外學者從多個方面對該模型進行了改進,形成多種形式的主動輪廓模型[10-19]。根據主動輪廓模型中停止項是否依賴于圖像的梯度可以將主動輪廓模型分為基于邊緣的主動輪廓模型和基于區(qū)域的主動輪廓模型。Ron Kimme[16]將邊緣與區(qū)域信息結合起來進行目標提取,提出了一種結合了魯棒對齊項(Robust Alignment Term)、加權面積(Weighted Area)、測地主動輪廓(Geodesic Active Contour)和最小方差(Minimal Variance)的主動輪廓模型。

1.2 掩膜主動輪廓模型

本文提出的掩膜主動輪廓模型,即一種帶有掩膜的主動輪廓模型。它主要依據Ron Kimme[16]所提出的主動輪廓模型,并根據實際需要進行了4個方面的改進:

(1)Ron Kimme模型只能從灰度圖像中提取目標,為了使模型能夠處理彩色或多光譜圖像,依據Chen等人的方法對最小方差項進行擴展[17-18],并對所有涉及圖像梯度的項,采用Di Zenzo梯度[20]代替?zhèn)鹘y(tǒng)的灰度梯度,從而將模型的適用范圍從灰度圖像擴展到彩色或多光譜圖像。

(2)根據Chunming Li[21]的不需要重新初始化的水平集思想,增加了一個懲罰項,使水平集函數接近符號距離函數。因此,可以允許采用比較大的時間步長。

(3)Ron Kimme模型及其它各種主動輪廓一般是在迭代開始前確定一個初始的時間步長,這個時間步長在各次迭代過程中保持不變。掩膜主動輪廓模型以一個逐漸收斂的可變時間步長代替固定時間步長,在迭代開始時可以將時間步長設置得比較大,隨著迭代的進行,步長逐漸縮小,從而在加快收斂速度的同時,保證了邊界提取的精度。

(4)在模型中加入了掩膜的概念。由于模型中最小方差項要求待分割的圖像是近似雙峰的,為了滿足這個條件,引入了掩膜的概念,從而得到掩膜主動輪廓模型。掩膜主動輪廓模型需要對待分割圖像和水平集函數同時做掩膜運算,即Im=I?Mask,φm=φ?Mask。初始掩膜取整個圖像范圍,在每次迭代過程中,取水平集內部區(qū)域(φ小于0的部分)為前景對象(水庫水面),對它進行數學形態(tài)學膨脹,從而得到當前迭代中使用的掩膜,代入水平集函數的演化方程進行下一次迭代。隨著迭代的進行,經過掩膜后的待分割圖像越來越接近雙峰圖像,從而滿足最小方差項對圖像的要求。同時,引入掩膜后,由于每次迭代都有一部分明顯不屬于水庫的像元不參加計算,從而計算的速度也大大加快。

改進后的掩膜主動輪廓模型可寫成如下形式:

式中:φm——經過掩膜后的水平集函數;Im——掩膜后的待分割的圖像,可以為灰度、彩色或多光譜圖像;α,β,γ,μ,ν——均為常系數;div ——散度 ;Δ——拉普拉斯操作符;gm(x,y)——邊緣指標函數(或稱邊緣停止函數),這里取

式中:n——迭代次數;tn——第n次迭代的時間步長。

為簡潔起見,記

則掩膜主動輪廓模型的迭代方程可以簡寫為如下形式:

式中:GAC——測地主動輪廓(Geodesic Active Contour)項,它對應的泛函∫g(Cs)ds可以理解為一個加權了的弧長。其中,Cs為參數化曲線,gm(x,y)為邊緣指標函數。當gm為常數1時該泛函表示曲線的弧長。作為對數據敏感的規(guī)則項,GAC用來控制主動輪廓線的光滑程度;AR——魯棒對齊項,它的作用是使圖像梯度方向與曲線法線方向盡可能一致,其中sign為符號函數,對于參數大于0、小于 0、等于0的情況,其返回值分別為1,-1,0;MV——最小方差項,它的作用是使區(qū)域內的方差盡可能小,其中cm1和cm2分別為閉合曲線(零水平集)內部和外部的圖像像元亮度中值(向量);WA——加權面積項,其對應的能量泛函EW(C)=∫∫Ω Cf(x,y)dxdy是曲線圍成的封閉圖形內部區(qū)域的加權面積,其中f(x,y)可以是任何標量函數,當簡單地取 f(x,y)=1時,加權面積項實際上變成了氣球力(Balloon Force)[16,21];P——懲罰項,用來描述一個函數φ接近符號距離函數的程度。根據掩膜主動輪廓模型的迭代公式(1),采用IDL(Interactive Data Language)語言[22]編程實現了掩膜主動輪廓模型,用以支持后續(xù)的試驗。

2 案例分析

2.1 試驗數據

雖然當前可供選擇的遙感數據源已經相當豐富,其中不乏分辨率在米級以及亞米級的數據,如RapidEye數據、Spot5 數據 、QuickBird 數據 、WorldView數據等。但高分辨率數據一般幅寬較窄,重復覆蓋較大區(qū)域的周期比較長,并且多數衛(wèi)星不屬于國內自主運營,數據的可獲取性較差,因此,難以滿足以月度或更短時間周期對密云水庫水面面積進行連續(xù)動態(tài)監(jiān)測的要求。北京一號小衛(wèi)星及運營系統(tǒng),是國家“十五”科技攻關計劃和高技術研究發(fā)展計劃(863計劃)聯(lián)合支持的研究成果,衛(wèi)星上裝有4 m全色和32 m多光譜雙傳感器,其32 m多光譜數據包括綠波段(520～600 nm),紅波段(630～690 nm),近紅外波段(760～900 nm)3個波段,幅寬達600 km,重訪周期3～5 d,能夠滿足對密云水庫進行月度監(jiān)測的需求。4 m全色數據成像幅寬24 km,雖然無法在每個月完全覆蓋密云水庫一次,但可選擇完全覆蓋密云水庫的4 m全色數據作為同期多光譜數據提取結果的驗證數據,驗證多光譜數據上水庫面積的提取精度。研究使用的遙感數據情況見表1。

表1 研究使用的遙感數據

2.2 模型參數設置

在模型迭代過程中,魯棒對齊項實際上對應著圖像中邊緣或梯度所起到的作用;最小方差項對應著對象內部同質性大小;測地主動輪廓項對應著曲線的光滑程度;加權面積項或氣球力起到的主要是加快曲線收縮或膨脹的作用。在這幾個參數中,最重要的是控制邊緣的魯棒對齊項和控制對象內部同質性的最小方差項在所有項中所占的權重,在確定模型參數時,可根據具體圖像的不同特征,給兩者賦予適當的權重系數。

通過取不同組權重系數進行反復試驗,發(fā)現模型中的魯棒對齊項系數、最小方差項系數、測地主動輪廓系數、氣球力系數和初始時間步長的參數值分別設置為5 000,100,100,100,0.01,對試驗中的絕大多數圖像可以取得比較好的提取結果,而且,除非各項系數的值發(fā)生非常大的變化,否則,系數的改變對結果的影響非常小,表現出非常強的穩(wěn)定性。因此,取這組參數作為模型的輸入參數,對包括上述6期多光譜影像在內的一系列多光譜影像進行水庫水面提取。以2007年6月26日數據為例,全色影像上人機交互提取結果和多光譜影像上模型法提取結果對比見圖1。

2.3 試驗結果及分析

基于上述6期北京一號小衛(wèi)星32 m多光譜影像數據,利用掩膜主動輪廓模型提取水庫水面邊界并統(tǒng)計水面面積。同時,采用人機交互解譯方法對同期的北京一號小衛(wèi)星4 m全色影像提取水庫水面邊界并統(tǒng)計水面面積,作為精度驗證的參考數據。將掩膜主動輪廓模型法提取的面積值(簡稱模型法面積)與4 m全色影像提取水庫水面面積(簡稱全色面積)對比,結果見表2。

圖1 2007年6月26日北京一號全色影像人機交互提取結果(a)與多光譜影像模型法提取結果(b)對比

表2 模型法面積與全色面積對比

從表2中可以看出,掩膜主動輪廓模型從6期影像上提取結果的平均精度可達97%以上。說明采用模型法提取的水庫水面面積具有較高的精度。通過計算精度的標準差可知,模型法精度的標準差為1.14%,說明模型法提取的水面面積比較穩(wěn)定,波動較小,比較適合連續(xù)監(jiān)測。另外,模型法面積與全色面積的差值全部為負值,說明模型中仍然存在著一定的系統(tǒng)誤差,使模型的結果比參考結果稍有偏小,經過目視檢查認為,該誤差主要來源于水庫邊界處的混合像元,當需要進一步提高水面面積監(jiān)測精度時,可以考慮采用混合像元分解方法對面積進行進一步修正。

將模型法面積與全色面積(驗證值)生成散點圖,并進行線性擬合,模型法面積與驗證值非常接近,其線性擬合方程為y=1.166x-15.13,R2為0.896 1。

2.4 模型應用

經過驗證,采用掩膜主動輪廓模型從多光譜數據上提取密云水庫水面面積能夠達到比較高的精度,滿足業(yè)務應用的精度要求,同時,各期提取結果的穩(wěn)定性較好。因此,采用該方法對2006年2月到2008年10月質量良好的北京一號多光譜影像進行處理,提取水庫水面邊界并統(tǒng)計面積,實現對密云水庫的月度動態(tài)變化遙感監(jiān)測,從而定量地了解其變化和發(fā)展趨勢,為密云水庫的科學管理、高效利用和有效保護提供理論依據和基礎數據。2006年2月到2008年10月各月份密云水庫水面面積變化曲線如圖2所示。從圖中可以看出密云水庫水面面積從2006年到2008年對應月份總的變化趨勢是逐年減少的。

圖2 2006年 2月至 2008年10月密云水庫各月份水面面積變化圖

3 結論

采用北京一號小衛(wèi)星多光譜數據等重訪周期較短的遙感數據對密云水庫水面面積進行連續(xù)密集監(jiān)測,要求相鄰時相的監(jiān)測結果具有較高的穩(wěn)定性和可比性。尤其是水體邊界變動比較小,要求水體提取方法的誤差相對于不同時相間水體變化必須足夠小,才能真實地反映密云水庫水體面積的變化情況。從試驗結果可以看出,掩膜主動輪廓模型法能比較準確地提取出水庫水面的邊界,滿足連續(xù)監(jiān)測對精度的要求;各期數據提取結果與參考數據相比,精度變化波動小,具有很高的相關性,因此能夠保證連續(xù)監(jiān)測的穩(wěn)定性。

掩膜主動輪廓模型提取水體,只需要在程序運行前設置好參數,提取過程中不需要人工干預,消減了人為誤差的干擾。采用本方法進行連續(xù)動態(tài)監(jiān)測,結果具有較高的穩(wěn)定性和可比性,能夠準確地反映水庫水面面積變化的趨勢,為水庫水面面積自動獲取和連續(xù)監(jiān)測方法研究探索了一條新路。另外,由于遙感圖像上不可避免地存在混合像元,會對提取面積產生一定的影響,如何采用混合像元分解方法進一步提高水面面積監(jiān)測精度是后續(xù)研究需要解決的問題。

[1]劉建波,戴昌達.TM圖像在大型水庫庫情監(jiān)測管理中的應用[J].環(huán)境遙感,1996,11(1):53-58.

[2]陸家駒,李士鴻.TM資料水體識別技術的改進[J].環(huán)境遙感,1992,7(1):17-23.

[3]楊存建,徐美.遙感信息機理的水體提取方法的探討[J].地理研究,1998,17(增刊):86-89.

[4]鐘春棋,曾從盛,柳錚錚.基于譜間特征與比值型指數的水體影像識別分析[J].地球信息科學,2008,10(5):663-669.

[5]王剛,李小曼,田杰.幾種TM影像的水體自動提取方法比較[J].測繪科學,2008,33(3):141-142.

[6]都金康,黃永勝,馮學智,等.衛(wèi)星影像的水體提取方法及分類研究[J].遙感學報,2001,5(3):214-219.

[7]鄧勁松,王珂,鄧艷華,等.SPOT-5衛(wèi)星影像中水體信息自動提取的一種有效方法[J].上海交通大學學報:農業(yè)科學版,2005,23(2):198-201.

[8]趙書河,馮學智,都金康.中巴資源一號衛(wèi)星水體信息提取方法研究[J].南京大學學報:自然科學版,2003,39(1):106-112.

[9]韓棟,楊曉梅,紀凱.小衛(wèi)星遙感影像自動提取水體方法研究[J].測繪科學,2008,33(1):51-54.

[10]Kass M,Witkin A,Terzopoulos D.Snake:Active contour models[J].International Journal of Computer Vision,1988,1(4):321-331.

[11]Cohen L D.On active contour models and balloons[J].CVGIP:Image Understanding,1991,53(2):211-218.

[12]Xu C,Prince J L.Snakes,shapes,and gradient vector flow[J].IEEE Trans.Imag.Proc.,1998,7(3):359-369.

[13]Caselles V,Catte F,Coll T,et al.A Geometric Model for Active Contours in Image Processing[J].Numerische Mathematik,1993,66(1):1-31.

[14]Caselles R V,Kimmel G.Sapiro,Geodesic active contours[J].International Journal of Computer Vision,1997,22(1):61-79.

[15]Chan T F,Vese L A.Active contours without edges,IEEE Trans[J].Image Processing,2001,10(2):266-277.

[16]Ron Kimme.Fast Edge Integration.Geometric Level Set Methods in Imaging[M].New York:Springer,2003:59-77.

[17]Chan T F,Sandberg B Y,Vese L A.Active contours without edges for vector-valued Images[J].Journal of Visual Communication and Image Representations,2000,11(2):130-141.

[18]Rousson M,Deriche R.A Variational Framework for Active and Adaptative Segmentation of Vector Valued Images[M].Geometric Level Set Methods in Imaging,Vision,and Graphics.New York:Springer,2003:195-205.

[19]Cohen L D.On active contour models and balloons[J].CVGIP:Image Understanding,1991,53(2):211-218.

[20]Silvano di Zenzo.A note on the gradient of a multi-image[J].Computer Vision,Graphics and Image Processing,1986,33(1):116-125.

[21]Li Chunming,Xu Chenyang,Gui Changfeng,et al.Level Set Evolution Without Reinitialization:A New Variational Formulation[C]//Proceedings of the 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition.Washington DC,USA:IEEE Computer Society,2005:430-436.

[22]IDL-Data Visualization Solutions.ITT Visual Information Solutions.[2009-12-25].http://www.ittvis.com/idl.