張永梅，潘振寬，曹叢華，端金鳴，逯京格

(1.國(guó)家海洋局北海預(yù)報(bào)中心，山東 青島 266061；2.山東省海洋生態(tài)環(huán)境與防災(zāi)減災(zāi)重點(diǎn)實(shí)驗(yàn)室, 山東 青島 266061；3.青島大學(xué) 計(jì)算機(jī)科學(xué)技術(shù)學(xué)院，山東 青島 266071)

基于變分水平集方法的滸苔綠潮面積信息提取

張永梅1,2，潘振寬3，曹叢華1,2，端金鳴3，逯京格3

(1.國(guó)家海洋局北海預(yù)報(bào)中心，山東 青島 266061；2.山東省海洋生態(tài)環(huán)境與防災(zāi)減災(zāi)重點(diǎn)實(shí)驗(yàn)室, 山東 青島 266061；3.青島大學(xué) 計(jì)算機(jī)科學(xué)技術(shù)學(xué)院，山東 青島 266071)

綠潮面積信息提取是綠潮遙感監(jiān)測(cè)中極其重要的環(huán)節(jié)，其結(jié)果將直接影響后續(xù)的統(tǒng)計(jì)分析和預(yù)測(cè)預(yù)警工作。目前科研人員一般基于傳統(tǒng)閾值方法對(duì)綠潮面積信息進(jìn)行提取，其提取結(jié)果具有不穩(wěn)定、效率低、人為因素影響大等缺點(diǎn)。針對(duì)上述問(wèn)題，本文基于變分水平集的對(duì)偶方法和分裂Bregman投影方法對(duì)綠潮面積信息進(jìn)行了提取，并提出一種對(duì)上述兩種方法提取出的綠潮面積信息結(jié)果進(jìn)行量化的新方法。分別基于傳統(tǒng)閾值方法、變分水平集的對(duì)偶方法和分裂Bregman投影方法進(jìn)行了3幅影像的綠潮信息提取實(shí)驗(yàn)，對(duì)提取結(jié)果進(jìn)行了比對(duì)分析。實(shí)驗(yàn)結(jié)果表明，對(duì)分辨率較高的衛(wèi)星遙感數(shù)據(jù)，無(wú)論從運(yùn)行效率還是從綠潮面積信息提取結(jié)果的精確性及穩(wěn)定性上，基于變分水平集的對(duì)偶方法和分裂Bregman投影方法均優(yōu)于基于傳統(tǒng)閾值方法。

綠潮；信息提?。粓D像分割；變分水平集方法

1 引言

滸苔是綠藻門石莼科的一屬，其在海上的大面積聚集被稱之為滸苔綠潮，簡(jiǎn)稱綠潮。大量滸苔漂浮聚集到岸邊，會(huì)阻塞航道，破壞海洋生態(tài)系統(tǒng)，其所覆蓋之處會(huì)造成海洋動(dòng)、植物的大量死亡，給近岸養(yǎng)殖戶帶來(lái)巨大的財(cái)產(chǎn)損失。綠潮也會(huì)對(duì)海邊旅游、度假、海上帆船比賽等造成嚴(yán)重影響，直接影響城市的旅游形象及政府的財(cái)政收入。國(guó)內(nèi)外海洋部門已經(jīng)把綠潮定義為一種海洋災(zāi)害。自2008年以來(lái)，每年春夏之間發(fā)生在黃海之濱的綠潮災(zāi)害受到了政府和相關(guān)部門的高度重視，政府和涉海部門都積極應(yīng)對(duì)，力爭(zhēng)把綠潮災(zāi)害造成的影響降到了最低。衛(wèi)星遙感綠潮監(jiān)測(cè)、預(yù)報(bào)可為政府及有關(guān)部門積極應(yīng)對(duì)綠潮災(zāi)害提供第一手材料，而綠潮預(yù)報(bào)的關(guān)鍵在于精確并快速地對(duì)綠潮面積(即綠潮實(shí)際覆蓋面積，下同)信息進(jìn)行提取。綠潮面積信息提取是衛(wèi)星遙感綠潮監(jiān)測(cè)中的重要環(huán)節(jié)，其結(jié)果將直接影響后期的統(tǒng)計(jì)分析及預(yù)測(cè)預(yù)警等工作。目前國(guó)內(nèi)在綠潮面積信息提取時(shí)，大多采用目視判讀、人工設(shè)定閾值的方法來(lái)獲得。這種方法主要依賴于專家知識(shí)，其優(yōu)點(diǎn)是操作簡(jiǎn)單，但具有效率低下、人為影響大等缺點(diǎn)，尤其在綠潮應(yīng)急期間其弊端更加明顯。如梁剛通過(guò)設(shè)定閾值的方式進(jìn)行滸苔信息的提取[1]；鐘山等通過(guò)選擇兩個(gè)閾值解決MODIS影像在綠潮監(jiān)測(cè)中提取出的綠潮實(shí)際面積的較大誤差問(wèn)題[2]；孫立娥對(duì)多源遙感數(shù)據(jù)通過(guò)設(shè)置合適閾值的方法進(jìn)行綠潮信息提取[3]；顧行發(fā)等通過(guò)人工輸入合適閾值的方法提取滸苔分布信息[4]。

綠潮面積信息提取時(shí)只需要關(guān)注感興趣區(qū)(即綠潮)的信息，以海水為背景、綠潮區(qū)為前景進(jìn)行提取，在數(shù)字圖像處理領(lǐng)域即為圖像分割?；谧兎址椒?、PDE方法和水平集方法的變分水平集模型，因其具有自適應(yīng)表達(dá)的拓?fù)浣Y(jié)構(gòu)、穩(wěn)定的計(jì)算方法、自動(dòng)的分割過(guò)程、多模型的集成能力、二/三維圖像分割表達(dá)的一致性等特點(diǎn)，已成為當(dāng)前圖像分割的熱點(diǎn)方法之一[5]。

目前，基于變分模型的圖像分割在很多領(lǐng)域取得了成功應(yīng)用，如物體識(shí)別、醫(yī)學(xué)影像的判讀、人體器官的三維重建、人臉識(shí)別、城市道路規(guī)劃、海洋溢油信息提取等，但還沒(méi)有將圖像分割的變分模型應(yīng)用在衛(wèi)星遙感綠潮監(jiān)測(cè)的信息提取工作中的研究報(bào)道。我們雖然基于變分模型在醫(yī)學(xué)圖像的分割與人體器官三維重建[9]等方面取得了成功應(yīng)用，但還沒(méi)有研究滸苔綠潮信息的提取，經(jīng)過(guò)預(yù)處理工作后的遙感影像和醫(yī)學(xué)影像一樣，只是一張黑白圖像。本文參照提取人體器官信息的原理，試圖把變分模型應(yīng)用到滸苔綠潮面積信息的提取方面。

Kass等是活動(dòng)輪廓的創(chuàng)始人，提出參數(shù)化Snake模型[6]，并在物體識(shí)別等領(lǐng)域做了應(yīng)用研究；1988年Oser和Sethian提出水平集方法[7]并被廣泛應(yīng)用于圖像分割，該方法在物體識(shí)別等領(lǐng)域做了應(yīng)用研究；2001年Chan和Vese提出基于變分水平集方法的兩相Chan-Vese模型[8]，該模型在含噪聲圖像下形狀識(shí)別、星空下雷區(qū)檢測(cè)以及星空下光源檢測(cè)等領(lǐng)域取得良好應(yīng)用。潘振寬等提出了一種新的變分水平集方法[9]，該方法面向三維圖像的多相分割，基于該方法在醫(yī)學(xué)領(lǐng)域做了大量數(shù)值實(shí)驗(yàn)，在醫(yī)學(xué)圖像分割(脂肪、骨骼、肌肉組織、下頜骨、牙齒信息的提取)、人體器官三維重建等方面取得良好的效果；王相海和李明提出了一種基于雙重輪廓演化曲線的圖像分割水平集模型[10]，在物體識(shí)別方面分割結(jié)果魯棒、準(zhǔn)確和穩(wěn)定；程相康等提出了一種快速水平集圖像分割算法，該算法在執(zhí)行效率和分割效果上均優(yōu)于傳統(tǒng)方法，并在醫(yī)學(xué)影像、物體識(shí)別等領(lǐng)域進(jìn)行了實(shí)驗(yàn)[11]。

針對(duì)綠潮災(zāi)害的特點(diǎn)和應(yīng)急減災(zāi)業(yè)務(wù)化需要，本文提出采用基于變分水平集模型的對(duì)偶方法(DM)和分裂Bregman投影方法(SBPM)進(jìn)行衛(wèi)星遙感影像綠潮信息提取，為使上述方法能夠成功應(yīng)用，需要解決的科學(xué)問(wèn)題是模型能否應(yīng)用到業(yè)務(wù)化工作中，模型能否自動(dòng)提取并得到精確的面積數(shù)據(jù)?；谔岢龅膯?wèn)題，首先研究了遙感數(shù)據(jù)的預(yù)處理方法，尤其是不規(guī)則影像裁剪技術(shù)的實(shí)現(xiàn)，標(biāo)志著上述方法能夠應(yīng)用到業(yè)務(wù)化工作中；其次提出一種對(duì)變分模型提取出的綠潮信息結(jié)果進(jìn)行量化的方法，得到了精確的面積數(shù)據(jù)，同時(shí)改變了傳統(tǒng)的人工目視比對(duì)多張圖像分割結(jié)果的弊端。與傳統(tǒng)閾值方法綠潮提取結(jié)果相比較，上述兩種方法無(wú)論從運(yùn)行效率還是從綠潮信息提取結(jié)果的精確性及穩(wěn)定性上均優(yōu)于傳統(tǒng)閾值方法。

2 遙感數(shù)據(jù)預(yù)處理

分別選用了2014年5月26日(實(shí)驗(yàn)區(qū)域一)、2015年5月26日(實(shí)驗(yàn)區(qū)域二)和2016年5月16日(實(shí)驗(yàn)區(qū)域三)3景國(guó)產(chǎn)環(huán)境一號(hào)(HJ-1)衛(wèi)星的數(shù)據(jù)進(jìn)行了分析研究。HJ-1衛(wèi)星是2008年9月發(fā)射成功的，現(xiàn)在雖然處在延壽期，但仍在不斷傳回?cái)?shù)據(jù)，其影像與正常壽命期影像各項(xiàng)指標(biāo)相接近，完全可以與正常壽命期影像一樣，在環(huán)境監(jiān)測(cè)等各方面應(yīng)用[12]。首先對(duì)環(huán)境衛(wèi)星數(shù)據(jù)進(jìn)行圖像鑲嵌、圖像去噪、大氣校正、不規(guī)則區(qū)域圖像裁剪等預(yù)處理工作。其中，不規(guī)則區(qū)域圖像裁剪是變分模型能夠成功應(yīng)用到綠潮提取的業(yè)務(wù)化工作、并能提高綠潮信息提取工作效率的重要工作。

在衛(wèi)星掃描地面時(shí)，一般難以穿透大氣，生成的衛(wèi)星遙感圖片(簡(jiǎn)稱“衛(wèi)片”，下同)含有大氣和光照等信息，為更好地提取綠潮信息，應(yīng)去除大氣和光照的影響，即進(jìn)行大氣校正工作。經(jīng)過(guò)數(shù)據(jù)定標(biāo)、波譜響應(yīng)函數(shù)的制作后利用ENVI平臺(tái)的FLAASH模塊進(jìn)行大氣校正。該模塊內(nèi)嵌目前精度較高的MODTRAN 4+模型，校正精度為像素級(jí)。

本文從NDVI值變化的情況驗(yàn)證大氣校正效果。依據(jù)滸苔的生物特性和光譜特性，常采用計(jì)算圖像的NDVI指數(shù)的方法進(jìn)行滸苔綠潮區(qū)域的判讀。植被指數(shù)有NDVI(歸一化植被指數(shù))、GVI(綠度植被指數(shù))、PVI(垂直植被指數(shù))、SAVI(土壤調(diào)節(jié)植被指數(shù))、EVI(增強(qiáng)型植被指數(shù))、DVI(差值植被指數(shù))及RVI(比值植被指數(shù))等多種。NDVI指數(shù)是最佳的表征陸地植被覆蓋度即生長(zhǎng)狀態(tài)的指示因子，在綠色植被遙感的判讀、解譯中應(yīng)用最廣泛[13]。NDVI值無(wú)單位，其值越高表示像元中綠色植被越多，即滸苔越密集，其值為負(fù)表示像元中滸苔較少。選取任意位置海水、陸地及滸苔區(qū)域像元的大氣校正前后NDVI之變化情況進(jìn)行分析，結(jié)果見(jiàn)表1。

表1 大氣校正前后像元NDVI值的變化情況

從表1可看出大氣校正后相同位置的NDVI值均有較大提高，NDVI的差值定義如下：

NDVIdiff=NDVIenter-NDVIsea，

(1)

式中，NDVIdiff為NDVI的差值結(jié)果，NDVIenter為滸苔的NDVI值，NDVIsea為海水的NDVI值，大氣校正前4個(gè)滸苔像元與第一個(gè)海水像元的差值NDVIdiff分別為0.383 3、0.476 157、0.728 037、0.969 664；大氣校正后4個(gè)滸苔像元與第一個(gè)海水像元的差值NDVIdiff值分別為0.513 738、0.577 918、0.741 091、0.872 315，除最后一個(gè)像元校正后差值比校正前稍小外，其他3個(gè)像元在大氣校正后差值均有明顯的增大。由此可見(jiàn)，大氣校正后增大了大多數(shù)滸苔NDVI值和海水NDVI值的差值，對(duì)滸苔密集度不是很高的像元影響明顯，但對(duì)滸苔密集度高的像元效果不明顯，有利于區(qū)分大部分滸苔信息和海水信息，有利于進(jìn)一步的傳統(tǒng)閾值方法和基于變分方法的綠潮面積信息提取。

一幅遙感影像通常包含陸地、云等非綠潮區(qū)。陸地與滸苔的歸一化植被指數(shù)(NDVI指數(shù))計(jì)算后結(jié)果相近，陸地像元將參與分割；云會(huì)掩蓋綠潮信息，這些因素都會(huì)造成對(duì)綠潮的提取結(jié)果不正確，需要剔除。另外，若圖像信息復(fù)雜則需要采用三相及以上變分模型進(jìn)行圖像分割，計(jì)算效率將受到影響，尤其在綠潮應(yīng)急期間。因此，對(duì)圖像進(jìn)行不規(guī)則區(qū)域裁剪是一項(xiàng)重要的預(yù)處理工作。

不規(guī)則綠潮區(qū)域裁剪可基于ENVI軟件的Subset Data via ROIS工具通過(guò)感興趣區(qū)構(gòu)建數(shù)據(jù)子集的方法實(shí)現(xiàn)，結(jié)果見(jiàn)圖1和圖2。

圖1 實(shí)驗(yàn)區(qū)域位置示意圖Fig.1 Location map of the experimental areas

圖2 最終不規(guī)則圖像裁切結(jié)果圖(實(shí)驗(yàn)區(qū)域一)Fig.2 The final cutting result of irregular image(area 1)

圖1為3景HJ-1A影像實(shí)驗(yàn)區(qū)域的位置示意圖。圖2為實(shí)驗(yàn)區(qū)域一的不規(guī)則圖像裁剪結(jié)果，圖2中的白色為綠潮信息，黑色為海水。經(jīng)過(guò)預(yù)處理工作后，我們得到了一幅質(zhì)量更好、更實(shí)用的遙感影像，尤其是不規(guī)則綠潮區(qū)域裁剪工作，既能夠加快綠潮信息提取速度，又是變分水平集方法能夠進(jìn)行綠潮信息提取業(yè)務(wù)化工作的關(guān)鍵之處。

3 兩相圖像分割的Chan-Vese模型及其快速算法

遙感影像經(jīng)過(guò)預(yù)處理工作后，只有綠潮和海水信息；綠潮信息解譯業(yè)務(wù)化工作對(duì)速度要求較高，尤其是綠潮災(zāi)害應(yīng)急期間。因此，本文基于效率更高的兩相圖像分割的變分水平集方法進(jìn)行綠潮信息提取。

水平集方法和變分水平集方法都可以用來(lái)進(jìn)行遙感影像綠潮信息提取，但因水平集方法分割弱邊緣圖像效果差、對(duì)圖像噪聲敏感及計(jì)算效率低等弊端，本文選用分割效果更好的基于變分水平集的兩相Chan-Vese模型及兩個(gè)快速算法進(jìn)行綠潮信息提取。

3.1 兩相圖像分割的Chan-Vese模型

Chan-Vese模型是基于區(qū)域的著名模型，該模型基于Mumford-Shah模型,但Mumford-Shah模型只是一個(gè)理論模型，因圖像和輪廓線的維度不一致，該模型不可解。為使其可解，研究人員提出改進(jìn)方法，其中，Chan-Vese提出了兩相圖相分割的Chan-Vese模型[8]，該模型引入水平集函數(shù)φ的Heaviside函數(shù)和Dirac函數(shù)，定義如下：

(2)

式中，f表示輸入的含噪聲圖像，u=(u1,u2)，u1、u2分別表示分段常值圖像的前景區(qū)域和背景區(qū)域的圖像灰度均值，|H(φ)|表示長(zhǎng)度項(xiàng)。該模型為分段常值近似模型，前兩項(xiàng)表示分段常值圖像與實(shí)際圖像之間的偏離程度，當(dāng)活動(dòng)輪廓線位于目標(biāo)曲線的邊界時(shí)，前兩項(xiàng)的和為最小。第三項(xiàng)為長(zhǎng)度項(xiàng)，表示活動(dòng)輪廓線的長(zhǎng)度，該值越小，輪廓線越光滑。當(dāng)式(2)能量泛函取得最小值時(shí)，即得到目標(biāo)曲線的理想輪廓，即：

(3)

(4)

公式(3)即為基于變分水平集的兩相圖相分割的Chan-Vese模型，公式(4)為其約束條件。對(duì)公式(3)進(jìn)行求解，首先可以利用交替優(yōu)化方法分別優(yōu)化φ、u1、u2，即當(dāng)φ固定時(shí)，

當(dāng)u1、u2固定時(shí)，

(5)

令R(u1,u2)=α1(u1-f)2-α2(u2-f)2， 對(duì)式(5)求解，求得其對(duì)應(yīng)的歐拉-拉格朗日方程，得到梯度下降流

(6)

上式可通過(guò)有限差分方法近似求解。

3.2 兩種快速算法

基于變分水平集的兩相圖像分割的Chan-Vese模型在曲線演化過(guò)程中為保持符號(hào)距離函數(shù)的特征定義一約束公式：

|φ(x,y,t)|=1，

(7)

為滿足約束條件，研究人員研究出多種快速方法，如分裂Bregman投影算法、對(duì)偶方法等。

3.2.1 對(duì)偶方法

對(duì)偶方法(DualMethod，DM)的主要思想是利用TV范數(shù)的對(duì)偶公式，改變?cè)瓉?lái)能量泛函的極小化處理方式，對(duì)偶公式定義為：

(8)

(9)

式中，R(u1,u2)=α1(u1-f)2-α2(u2-f)2。

3.2.2 分裂Bregman投影算法(SplitBregmanProjectionMethod，SBPM)

(10)

(11)

(12)

(13)

公式(12)可采用半隱式高斯賽德?tīng)?Semi-implicitGaussSeideliteration)迭代方法進(jìn)行求解,得到

(14)

式中，α=φi+1,j+φi-1,j+φi,j+1+φi,j-1。對(duì)公式(14)

進(jìn)行凸松弛，得到

(15)

為提高計(jì)算效率，式(13)上面的方程可采用解析廣義軟閾值公式進(jìn)行求解，得到

(16)

式(13)下面的約束公式可采用投影公式進(jìn)行計(jì)算，即

(17)

ε為正的小量，以防止分母為零溢出。

4 數(shù)值實(shí)驗(yàn)及分析

分別基于傳統(tǒng)閾值方法和變分水平集方法進(jìn)行綠潮信息提取實(shí)驗(yàn)。

遙感數(shù)據(jù)預(yù)處理及傳統(tǒng)閾值方法提取綠潮信息運(yùn)行在美國(guó)ESRI公司的ENVI4.5軟件平臺(tái)上。經(jīng)過(guò)對(duì)圖像進(jìn)行NDVI指數(shù)計(jì)算、對(duì)不規(guī)則綠潮區(qū)域進(jìn)行掩膜、NDVI結(jié)果與掩膜結(jié)果進(jìn)行波段計(jì)算等方法進(jìn)行綠潮信息提取，最終的提取結(jié)果通過(guò)ROI工具輸出保存；基于變分水平集方法的數(shù)值實(shí)驗(yàn)基于Matlab 2010b，將DM和SBPM方法通過(guò)Matlab語(yǔ)言編程實(shí)現(xiàn)，并設(shè)計(jì)一量化公式(詳見(jiàn)公式(18))對(duì)圖像結(jié)果進(jìn)行量化，提取出綠潮精確的面積數(shù)據(jù)，輸入的是不規(guī)則綠潮區(qū)域圖像文件，輸出的是最終提取出的綠潮面積數(shù)據(jù)和運(yùn)行時(shí)間。硬件環(huán)境是Intel(R)Xeon(R)CPU E5-2620 v2@2.10 GHz(雙處理器)，32 GB內(nèi)存的惠普工作站。經(jīng)過(guò)遙感影像預(yù)處理后得到的3個(gè)實(shí)驗(yàn)區(qū)域的圖像大小分別為1 823×1 124、1 109×795、772×620。

4.1 傳統(tǒng)閾值方法

經(jīng)過(guò)人工目視判讀，選取部分閾值進(jìn)行綠潮信息提取，提取部分結(jié)果見(jiàn)表2，表2和表3中的“時(shí)間”均指從開(kāi)始進(jìn)行綠潮信息提取到得到精確的提取面積的工作時(shí)間，不包含前期的遙感數(shù)據(jù)預(yù)處理時(shí)間。選取表2中部分圖像結(jié)果進(jìn)行展示，見(jiàn)圖3至圖8。

表2 傳統(tǒng)閾值方法綠潮信息提取結(jié)果

續(xù)表2

圖3至圖8均為圖像NDVI指數(shù)計(jì)算結(jié)果與基于傳統(tǒng)閾值方法提取出的綠潮信息結(jié)果的疊加，只對(duì)選定的不規(guī)則區(qū)域進(jìn)行提取。為了更清楚地展示提取結(jié)果，對(duì)原NDVI結(jié)果圖像進(jìn)行了反相處理，并在每張圖上對(duì)小紅框表示的局部區(qū)域進(jìn)行了右側(cè)大紅框的3倍放大顯示。圖中紅色為提取出的綠潮信息結(jié)果，黑色為綠潮原信息。閾值的高低表示像元中所含滸苔密集度的高低，實(shí)驗(yàn)區(qū)域一和實(shí)驗(yàn)區(qū)域二的影像均為滸苔生長(zhǎng)茂盛期數(shù)據(jù)，其滸苔密集度高，閾值的選取較大；實(shí)驗(yàn)區(qū)域三影像為滸苔綠潮初生階段數(shù)據(jù)，其滸苔密集度較低，閾值的選取較小。經(jīng)過(guò)人工觀察所有圖像結(jié)果并分析提取出的面積數(shù)據(jù)(即表2)可得出，3個(gè)實(shí)驗(yàn)區(qū)域均隨著閾值的不斷減小，提取出的綠潮信息越來(lái)越精細(xì)，得到的綠潮面積也越來(lái)越大。但并不是閾值可以無(wú)限減小，閾值過(guò)小會(huì)導(dǎo)致微弱的噪聲被提取出來(lái)，閾值的多少，需要人工判讀、設(shè)定。在實(shí)驗(yàn)區(qū)域一中當(dāng)閾值等于-0.07至-0.15時(shí)有較多的滸苔信息沒(méi)有被提取出來(lái)，尤其是當(dāng)閾值為-0.07時(shí)，而當(dāng)閾值等于-0.16、-0.17時(shí)提取出的綠潮信息較多，人工目視判讀也較精確，其中閾值等于-0.17、綠潮面積為51.067 8 km2時(shí)提取的效果最好，基本上綠潮信息都被提取出來(lái)，并且沒(méi)有對(duì)噪聲信息進(jìn)行誤提取。實(shí)驗(yàn)區(qū)域二當(dāng)閾值為-0.2，綠潮面積為18.09 km2時(shí)提取的綠潮信息最好，而當(dāng)閾值為-0.21時(shí)，出現(xiàn)噪聲；實(shí)驗(yàn)區(qū)域三當(dāng)閾值為-0.37、綠潮面積為5.252 4 km2時(shí)提取出的綠潮信息最好，而當(dāng)閾值為-0.38時(shí)，出現(xiàn)噪聲。當(dāng)出現(xiàn)噪聲時(shí)，提取結(jié)果不可用，則停止提取。

圖3 閾值為-0.07時(shí)，綠潮面積為30.301 2 km2(實(shí)驗(yàn)區(qū)域一)Fig.3 Threshold is -0.07，the area of green tide is 30.301 2 km2 (area 1)

圖4 閾值為-0.17時(shí)，綠潮面積為51.067 8 km2(實(shí)驗(yàn)區(qū)域一)Fig.4 Threshold is -0.17，the area of green tide is 51.067 8 km2 (area 1)

圖5 閾值為-0.13時(shí)，綠潮面積為13.337 1 km2(實(shí)驗(yàn)區(qū)域二)Fig.5 Threshold is -0.13，the area of green tide is 13.337 1 km2 (area 2)

圖6 閾值為-0.2時(shí)，綠潮面積為18.09 km2(實(shí)驗(yàn)區(qū)域二)Fig.6 Threshold is -0.2，the area of green tide is 18.09 km2(area 2)

圖7 閾值為-0.3時(shí)，綠潮面積為1.259 1 km2(實(shí)驗(yàn)區(qū)域三)Fig.7 Threshold is -0.3，the area of green tide is 1.259 1 km2 (area 3)

圖8 閾值為-0.37時(shí)，綠潮面積為5.252 4 km2(實(shí)驗(yàn)區(qū)域三)Fig.8 Threshold is -0.37，the area of green tide is 5.252 4 km2 (area 3)

4.2 變分水平集方法

本文基于對(duì)偶方法(DM)和分裂Bregman投影方法(SBPM)進(jìn)行衛(wèi)星遙感影像綠潮信息提取數(shù)值實(shí)驗(yàn)，為與傳統(tǒng)閾值方法提取綠潮信息進(jìn)行比較，對(duì)基于變分水平集方法的綠潮信息提取結(jié)果成功進(jìn)行量化，得出最終的綠潮提取面積，單位是平方千米。

4.2.1 基于變分水平集方法的綠潮信息提取結(jié)果量化方法

當(dāng)研究區(qū)域是規(guī)則圖像時(shí)，其最終水平集函數(shù)分割結(jié)果是前景像素值是1，背景是0，不規(guī)則圖像的分割結(jié)果圖像內(nèi)與規(guī)則圖像相反，圖像內(nèi)前景像素值是0，背景是1，圖像外的區(qū)域?yàn)?。根據(jù)上述規(guī)律并與ENVI軟件相結(jié)合可設(shè)計(jì)一量化公式得出提取出的綠潮區(qū)域總面積：

Last_area=(clip_pixel-last_sum)/

clip_pixel×clip_area

(18)

式中，Last_area為最終提取出的綠潮信息的總面積，clip_pixel為不規(guī)則圖像的像素點(diǎn)總個(gè)數(shù)，last_sum為不規(guī)則圖像內(nèi)的背景像素點(diǎn)總個(gè)數(shù)，clip_area為不規(guī)則圖像的總面積。

4.2.2 數(shù)值實(shí)驗(yàn)

本文采用一個(gè)水平集函數(shù)，DM方法和SBPM方法部分參數(shù)初始值相同：γ=1，σ=1，θ=3 000，DM方法中的τ=0.125。當(dāng)上一次的綠潮提取面積與本次的提取面積的差值小于某最小值ε時(shí)，程序的循環(huán)次數(shù)即為迭代次數(shù)，ε的單位為平方千米。經(jīng)過(guò)多組試驗(yàn)，并和傳統(tǒng)閾值方法提取結(jié)果比較分析發(fā)現(xiàn)，當(dāng)最小值ε為0.5時(shí)，提取出的綠潮面積結(jié)果不論從程序提取面積的精度上、人工目視查看圖像的結(jié)果，還是從運(yùn)行時(shí)間的效率上都能達(dá)到最優(yōu)。最終提取出的綠潮信息部分?jǐn)?shù)據(jù)結(jié)果見(jiàn)表3，部分圖像結(jié)果見(jiàn)圖9至圖16。

表3 DM方法和SBPM方法綠潮信息提取結(jié)果

圖9 SBPM方法初始化(實(shí)驗(yàn)區(qū)域一)Fig.9 The image processing result by the SBPM method initializing(area 1)

圖10 DM方法初始化(實(shí)驗(yàn)區(qū)域一)Fig.10 The image processing result by the DM method initializing(area 1)

圖11 區(qū)域一(DM方法：α=0.1, 迭代次數(shù)為20, 綠潮面積為51.640 9 km2)Fig.11 Area 1 (the DM method: α=0.1, the number of iterations is 20, the area of green tide is 51.640 9 km2)

圖12 區(qū)域一(SBPM方法：α=1, 迭代次數(shù)為7, 綠潮面積為50.412 7 km2)Fig.12 Area 2 (the SBPM method：α=1, the number of iterations is 7, the area of green tide is 50.412 7 km2)

圖13 區(qū)域二(DM方法：α=0.1, 迭代次數(shù)為22, 綠潮面積為18.075 6 km2)Fig.13 Area 2 (the DM method: α=0.1, the number of iterations is 22, the area of green tide is 18.075 6 km2)

圖14 區(qū)域二(SBPM方法：α=1, 迭代次數(shù)為9，綠潮面積為18.167 4 km2)Fig.14 Area 2 (the SBPM method: α=1, the number of iterations is 9,the area of green tide is 18.167 4 km2)

圖15 區(qū)域三(DM方法：α=1, 迭代次數(shù)為5, 綠潮面積為5.299 2 km2)Fig.15 Area 3 (the DM method: α=1, the number of iterations is 5, the area of green tide is 5.299 2 km2)

圖16 區(qū)域三(SBPM方法：α=0.1, 迭代次數(shù)為22，綠潮面積為4.756 5 km2)Fig.16 Areas 3 (the SBPM method: α=0.1, the number of iterations is 22, the area of green tide is 4.756 5 km2)

圖17 局部放大結(jié)果(實(shí)驗(yàn)區(qū)域1——傳統(tǒng)方法：閾值為-0.17時(shí)，綠潮面積為51.067 8 km2)Fig.17 Local amplification resuls(area 1- the traditional method: threshold= -0.17，the area of green tide =51.067 8 km2)

圖18 局部放大結(jié)果(實(shí)驗(yàn)區(qū)域1—變分方法：DM方法α=0.1, 迭代次數(shù)為20，綠潮面積為51.640 9 km2)Fig.18 Local amplification resuls(area 1- variational method: the DM method: α=0.1, the number of iterations is 20, the area of green tide is 51.640 9 km2)

圖9至圖16為基于變分水平集方法的綠潮信息提取的部分結(jié)果，為增強(qiáng)顯示，對(duì)原NDVI結(jié)果進(jìn)行了圖像反相處理，并在每張圖上對(duì)小紅框表示的局部區(qū)域進(jìn)行了右側(cè)大紅框的放大顯示，除圖11、圖12對(duì)局部放大了2倍外，其他局部均放大3倍，藍(lán)色為圖像外區(qū)域，紅色為提取出的綠潮信息結(jié)果，黑色為綠潮原信息。圖17、圖18為基于傳統(tǒng)方法和變分方法提取綠潮面積值基本接近時(shí)圖像的細(xì)節(jié)對(duì)比圖。從表3可以看出，初始參數(shù)不同，提取出的綠潮面積略有不同。從對(duì)傳統(tǒng)閾值方法進(jìn)行綠潮信息提取的分析中可得出實(shí)驗(yàn)區(qū)域一當(dāng)閾值等于-0.17時(shí)綠潮信息提取得最好(提取的綠潮面積是51.067 8 km2)、實(shí)驗(yàn)區(qū)域二當(dāng)閾值等于-0.2時(shí)綠潮信息提取得最好(提取的綠潮面積是18.09 km2)、實(shí)驗(yàn)區(qū)域三當(dāng)閾值等于-0.37時(shí)綠潮信息提取得最好(提取的綠潮面積是5.252 4 km2)，在基于變分的綠潮信息提取實(shí)驗(yàn)中(表3)，實(shí)驗(yàn)區(qū)域一除第5、7、11及12次實(shí)驗(yàn)與閾值等于-0.17時(shí)的結(jié)果有一定的偏差外，其他實(shí)驗(yàn)結(jié)果均與傳統(tǒng)方法閾值為-0.17的綠潮信息提取最好時(shí)的面積接近；實(shí)驗(yàn)區(qū)域二除第10及14次實(shí)驗(yàn)與閾值等于-0.2時(shí)的結(jié)果有一定的偏差外，其他實(shí)驗(yàn)結(jié)果均與傳統(tǒng)方法閾值為-0.2的綠潮信息提取最好時(shí)的面積接近；實(shí)驗(yàn)區(qū)域三除第13及14次實(shí)驗(yàn)與閾值等于-0.37時(shí)的結(jié)果有一定的偏差外，其他實(shí)驗(yàn)結(jié)果均與傳統(tǒng)方法閾值為-0.37的綠潮信息提取最好時(shí)的面積接近?；谧兎址椒ǖ奶崛〗Y(jié)果精度高且穩(wěn)定性好。從表2和表3還可以看出兩者在計(jì)算時(shí)間上差別較大，實(shí)驗(yàn)區(qū)域一基于傳統(tǒng)閾值方法的運(yùn)行時(shí)間最短需要1 min 33 s、最長(zhǎng)需要1 min 59 s，而基于變分水平集的運(yùn)行時(shí)間最短僅為15.366 1 s，最長(zhǎng)也不過(guò)為73.211 3 s(即1 min 13 s)；實(shí)驗(yàn)區(qū)域二基于傳統(tǒng)閾值方法的運(yùn)行時(shí)間最短需要1 min、最長(zhǎng)需要1 min 24 s，而基于變分水平集的運(yùn)行時(shí)間均不超過(guò)1 min，最短僅為7.082 4 s，最長(zhǎng)也不過(guò)為45.957 9 s；實(shí)驗(yàn)區(qū)域三基于傳統(tǒng)閾值方法的運(yùn)行時(shí)間最短需要53 s、最長(zhǎng)需要1 min 12 s，而基于變分水平集的運(yùn)行時(shí)間均不超過(guò)1 min，最短僅為3.229 2 s，最長(zhǎng)也不過(guò)為10.608 1 s，在運(yùn)行時(shí)間上，基于變分的方法明顯優(yōu)于傳統(tǒng)閾值方法。

5 結(jié)論

目前科研人員一般基于傳統(tǒng)閾值方法對(duì)綠潮信息進(jìn)行提取，其提取結(jié)果具有人為因素影響大、不穩(wěn)定、效率低等缺點(diǎn)，針對(duì)這些缺點(diǎn)，本文基于圖像分割的變分水平集方法進(jìn)行了綠潮信息提取研究，重點(diǎn)研究了基于變分水平集的兩相Chan-Vese模型及SBPM方法和DM方法等快速算法，提出了一種對(duì)這兩種方法提取的綠潮信息結(jié)果進(jìn)行量化的方法，并分別進(jìn)行了多景遙感數(shù)據(jù)的傳統(tǒng)閾值方法和基于變分模型的SBPM方法和DM方法提取綠潮信息的數(shù)值實(shí)驗(yàn)。實(shí)驗(yàn)結(jié)果表明，圖像分割的變分模型能夠應(yīng)用到綠潮衛(wèi)星遙感監(jiān)測(cè)業(yè)務(wù)化工作中；對(duì)分辨率較高的衛(wèi)星遙感數(shù)據(jù)，基于變分水平集方法提取滸苔綠潮信息面積提取精度高且穩(wěn)定性好，在時(shí)間效率上有較大的優(yōu)勢(shì)。

本文僅針對(duì)HJ-1A數(shù)據(jù)、基于SBPM方法和DM方法提取綠潮信息，未來(lái)可針對(duì)其他數(shù)據(jù)源、研究其他模型進(jìn)行綠潮信息提取，以提高衛(wèi)星遙感業(yè)務(wù)化監(jiān)測(cè)綠潮、解譯綠潮的水平。

[1] 梁剛. 大型藻類遙感監(jiān)測(cè)方法研究[D]. 大連: 大連海事大學(xué), 2011.

Liang Gang. The method of monitoring macroalgae by remote sensing[D]. Dalian: Dalian Maritime University, 2011.

[2] 鐘山, 丁一, 李振. MODIS滸苔遙感監(jiān)測(cè)誤差分析研究[J]. 遙感信息, 2013, 28(1):38-42.

Zhong Shan, Ding Yi, Li Zhen. Error analysis on enteromorpha prolifera monitoring using MODIS data[J]. Remote Sensing Information, 2013, 28(1):38-42.

[3] 孫立娥. 綠潮遙感探測(cè)影響因素分析[D]. 青島: 青島大學(xué), 2012.

Sun Lie. Analysis on influences of the green tide remote sensing detecting[D]. Qingdao: Qingdao University, 2012.

[4] 顧行發(fā), 陳興峰, 尹球, 等. 黃海滸苔災(zāi)害遙感立體監(jiān)測(cè)[J]. 光譜學(xué)與光譜分析, 2011, 31(6):1627-1632.

Gu Xingfa, Chen Xingfeng, Yin Qiu, et al. Stereoscopic remote sensing used in enteromorpha prolifra disaster in Chinese Yellow Sea[J]. Spectroscopy and Spectral Analysis, 2011, 31(6):1627-1632.

[5] Chan T, Sandberg B, Moelich M. Some recent developments in variational image segmentation[M]//Image Processing Based on Partial Differential Equations, New York: Springer-Verlag, 2007, Part Ⅲ:175-210.

[6] Kass M，Witriw A，Terzopoulos D. Snakes: Active contour models[J]. International Journal of Computer Vision,1988,1:321-331.

[7] Osher S, Sethian J A. Fronts propagating with curvature-dependent speed: algorithms based on hamilton-jacobi formulations[J]. Journal of Computational Physics, 1988, 79(1):12-49.

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

[9] 潘振寬, 李華, 魏偉波, 等. 三維圖像多相分割的變分水平集方法[J]. 計(jì)算機(jī)學(xué)報(bào), 2009, 32(12): 2464-2474.

Pan Zhenkuan, Li Hua, Wei Weibo, et al. A variational level set method of multiphase segmentation for 3D images[J]. Chinese Journal of Computer, 2009, 32(12):2464-2474.

[10] 王相海, 李明. 矢量C-V模型的高光譜遙感影像分割[J]. 遙感學(xué)報(bào), 2015, 19(3):443-450.

Wang Xianghai, Li Ming. Segmentation of hyperspectral remote sensing image using vector C-V model[J]. Journal of Remote Sensing, 2015, 19(3):443-450.

[11] 程相康，朱宏擎. 一種基于高斯混合模型的快速水平集圖像分割方法[J]. 華東理工大學(xué)學(xué)報(bào)(自然科學(xué)版)，2015, 41(6): 808-814.

Cheng Xiangkang, Zhu Hongqing. A fast level set method for image segmentation based on Gaussian Mixture Models[J].Journal of East China University of Science and Technology:Natural Science Edition, 2015，41(6)：808-814.

[12] 蘇曉慧, 張曉東, 蘇偉, 等. HJ-1衛(wèi)星延壽期的CCD影像質(zhì)量評(píng)價(jià)與可用性分析[J]. 農(nóng)業(yè)工程學(xué)報(bào), 2012, 28(12):164-170.

Su Xiaohui, Zhang Xiaodong, Su Wei, et al. Quality evaluation and usability annlysis of CCD image in life extension period for HJ-1 satellite[J]. Transactions of the Chinese Society of Agricultural Engineering, 2012, 28(12):164-170.

[13] 趙英時(shí)，陳冬梅，李小明. 遙感應(yīng)用分析原理與方法[M]. 北京: 科學(xué)出版社, 2003.

Zhao Yingshi, Chen Dongmei, Li Xiaoming. Analysis Principle and Method of Remote Sensing Applications[M]. Beijing: Science Press, 2003.

Information extraction of enteromorpha green tide area based on variational level set method

Zhang Yongmei1,2, Pan Zhenkuan3, Cao Conghua1,2, Duan Jinming3, Lu Jingge3

(1.NorthChinaSeaMarineForecastingCenterofStateOceanicAdministration,Qingdao266061,China;2.ShandongProvincialLaboratoryofMarineEcologyandEnvironment&DisasterPreventionandMitigation,Qingdao266061,China;3.SchoolofComputerScience&Technology,QingdaoUniversity,Qingdao266071,China)

Green tide area information extraction is an important link of the remote sensing monitoring, its result will directly affect the subsequent statistical analysis and the early warning prediction. Now researchers generally extract the green tide area information based on the traditional threshold method, and this approach has many disadvantages such as low efficiency, unstable result and human factors. Against to above-mentioned problem, the dual method and split Bregman projection method based on variational level set method was studied used to the green tide area information extraction in this paper. A new quantization method was proposed, which was used to deal with the green tide information result extracted by the two mentioned methods. Based on the traditional threshold method and the dual method and split Bregman projection method, the experiments of three images were respectively carried out, and extraction results were compared and analyzed. To higher resolution satellite remote sensing data, the experiment results show that not only the extraction efficiency but the accuracy and stability based on variational level set method are all superior to the traditional threshold method.

greed tide; information extraction; image segmentation; variational level set method

10.3969/j.issn.0253-4193.2017.09.012

2016-09-25；

2016-10-30。

國(guó)家自然科學(xué)基金項(xiàng)目(41306028)。

張永梅(1971—),女,山東省淄博市人，高級(jí)工程師,主要研究方向?yàn)檫b感圖像處理、變分圖像分割。E-mail:zhyongmei82@126.com

TP79

A

0253-4193(2017)09-0121-12

張永梅，潘振寬，曹叢華, 等. 基于變分水平集方法的滸苔綠潮面積信息提取[J]. 海洋學(xué)報(bào), 2017, 39(9): 121-132,

Zhang Yongmei, Pan Zhenkuan, Cao Conghua, et al. Information extraction of enteromorpha green tide area based on variational level set method[J]. Haiyang Xuebao, 2017, 39(9): 121-132, doi:10.3969/j.issn.0253-4193.2017.09.012