欒奎峰，劉 帥，潘與佳，朱衛(wèi)東*，李丕學,裘 誠，邱振戈，沈 蔚，3，王 潔，王振華

(1.上海海洋大學 海洋科學學院，上海 201306；2.上海市海洋監(jiān)測預報中心，上海 200062；3.上海河口海洋測繪工程技術研究中心，上海 201306)

0 引言

海岸線是海洋與陸地的分界線[1]。如何快速、準確提取海岸線具有重要意義。遙感技術具有全天候、全天時、覆蓋面積廣和數(shù)據(jù)來源豐富等優(yōu)點，已成為海岸線提取的重要手段[2-3]。現(xiàn)有基于遙感手段的海岸線提取方法主要有：閾值分割法[4]、邊緣檢測法[5]、活動輪廓模型法[6]、面向?qū)ο蠓╗7]、區(qū)域生長法[8]和元胞自動機法[9]等。分水嶺算法是一種基于形態(tài)學的圖像分割算法，定位準確，適用于分割粘連的目標，得到單個像素寬度和連續(xù)封閉分割邊緣，適用于海岸線的提取。VINCENT et al[10]最早提出一種比較經(jīng)典的分水嶺檢測算法，以梯度圖像作為輸入進行分割；GAO et al[11]先采用數(shù)學形態(tài)學中的面積算子對圖像進行簡化，然后再進行分水嶺分割；NING et al[12]在Mean-Shift聚類分割后進行分水嶺分割，并通過人工干預進行區(qū)域合并得到分割結(jié)果；林振榮 等[13]先采用形態(tài)學濾波進行降噪處理，然后在分水嶺分割后進行區(qū)域合并；PUISSANT et al[14]基于形態(tài)學工具“擊中擊不中變換”定義模板，通過雙重閾值對結(jié)果進行組合，再進行標記的分水嶺變換，得出準確的海岸線。隨著國產(chǎn)高分辨率多光譜衛(wèi)星數(shù)據(jù)的發(fā)展，影像上像元的純凈度提高，地物更加清晰，但細節(jié)也被放大，噪聲變得更多，同物異譜和異物同譜現(xiàn)象大量存在，使影像中出現(xiàn)許多局部灰度極值區(qū)域，原有的標記分水嶺算法存在過分割現(xiàn)象，產(chǎn)生許多無效數(shù)據(jù)。

針對分水嶺方法在高分辨率多光譜數(shù)據(jù)提取岸線中存在的過分割和局部極值問題，本文提出了基于擴展極大值和極小值變換的標記分水嶺水邊線提取方法，采用形態(tài)學重建和擴展極值變換等方法抑制噪聲點的影響，以達到較好的水邊線提取效果。

1 改進的擴展極值變換標記分水嶺水邊線提取方法

改進的擴展極值變換標記分水嶺水邊線提取方法主要步驟如圖1所示：首先對待分割的高分辨率多光譜影像進行輻射定標和大氣校正等預處理，其次對預處理后的影像進行歸一化水體指數(shù)(Normalized Difference Water Index，NDWI)計算和梯度圖像計算，并對NDWI計算后的影像進行形態(tài)學重建。然后，使用本文提出的改進的擴展極值變換方法對圖像進行變換，獲得前景和背景標記，并對原梯度圖像進行強制最小值標定處理，修正梯度圖像。最后，對修正后的梯度圖像進行分水嶺變換，完成水邊線的提取。

圖1 算法流程圖Fig.1 Algorithm flow chart

1.1 形態(tài)學重建濾波

分水嶺變換是利用圖像灰度梯度幅值作為分割函數(shù)，對每個像素點的灰度級進行排序，然后按照由低到高的順序進行處理[15]。為突出水體信息、增強水體與陸地部分像元灰度值的差異，采用NDWI計算方法，對高分辨率多光譜影像數(shù)據(jù)進行波段運算。之后，利用形態(tài)學重建濾波方法對處理后的影像進行平滑，去除影像中的噪聲。形態(tài)學基本運算包括腐蝕、膨脹、開運算和閉運算，開運算是先腐蝕后膨脹，閉運算是先膨脹后腐蝕。膨脹和腐蝕運算定義如下

(g⊕l)(x,y)=max{g(x+i,y+j)-
l(i,j)|(x+i,y+j)∈Dg,(i,j)∈Dl}

(1)

(gΘl)(x,y)=min{g(x+i,y+j)-
l(i,j)|(x+i,y+j)∈Dg,(i,j)∈Dl}

(2)

式中：g為待處理圖像，l為結(jié)構(gòu)元素，(x,y)為待處理圖像中運算的原點，(i,j)表示運算過程中結(jié)構(gòu)元素范圍內(nèi)像元的位置，Dg表示待處理圖像的范圍，Dl表示結(jié)構(gòu)元素的范圍，⊕、Θ分別表示膨脹和腐蝕運算。

(3)

(4)

結(jié)構(gòu)元素l的選取影響著影像形態(tài)學重建后的結(jié)果。在常用的結(jié)構(gòu)元素中，圓盤型結(jié)構(gòu)元素對各種地物的處理效果都比較均衡，具有各向同性，不會對圖像特征值造成畸變[17]。結(jié)構(gòu)元素大小m應大于目標區(qū)域內(nèi)部單體地物，同時小于目標區(qū)域的最大內(nèi)切圓半徑，目的是使目標區(qū)域內(nèi)部更為平滑，同時保持其整體的外形輪廓。通過先進行開重建運算再進行閉重建運算的二次重建濾波，可有效平滑圖像和削弱噪聲。

1.2 基于擴展極值變換的形態(tài)學標記

經(jīng)過形態(tài)學重建處理的多光譜影像數(shù)據(jù)，雖然大部分局部噪聲被去除，但仍存在一些對影像梯度修正的準確性造成干擾的噪聲。擴展極值變換[18]是一種形態(tài)學閾值算子，是通過設定閾值，限制灰度局部極值點的數(shù)量，提取出顯著的灰度極值區(qū)域，將大多數(shù)無關區(qū)域標記為0。一般的研究是在影像分割過程中僅采用擴展極小值變換對影像中的灰度極小值區(qū)域進行限制，然后進行標記提取[17,19-20]。但由于高分辨多光譜衛(wèi)星影像細節(jié)豐富，同物異譜和異物同譜等因素，影像中不僅灰度極小值區(qū)域會對分割結(jié)果造成影響，過多的灰度極大值點和區(qū)域也會對分割結(jié)果造成一定程度的干擾。因此，本文提出了基于擴展極小值和極大值變換的標記分水嶺方法，以此對形態(tài)學重建后的圖像進行重構(gòu)，從而解決高分辨多光譜數(shù)據(jù)的過分割問題。

形態(tài)學重建濾波處理后的圖像，陸地范圍的像元灰度值低于水體像元的灰度值，則采用擴展極小值變換標記陸地，其公式為

If=Hmin(M|h1)

(5)

式中：If為擴展極小值變換后的圖像，Hmin為形態(tài)學擴展極小值變換運算，M為形態(tài)學重建后的圖像，h1表示擴展極小值變換運算的灰度閾值。擴展極小值變換后的圖像中陸地區(qū)域標記為1，其余部分標記為0，消除陸地部分低于閾值的局部灰度極小值區(qū)域，并與已標記為陸地的區(qū)域合并。

為去除水體區(qū)域中存在的局部噪聲點，本文采用擴展極大值變換的方法，對圖像進行擴展極值變換，其公式為

Ib=Hmax(M|h2)

(6)

式中：Ib為擴展極大值變換后的圖像，Hmax表示形態(tài)學的擴展極大值變換運算，h2為擴展極大值變換運算的灰度閾值。擴展極大值變換后的圖像中水體區(qū)域標記為1，其余部分標記為0，消除水域部分高于閾值的局部灰度極大值區(qū)域，并與已標記為水域部分的區(qū)域進行合并。

擴展極小值和極大值變換過程中，h1和h2的設定是擴展極值變換的關鍵，閾值過小則處理后的圖像仍會存在過分割現(xiàn)象，閾值過大會使得最終的結(jié)果出現(xiàn)錯誤分割。本文采用基于最大類間方差法[21]設定上述兩個閾值的初值，其公式為

h1=t-Imin

(7)

h2=Imax-t

(8)

式中：t為采用最大類間方差法計算的圖像分割閾值，Imin和Imax分別為全圖像中像元灰度的最小值和最大值。由于高分辨率多光譜數(shù)據(jù)的差異，h1和h2的值在實際計算中需要根據(jù)標記計算的效果再進行細微調(diào)整。

本文對擴展極值變換后的圖像進行形態(tài)學開、閉運算和腐蝕運算的處理，基于前景和背景標記的結(jié)果，采用形態(tài)學強制最小值標定的方法修正原始梯度圖像。修正后的梯度圖像f1為[20]

f1=IMmin(f|I)

(9)

式中：IMmin表示形態(tài)學強制最小值標定運算，f表示原始梯度圖像，I表示前景和背景標記圖像。f1對應于標記圖像I中為1的點，均被標定為0，其余像元點的灰度值與原始梯度圖像保持不變。修正后的梯度圖像中保留邊緣輪廓信息，其余區(qū)域均標記為0，完成局部噪聲的去除。在此基礎上再開展分水嶺變換的處理，得到單像素寬的準確水邊線[22]。

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

2.1 研究區(qū)概況及數(shù)據(jù)源

本文選取南海海域的南薰島和鴻庥島作為2種岸線類型的研究區(qū)(圖2)。南薰島是吹沙填海而來的人工島，島嶼面積約為0.18 km2，岸線類型主要為人工岸線。鴻庥島為自然形成的小島嶼，面積約為0.08 km2，外觀呈橢圓形，岸線類型主要為沙質(zhì)岸線。兩個島嶼海岸線類型特征明顯。

圖2 研究區(qū)地理位置示意圖Fig.2 Geographical position of the study area

遙感數(shù)據(jù)選用GF-2衛(wèi)星遙感影像，2個島嶼的衛(wèi)星影像獲取時間分別為2017年5月29日和2017年3月21日。高分二號全色影像空間分辨率為1 m，多光譜影像空間分辨率為4 m。其中全色影像波段為450～900 nm；多光譜影像分為4個波段，分別是450～520 nm(藍波段)、520～590 nm(綠波段)、630～690 nm(紅波段)、770～890 nm(近紅外波段)[23]。本文使用其多光譜影像數(shù)據(jù)進行實驗。

2.2 實驗結(jié)果與分析

2.2.1 數(shù)據(jù)處理

實驗首先對南薰島原始遙感圖像進行輻射定標、大氣校正等預處理(圖3a)，再進行NDWI運算，并建立Sobel梯度圖像(圖3b)，圖像邊緣信息得到了加強和突出，但由于存在過多的噪聲點，會出現(xiàn)嚴重的過分割現(xiàn)象(圖3c)。在NDWI運算后進行形態(tài)學重建濾波(圖3d)，圖像中海陸像元灰度值對比差異顯著，噪聲在一定程度上得到消除，同時邊緣輪廓信息也能得到較好的保留。其次，針對陸地部分的噪聲，在形態(tài)學重建后進行擴展極小值變換運算，為保證能得到分割結(jié)果，再進行形態(tài)學開、閉運算和腐蝕運算來收縮邊緣，得到陸地標記結(jié)果(圖3e)，圖像中陸地部分的灰度極小值點被消除，但水體部分由于存在過多的灰度極大值點，會對分割結(jié)果造成干擾。針對水體部分存在過多的灰度極大值點，采用本文提出的擴展極大值變換運算方法，在形態(tài)學重建后進行擴展極大值變換運算，同時收縮邊緣，得到水體標記結(jié)果(圖3f)，可發(fā)現(xiàn)圖像中水體部分的灰度極大值點和區(qū)域得到有效的去除，水體部分過分割現(xiàn)象基本消除。結(jié)合前景(陸地)、背景(水體)標記，采用形態(tài)學強制最小值標定方法對圖3b進行修正，得到修正后梯度圖像(圖3g)，并在此基礎上進行分水嶺變換，圖3h為分割結(jié)果，中間單像素寬的分割線即為水邊線，疊加在原圖上如圖3i所示。

(a)預處理后圖像

(d)形態(tài)學重建圖像

(g)修正后梯度圖像

對比本文方法提取的水邊線和圖像中的水邊線，輪廓封閉，形狀較為規(guī)則，對于伸向海洋的碼頭邊緣定位也較準確，兩者基本重合，說明本文所用算法對南薰島的陸地與海水部分能夠較為準確地分割與提取。

同樣，用本文算法對以沙質(zhì)岸線為主的鴻庥島進行實驗，實驗結(jié)果如圖4所示。在對預處理后圖像進行NDWI運算和形態(tài)學重建后，圖像上的大部分噪聲點已經(jīng)消除，但依舊存在部分灰度極值區(qū)域?qū)τ跋裉荻刃拚臏蚀_性造成干擾。采用擴展極小值變換對陸地部分的灰度極小值區(qū)域進行限制，同時采用擴展極大值變換對水體部分的灰度極大值區(qū)域進行限制，完成了前景和背景標記，圖像中陸地和水體內(nèi)部的灰度極值區(qū)域分別得到有效的去除。然后，采用形態(tài)學強制最小值標定方法和分水嶺變換后得到單像素寬的水邊線。從結(jié)果中發(fā)現(xiàn)，沙質(zhì)岸線同樣可以得到較為準確的海陸分割和水邊線提取結(jié)果，對沙灘能夠進行準確識別，水邊線提取結(jié)果與圖像中的水邊線基本吻合。

(a)預處理后圖像

(c)過分割圖像

(e)擴展極小值變換后陸地標記

(g)修正后梯度圖像

(h)分割結(jié)果疊加圖像

擴展極值變換運算中，對于h1和h2的設定是變換的關鍵所在。根據(jù)公式(7)和(8)的閾值確定方法，南薰島影像的h1和h2值分別為55和60。在上述閾值基礎上，分別設定過小和過大的閾值進行實驗，實驗結(jié)果如圖5所示。從結(jié)果中發(fā)現(xiàn)，當h1和h2均設定過小時(圖5a)，圖像中仍會存在過多的灰度極值區(qū)域，最終結(jié)果會出現(xiàn)過分割現(xiàn)象；當h1設定過小、h2設定過大(圖5b)或h1設定過大、h2設定過小(圖5c)時，會使得圖像中的灰度極值區(qū)域被過分消除或合并，最終均會出現(xiàn)錯誤分割的結(jié)果；當h1和h2均設定過大時，會使得圖像中的灰度極值區(qū)域被完全消除，最終無法得到分割結(jié)果，使整幅影像呈現(xiàn)單一灰度值。由此說明擴展極值變換閾值設定得過大或過小均會使得最終結(jié)果出現(xiàn)過分割、錯誤分割現(xiàn)象，甚至無法得到分割結(jié)果。

在以上實驗過程中，由于部分水陸邊緣光譜信息較復雜，水陸差異對比不顯著，導致分水嶺變換時部分海陸邊界定位不夠準確，對水邊線的提取精度存在一定影響。

2.2.2 精度評價

為定量評價本文方法提取水邊線的精度和可靠性，以目視解譯提取的水邊線作為參考，采用緩沖區(qū)分析法[24]開展精度驗證。該方法中，先對參考水邊線以一定的像素寬度建立緩沖區(qū)，提取出的水邊線如果落入到該緩沖區(qū)內(nèi)，則認為該段水邊線提取結(jié)果與參考水邊線相匹配，其長度記為TP1；如若不在緩沖區(qū)內(nèi)，則認為不相匹配，長度記為FP。對采用本文方法提取出的水邊線也建立相同像素寬度的緩沖區(qū)，如果參考水邊線落入到緩沖區(qū)內(nèi)，則認為二者相匹配，記其長度為TP2；如果不在緩沖區(qū)內(nèi)，記其長度為FN?；谏鲜龇治?，建立完整度(Complete)、正確度(Correct)和質(zhì)量(Quality)這三個評價指標：

Complete=TP2/(TP2+FN)

(10)

Correct=TP1/(TP1+FP)

(11)

(12)

長度相對誤差也常用于評價海岸線提取精度，其計算方法為：(提取的海岸線長度—參考海岸線長度)/參考海岸線長度。該值越低，說明長度誤差越小，精度越高。

通過以上4個評價指標驗證本方法提取水邊線的精度，其結(jié)果如表1和表2所示。

表1 南薰島提取結(jié)果精度分析Tab.1 Accuracy analysis of Nanxun Dao extraction results

表2 鴻庥島提取結(jié)果精度分析Tab.2 Accuracy analysis of Hongxiu Dao extraction results

從上述結(jié)果可以發(fā)現(xiàn)，南薰島的長度相對誤差為0.57%，鴻庥島的長度相對誤差為-1.09%，說明兩個島嶼水邊線提取結(jié)果的長度精度均較高。對于南薰島，當緩沖半徑為0.5個像素時，完整度為92.08%，正確度為91.55%，質(zhì)量為84.87%，提取效果良好；當緩沖半徑為1個像素時，完整度為96.67%，正確度為96.00%，質(zhì)量為92.93%，提取的準確性較好，定位精度較高，說明該方法對于人工岸線的提取精度能夠達到1個像素。對于鴻庥島，當緩沖半徑為0.5個像素時，質(zhì)量僅有69.09%，提取效果較差；當緩沖半徑為1個像素時，正確度達到90%以上，說明提取結(jié)果較為準確；當緩沖半徑為1.5個像素時，完整度、正確度和質(zhì)量均在90%以上，此時提取效果較好，說明該方法提取沙質(zhì)岸線的精度能達到1.5個像素。

3 結(jié)論

針對分水嶺變換方法在高分辨率多光譜衛(wèi)星影像數(shù)據(jù)應用中存在同物異譜和異物同譜等因素引起過分割的問題，本文提出了一種基于擴展極值變換標記分水嶺的海島水邊線提取算法，并采用高分二號衛(wèi)星多光譜影像，在南海海域的南薰島和鴻庥島進行了水邊線提取的方法驗證，在現(xiàn)有的標記分水嶺算法基礎上，增加了擴展極大值變換的影像數(shù)據(jù)背景標記，再經(jīng)過強制最小值標定和分水嶺變換后提取海島水邊線信息。實驗中，人工岸線為主的海島水邊線提取精度在1個像素之內(nèi)高于90%，沙質(zhì)岸線為主的海島水邊線提取精度在1.5個像素之內(nèi)高于90%。實驗結(jié)果表明：本文所用的改進算法能夠較好抑制過分割現(xiàn)象，對于人工和沙質(zhì)岸線提取的精度均較高，對于面積較小的島嶼也能準確提取出其水邊線信息。部分區(qū)域受自然因素和影像光譜信息影響，使得分水嶺變換后定位不準確，存在水邊線提取結(jié)果出現(xiàn)細微偏差的現(xiàn)象。同時，本研究缺少海島的多期觀測數(shù)據(jù)和潮汐數(shù)據(jù)，未對其提取的結(jié)果進行潮汐改正計算，僅為海島水邊線。