徐一凱，胡 松，李陽(yáng)東，朱宇航

（上海海洋大學(xué)海洋科學(xué)學(xué)院，上海 201306）

海岸帶是沿海地區(qū)地形變化最快的地區(qū)之一[1-2]。海岸線變化原因主要分為自然因素和人為因素[3]，獲得岸線信息對(duì)于研究自然與人類環(huán)境尤為重要。提取海岸線可通過(guò)實(shí)地調(diào)查獲得[4]，該方法精度較高，但進(jìn)行大規(guī)模海岸線測(cè)量非常昂貴和耗時(shí)。此外，受到地理?xiàng)l件限制，有時(shí)很難對(duì)無(wú)法進(jìn)入的地區(qū)進(jìn)行測(cè)量。遙感（Remote Sensing）技術(shù)則可提供實(shí)地調(diào)查無(wú)法獲得的即時(shí)、大規(guī)模的海岸線信息[5-7]。

利用遙感技術(shù)提取海岸線的方法主要有目視解譯和自動(dòng)提取。目視解譯依據(jù)專業(yè)人員的主觀判斷，通過(guò)目視觀察提取岸線，操作簡(jiǎn)單、提取岸線精度高，但工作效率低、工作量大，不適合快速提取大范圍遙感影像數(shù)據(jù)。與目視解譯相比，自動(dòng)提取具有提取速度快、效率高的優(yōu)勢(shì)，但存在同物異譜和同譜異物的問(wèn)題，因此，需要設(shè)計(jì)科學(xué)的提取算法。常用的自動(dòng)提取方法有：基于邊緣檢測(cè)的方法[8]、基于指數(shù)分析的方法[9]、基于閾值分割的方法[10]、基于區(qū)域生長(zhǎng)的方法[11]、基于神經(jīng)網(wǎng)絡(luò)的方法[12]和亞像素方法[13]。然而，在進(jìn)行自動(dòng)提取時(shí)，往往利用單一算法對(duì)岸線進(jìn)行提取[14]，未考慮單一算法對(duì)不同類型岸線提取的適用性，所以本文針對(duì)不同海岸類型評(píng)測(cè)不同方法對(duì)其的適用性，采取適用方法進(jìn)行該海岸類型的岸線提取，從而提高岸線提取精度，并對(duì)提取出的岸線進(jìn)行分析，總結(jié)研究區(qū)近年來(lái)岸線變化的原因。

1 研究區(qū)域及數(shù)據(jù)方法

1.1 研究區(qū)域概況

研究區(qū)域象山港位于寧波市東南部（29?24′N—29?46′ N，121?25′ E—122?00′ E），地處穿山半島和象山半島之間，沿東北—西南走向?yàn)橐粋€(gè)狹長(zhǎng)的半封閉式海灣。全港長(zhǎng)達(dá)60 多千米，水深一般10~15 m，海岸線蜿蜒曲折。港內(nèi)有西滬港、黃墩港、鐵港3 個(gè)分港。

1.2 數(shù)據(jù)獲取

1.2.1 岸線提取數(shù)據(jù)

鑒于象山港區(qū)域近15 年間岸線變化較為劇烈，研究其岸線變遷需要進(jìn)行較長(zhǎng)時(shí)間跨度的對(duì)比分析，故本文從地理空間數(shù)據(jù)云（http：//www.gscloud.cn）選取了2002 年1 月3 日的Landsat 5 TM、2010年11 月1 日和2017 年3 月9 日Landsat 7 ETM 影像數(shù)據(jù)。其中，Landsat TM 影像包含7 個(gè)波段，Landsat ETM+影像數(shù)據(jù)包括8 個(gè)波段（表1）。

表1 Landsat TM 波段和Landsat ETM+波段

1.2.2 岸線分類數(shù)據(jù)

為比較不同的岸線提取方法在提取不同岸線時(shí)的精確程度，本文采用2017 年9 月的谷歌地球（Google Earth）全球高清影像進(jìn)行岸線分類。通過(guò)目視解譯方法，以離岸線500 m 以內(nèi)區(qū)域的主要用地類型作為該段岸線的類型進(jìn)行手動(dòng)劃分。

1.3 數(shù)據(jù)方法

不同類型的海岸有不同的地面特征，單一的算法通常難以保證岸線提取的準(zhǔn)確性，故本文采用了不同的閾值分割法對(duì)象山港區(qū)域影像進(jìn)行了圖像二值化處理，再對(duì)二值化后的圖像利用Canny 邊緣檢測(cè)算子進(jìn)行了邊緣線提取。具體技術(shù)流程如圖1所示。

圖1 岸線提取技術(shù)流程

1.3.1 圖像二值化

圖像二值化就是利用閾值將原始圖像分割成前景和背景兩幅圖像。因此，圖像二值化的關(guān)鍵是選擇最優(yōu)閾值，在取閾值時(shí)，背景與前景的差異應(yīng)最大。

（1）大津法（OTSU 法）

最優(yōu)閾值是衡量差值的標(biāo)準(zhǔn)，而在OTSU 算法中，衡量差值的標(biāo)準(zhǔn)是最大類間方差。

OTSU 算法的原理是將圖像中每個(gè)像素點(diǎn)的灰度值與假定的閾值T進(jìn)行比較，將像素灰度值小于閾值T的歸類為C1，大于閾值T的歸類為C2。假設(shè)歸類為C1和C2的均值分別為m1、m2，圖像全局均值為mG，而此時(shí)像素被分為C1和C2類的概率分別為p1、p2，計(jì)算公式如下所示。

類間方差σ2表達(dá)式如下所示。

把上式化簡(jiǎn)，可得

然后，遍歷所有灰度級(jí)，求出使σ2最大的灰度值k就是OTSU 閾值。

（2）改進(jìn)的歸一化差異水體指數(shù)法（MNDWI 法）

遙感圖像對(duì)于不同類型的地物反射率不同，從而呈現(xiàn)出的亮度也不同。如砂質(zhì)海岸的干燥灘面反射率較高，在遙感圖像上表現(xiàn)為較亮區(qū)域；與之相反，潮濕灘面反射率較低，在遙感圖像上表現(xiàn)為較暗區(qū)域；由于海水反射率最低，在遙感圖像上表現(xiàn)為最暗區(qū)域。遙感圖像瞬時(shí)水邊線應(yīng)為潮濕灘面和干燥灘面的分界線，由于此分界線并不明顯，若直接提取，則精度不高[15]。

針對(duì)此問(wèn)題，本文采用MNDWI 來(lái)增大潮濕灘面和干燥灘面的差異，并進(jìn)行岸線提取。其公式如下。

式中，P（Green）和P（MIR）分別代表綠波段（Landsat TM/ETM 中的波段2） 和中紅外波段（Landsat TM/ETM 中的波段5）的亮度值。

利用此公式將水體凸顯出來(lái)，再設(shè)置閾值來(lái)提取水體，進(jìn)行水陸分離，實(shí)現(xiàn)圖像二值化。

1.3.2 Canny 算子提取水邊線

Canny 邊緣檢測(cè)算子是一種多級(jí)邊緣檢測(cè)算法。盡管Canny 算子的操作較為復(fù)雜，但相比Sobel 邊緣檢測(cè)算子，其能處理包含噪聲較多的圖像，且具有Robert 算子定位精確的優(yōu)點(diǎn)。Canny 算子具有提取位置精度高，包含信息多的優(yōu)點(diǎn)，近年來(lái)已被廣泛用于邊緣提取等圖形檢測(cè)中[16]。Canny 算子原理包括以下4 個(gè)步驟：①用高斯濾波器平滑處理原圖像；②梯度幅值和方向的計(jì)算；③消除冗余窗口；④雙閾值算法檢測(cè)和連接邊緣。

2 數(shù)據(jù)處理與步驟

2.1 岸線分類步驟

本文根據(jù)海岸地理環(huán)境特點(diǎn)，將海岸類型分為人工海岸和自然海岸兩個(gè)一級(jí)類別，其中人工海岸分為養(yǎng)殖區(qū)、農(nóng)田、碼頭、建設(shè)用地、裸地5 個(gè)二級(jí)類別，自然海岸分為淤泥質(zhì)海岸、砂質(zhì)海岸、基巖海岸3 個(gè)二類級(jí)別，具體如表2 所示。

表2 海岸分類

參考文獻(xiàn)[11]，根據(jù)象山港特點(diǎn)，本文對(duì)于淤泥質(zhì)海岸和砂質(zhì)海岸，非海水部分在漲潮作用下可以被海水覆蓋，因此上述海岸的非海水部分應(yīng)與海水一起標(biāo)記為“可漲潮區(qū)域”，另一方面，基巖海岸在漲潮作用下不被海水覆蓋，應(yīng)標(biāo)為“不可漲潮區(qū)域”。人工海岸是海岸的一種改良形式，以養(yǎng)殖區(qū)、農(nóng)田、碼頭、施工圍堰等為特征，在漲潮時(shí)不被海水覆蓋，應(yīng)標(biāo)為“不可漲潮區(qū)域”，對(duì)于浮閥養(yǎng)殖區(qū)進(jìn)行水下潛水或在漲潮時(shí)可被海水覆蓋的區(qū)域，應(yīng)標(biāo)為“可漲潮區(qū)域”[17]。

在進(jìn)行岸線劃分時(shí)，Google Earth 高清影像部分區(qū)域顏色相近，如農(nóng)田和圍塘養(yǎng)殖（圖2 框示區(qū)域），無(wú)法直接進(jìn)行準(zhǔn)確劃分，所以對(duì)象山港整個(gè)海岸線進(jìn)行了實(shí)地考察，并最終得出了較為準(zhǔn)確的岸線分布情況。

本文利用了Google Earth 全球真彩色影像，結(jié)合上述實(shí)地考察，對(duì)岸線類型進(jìn)行目視解譯并提取岸線，將其導(dǎo)入ArcGIS 進(jìn)行要素的整合，獲得岸線類型及分布如圖2 所示。

圖2 象山港岸線類型及分布

2.2 岸線提取步驟

2.2.1 圖像預(yù)處理

為獲得象山港區(qū)域的完整圖像，將同時(shí)期的兩幅Landsat 影像并進(jìn)行拼接處理，并按照象山港的經(jīng)緯度范圍進(jìn)行裁剪。由于2003 年Landsat7 ETM+機(jī)載掃描行校正器故障，導(dǎo)致其之后圖像會(huì)出現(xiàn)條帶，所以本文使用了條帶插件將2010 年及2017 年的影像進(jìn)行修復(fù)。

由于需要對(duì)不同時(shí)間、不同傳感器獲取的圖像進(jìn)行比較，對(duì)去條帶后的遙感圖像進(jìn)行了輻射定標(biāo)，將圖像的亮度灰度值轉(zhuǎn)換為絕對(duì)的輻射亮度。

然后，對(duì)輻射定標(biāo)之后的圖像進(jìn)行了大氣校正，以消除大氣影響所造成的輻射誤差，將輻射定標(biāo)值反演為地表真實(shí)信息，從而保證地物波譜信息的準(zhǔn)確性。

最后，對(duì)各時(shí)期圖像進(jìn)行幾何校正，消除及減弱由各種內(nèi)外因素造成的遙感圖像的幾何畸變。

2.2.2 OTSU 結(jié)果

由于基巖海岸和人工海岸在遙感圖像中都具有較明顯的水陸分界，利用邊緣檢測(cè)便可提取水邊線，其因不考慮海岸線的背景差異，因此獲得的海岸線位置一般較為準(zhǔn)確。

然而，在遙感圖像采集過(guò)程中，周期性傳感器偏移或電磁干擾引起的噪聲會(huì)降低圖像質(zhì)量，因此在進(jìn)行邊緣檢測(cè)之前，需要對(duì)圖像進(jìn)行去噪處理。因?yàn)橹屑t外波段對(duì)于水陸分離情況較好，所以將經(jīng)過(guò)大氣校正后的圖像中的MIR 灰度圖像進(jìn)行中值濾波以減少噪聲對(duì)結(jié)果的影響，然后采用OTSU 法對(duì)該圖像進(jìn)行二值化處理。

2.2.3 MNDWI 結(jié)果

本文將影像導(dǎo)入ENVI 中，利用公式（5）將水體凸顯出來(lái)，再嘗試并最終設(shè)置了合適的閾值來(lái)提取水體，進(jìn)行水陸分離，實(shí)現(xiàn)圖像二值化。

2.2.4 Canny 算子提取水邊線結(jié)果

基于在上文中提及的Canny 算子步驟，本文在OTSU 算法所得二值圖像的基礎(chǔ)上對(duì)其進(jìn)行形態(tài)學(xué)操作，然后通過(guò)Matlab 編程利用Canny 算子分別對(duì)兩種二值化方法得到的結(jié)果圖進(jìn)行邊緣檢測(cè)并消除孔洞，以保留主要岸線及島嶼輪廓。以2017 年為例，結(jié)果如圖3 所示。

如圖3 所示，利用Canny 算子提取出的水邊線連續(xù)、清晰且較為平滑，虛假邊緣被很好地抑制，受噪聲干擾小，提取效果較好。將Canny 算子計(jì)算過(guò)后的圖像導(dǎo)入ArcGIS 中并進(jìn)行柵格數(shù)據(jù)轉(zhuǎn)矢量數(shù)據(jù)處理得到矢量海岸線，并將不同年份的岸線進(jìn)行疊合以進(jìn)行后續(xù)的觀察和分析。

圖3 2017 年二值化后Canny 算子運(yùn)算結(jié)果

3 結(jié)果分析

3.1 算法比較

如上所述，本文分別采用了OTSU 法與MNDWI 法對(duì)圖像進(jìn)行了二值化處理，使用MNDWI法提取出的岸線比使用OTSU 法提取的岸線更向海推進(jìn)，尤其是箭頭指向區(qū)域。將真彩色影像進(jìn)行疊加分析，發(fā)現(xiàn)相差部分為淤泥質(zhì)和砂質(zhì)海岸（圖4區(qū)域1 和區(qū)域2），從而驗(yàn)證了MNDWI 法相比OTSU 法更適合提取砂質(zhì)和淤泥質(zhì)海岸。而在提取基巖岸線及人工岸線時(shí)，OTSU 法更加精細(xì)及平滑（圖4 區(qū)域3）。

基于上述情況，本文將砂質(zhì)岸線、淤泥質(zhì)岸線及浮筏養(yǎng)殖區(qū)按照MNDWI 法提取，將人工岸線和基巖岸線按照OTSU 法提取，將兩者進(jìn)行矢量裁剪并拼接，得到最終的象山港岸線。

3.2 岸線類型變化分析

岸線類型與經(jīng)濟(jì)發(fā)展、航行安全、生態(tài)保護(hù)有密切聯(lián)系，不同的岸線類型往往對(duì)應(yīng)了不同的用海功能，港口的開(kāi)發(fā)往往選擇基巖岸線，而淤泥質(zhì)岸線地形可以用作天然的鹽場(chǎng)。象山港地區(qū)在進(jìn)行經(jīng)濟(jì)開(kāi)發(fā)的同時(shí)，也秉著生態(tài)保護(hù)的理念，保證經(jīng)濟(jì)發(fā)展的同時(shí)，打造多元化的海岸區(qū)域[18]。

為研究各岸線類型的占比變化，本文利用ArcGIS 分別對(duì)2002 年、2010 年及2017 年各岸線類型進(jìn)行了長(zhǎng)度計(jì)算，結(jié)果如表3 至表5 及圖5 所示。

圖5 2002—2017 年各岸線占比變化

表3 2002 年各類岸線統(tǒng)計(jì)

表5 2017 年各類岸線占比

研究結(jié)果顯示，2002—2017 年，各類自然岸線占比均呈現(xiàn)持續(xù)下降趨勢(shì)，總占比由2002 年的45.83%下降至2010 年的38.30%，在2017 年僅剩31.94%；而各類人工岸線占比卻在持續(xù)上升，總占比由2002 年的54.17%上升至2010 年的61.70%，在2017 年達(dá)到了68.06%。其中，自然岸線中的基巖岸線變化較小，而淤泥砂質(zhì)岸線由2002 年至2017 同比下降44.79%；而人工岸線中的養(yǎng)殖區(qū)由2002 年至2017 年同比上升43.05%。

3.3 岸線變遷變化分析

表4 2010 年各類岸線占比

為具體分析其變化原因，本文按照行政區(qū)劃圖，把象山港分為郭巨—大嵩、大嵩—桐照、鐵港、黃墩港、峽山—烏沙、西滬港、西澤—錢倉(cāng)七個(gè)區(qū)域，然后根據(jù)上文所做出的2002 年、2010 年、2017 年岸線變遷圖（圖6），分析15 年間象山港岸線的變化情況?？梢钥吹?，2002—2017 年間，郭巨—大嵩、大嵩—桐照末尾、鐵港、黃墩港、西滬港變化較大，其他地區(qū)變化較小。本文將變化較大區(qū)域用箭頭標(biāo)明，并對(duì)照高清衛(wèi)星圖對(duì)該區(qū)域變化原因進(jìn)行分析。

圖6 2002—2017 年岸線變遷圖

主要變化如下：郭巨—大嵩：2002—2010 年大嵩附近岸線出現(xiàn)較大的擴(kuò)張（區(qū)域2），2010—2017年出現(xiàn)巨大的擴(kuò)張變化（區(qū)域1），查閱資料得知，2002—2010 年的岸線變化是由于圍海造陸所致，而2010—2017 年的變化是由于凸出的部分是原有島嶼（梅山島），在2010—2017 年間梅山水庫(kù)的建設(shè)將島嶼和大陸連接，從而致使岸線出現(xiàn)了巨大變化（圖7（a））。

大嵩—桐照末尾：2002—2010 年桐照西側(cè)岸線出現(xiàn)較大擴(kuò)張（區(qū)域4），2010—2017 年桐照東側(cè)岸線出現(xiàn)較大擴(kuò)張（區(qū)域3），2002—2010 年桐照西側(cè)區(qū)域進(jìn)行了圍海造陸，而桐照東側(cè)出現(xiàn)岸線變化原因如區(qū)域1 相似，即將原有島嶼規(guī)劃進(jìn)陸地（圖7（b））。

圖7 2002—2017 年岸線分析

4 結(jié) 論

本文以象山港為研究區(qū)域，結(jié)合實(shí)地考察進(jìn)行遙感圖像岸線自動(dòng)提取?？紤]到不同方法對(duì)于不同岸線類型的適用度不同，故本文基于象山港各岸線類型，使用不同方法進(jìn)行提取，研究結(jié)論如下。

（1）基于Google Earth 全球高清影像并結(jié)合實(shí)地考察，將象山港岸線類型劃分為人工海岸和自然海岸兩個(gè)一級(jí)類，其中自然海岸包含淤泥質(zhì)岸線、砂質(zhì)岸線和基巖岸線類型，人工岸線包含農(nóng)田、碼頭、建設(shè)用地和裸地類型。

（2）采用Landsat 影像，結(jié)合岸線類型提出水邊線綜合提取策略。與真彩色影像疊合比較發(fā)現(xiàn)，按照基巖海岸和人工海岸采取OTSU 方法，砂質(zhì)海岸和淤泥質(zhì)海岸采用MNDWI 方法，再使用改進(jìn)的Canny 邊緣檢測(cè)算子提取水邊線，此方法提取出的水邊線較為貼合真實(shí)水邊線。

（3）2002—2017 年，象山港區(qū)域各類自然岸線占比均呈現(xiàn)持續(xù)下降趨勢(shì)；反之，各類人工岸線占比均呈現(xiàn)持續(xù)上升趨勢(shì)。其中，對(duì)比2017 年與2002 年數(shù)據(jù)，自然岸線中的淤泥砂質(zhì)岸線同比下降，人工岸線中的養(yǎng)殖區(qū)同比上升。

（4）由于人工建設(shè)將自然島嶼與大陸相連并修建水庫(kù)，致使郭巨—大嵩、大嵩—桐照末尾的岸線出現(xiàn)了巨大變化；鐵港、黃墩港變化較大則是由于圍塘養(yǎng)殖。

本文仍存在很多不足之處，由于未進(jìn)行潮位校正，導(dǎo)致部分砂質(zhì)及淤泥質(zhì)提取岸線與真實(shí)岸線仍存在誤差，例如泥沙淤積較為嚴(yán)重的西滬港區(qū)域由于受季節(jié)及潮位的影響，2010 年的岸線與2002 年和2017 年的岸線有所差異（圖6 區(qū)域5），后續(xù)將開(kāi)展海洋數(shù)值模型潮位模擬研究用于校正。