苗立新,唐守正,李 霞,紀(jì)中奎,文 強(qiáng),熊志明

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

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

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

1 掩膜主動(dòng)輪廓模型

1.1 主動(dòng)輪廓模型

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

1.2 掩膜主動(dòng)輪廓模型

本文提出的掩膜主動(dòng)輪廓模型,即一種帶有掩膜的主動(dòng)輪廓模型。它主要依據(jù)Ron Kimme[16]所提出的主動(dòng)輪廓模型,并根據(jù)實(shí)際需要進(jìn)行了4個(gè)方面的改進(jìn):

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

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

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

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

改進(jìn)后的掩膜主動(dòng)輪廓模型可寫成如下形式:

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

式中:n——迭代次數(shù);tn——第n次迭代的時(shí)間步長(zhǎng)。

為簡(jiǎn)潔起見,記

則掩膜主動(dòng)輪廓模型的迭代方程可以簡(jiǎn)寫為如下形式:

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

2 案例分析

2.1 試驗(yàn)數(shù)據(jù)

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

表1 研究使用的遙感數(shù)據(jù)

2.2 模型參數(shù)設(shè)置

在模型迭代過(guò)程中,魯棒對(duì)齊項(xiàng)實(shí)際上對(duì)應(yīng)著圖像中邊緣或梯度所起到的作用;最小方差項(xiàng)對(duì)應(yīng)著對(duì)象內(nèi)部同質(zhì)性大小;測(cè)地主動(dòng)輪廓項(xiàng)對(duì)應(yīng)著曲線的光滑程度;加權(quán)面積項(xiàng)或氣球力起到的主要是加快曲線收縮或膨脹的作用。在這幾個(gè)參數(shù)中,最重要的是控制邊緣的魯棒對(duì)齊項(xiàng)和控制對(duì)象內(nèi)部同質(zhì)性的最小方差項(xiàng)在所有項(xiàng)中所占的權(quán)重,在確定模型參數(shù)時(shí),可根據(jù)具體圖像的不同特征,給兩者賦予適當(dāng)?shù)臋?quán)重系數(shù)。

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

2.3 試驗(yàn)結(jié)果及分析

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

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

表2 模型法面積與全色面積對(duì)比

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

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

2.4 模型應(yīng)用

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

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

3 結(jié)論

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

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

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

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

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

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

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

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

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

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

[9]韓棟,楊曉梅,紀(jì)凱.小衛(wèi)星遙感影像自動(dòng)提取水體方法研究[J].測(cè)繪科學(xué),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.