基于不同岸線類型的象山港岸線提取及變化研究

徐一凱，胡 松，李陽東，朱宇航

（上海海洋大學海洋科學學院，上海 201306）

海岸帶是沿海地區地形變化最快的地區之一[1-2]。海岸線變化原因主要分為自然因素和人為因素[3]，獲得岸線信息對于研究自然與人類環境尤為重要。提取海岸線可通過實地調查獲得[4]，該方法精度較高，但進行大規模海岸線測量非常昂貴和耗時。此外，受到地理條件限制，有時很難對無法進入的地區進行測量。遙感（Remote Sensing）技術則可提供實地調查無法獲得的即時、大規模的海岸線信息[5-7]。

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

1 研究區域及數據方法

1.1 研究區域概況

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

1.2 數據獲取

1.2.1 岸線提取數據

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

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

1.2.2 岸線分類數據

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

1.3 數據方法

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

圖1 岸線提取技術流程

1.3.1 圖像二值化

圖像二值化就是利用閾值將原始圖像分割成前景和背景兩幅圖像。因此，圖像二值化的關鍵是選擇最優閾值，在取閾值時，背景與前景的差異應最大。

（1）大津法（OTSU 法）

最優閾值是衡量差值的標準，而在OTSU 算法中，衡量差值的標準是最大類間方差。

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

類間方差σ2表達式如下所示。

把上式化簡，可得

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

（2）改進的歸一化差異水體指數法（MNDWI 法）

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

針對此問題，本文采用MNDWI 來增大潮濕灘面和干燥灘面的差異，并進行岸線提取。其公式如下。

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

利用此公式將水體凸顯出來，再設置閾值來提取水體，進行水陸分離，實現圖像二值化。

1.3.2 Canny 算子提取水邊線

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

2 數據處理與步驟

2.1 岸線分類步驟

本文根據海岸地理環境特點，將海岸類型分為人工海岸和自然海岸兩個一級類別，其中人工海岸分為養殖區、農田、碼頭、建設用地、裸地5 個二級類別，自然海岸分為淤泥質海岸、砂質海岸、基巖海岸3 個二類級別，具體如表2 所示。

表2 海岸分類

參考文獻[11]，根據象山港特點，本文對于淤泥質海岸和砂質海岸，非海水部分在漲潮作用下可以被海水覆蓋，因此上述海岸的非海水部分應與海水一起標記為“可漲潮區域”，另一方面，基巖海岸在漲潮作用下不被海水覆蓋，應標為“不可漲潮區域”。人工海岸是海岸的一種改良形式，以養殖區、農田、碼頭、施工圍堰等為特征，在漲潮時不被海水覆蓋，應標為“不可漲潮區域”，對于浮閥養殖區進行水下潛水或在漲潮時可被海水覆蓋的區域，應標為“可漲潮區域”[17]。

在進行岸線劃分時，Google Earth 高清影像部分區域顏色相近，如農田和圍塘養殖（圖2 框示區域），無法直接進行準確劃分，所以對象山港整個海岸線進行了實地考察，并最終得出了較為準確的岸線分布情況。

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

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

2.2 岸線提取步驟

2.2.1 圖像預處理

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

由于需要對不同時間、不同傳感器獲取的圖像進行比較，對去條帶后的遙感圖像進行了輻射定標，將圖像的亮度灰度值轉換為絕對的輻射亮度。

然后，對輻射定標之后的圖像進行了大氣校正，以消除大氣影響所造成的輻射誤差，將輻射定標值反演為地表真實信息，從而保證地物波譜信息的準確性。

最后，對各時期圖像進行幾何校正，消除及減弱由各種內外因素造成的遙感圖像的幾何畸變。

2.2.2 OTSU 結果

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

然而，在遙感圖像采集過程中，周期性傳感器偏移或電磁干擾引起的噪聲會降低圖像質量，因此在進行邊緣檢測之前，需要對圖像進行去噪處理。因為中紅外波段對于水陸分離情況較好，所以將經過大氣校正后的圖像中的MIR 灰度圖像進行中值濾波以減少噪聲對結果的影響，然后采用OTSU 法對該圖像進行二值化處理。

2.2.3 MNDWI 結果

本文將影像導入ENVI 中，利用公式（5）將水體凸顯出來，再嘗試并最終設置了合適的閾值來提取水體，進行水陸分離，實現圖像二值化。

2.2.4 Canny 算子提取水邊線結果

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

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

圖3 2017 年二值化后Canny 算子運算結果

3 結果分析

3.1 算法比較

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

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

3.2 岸線類型變化分析

岸線類型與經濟發展、航行安全、生態保護有密切聯系，不同的岸線類型往往對應了不同的用海功能，港口的開發往往選擇基巖岸線，而淤泥質岸線地形可以用作天然的鹽場。象山港地區在進行經濟開發的同時，也秉著生態保護的理念，保證經濟發展的同時，打造多元化的海岸區域[18]。

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

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

表3 2002 年各類岸線統計

表5 2017 年各類岸線占比

研究結果顯示，2002—2017 年，各類自然岸線占比均呈現持續下降趨勢，總占比由2002 年的45.83%下降至2010 年的38.30%，在2017 年僅剩31.94%；而各類人工岸線占比卻在持續上升，總占比由2002 年的54.17%上升至2010 年的61.70%，在2017 年達到了68.06%。其中，自然岸線中的基巖岸線變化較小，而淤泥砂質岸線由2002 年至2017 同比下降44.79%；而人工岸線中的養殖區由2002 年至2017 年同比上升43.05%。

3.3 岸線變遷變化分析

表4 2010 年各類岸線占比

為具體分析其變化原因，本文按照行政區劃圖，把象山港分為郭巨—大嵩、大嵩—桐照、鐵港、黃墩港、峽山—烏沙、西滬港、西澤—錢倉七個區域，然后根據上文所做出的2002 年、2010 年、2017 年岸線變遷圖（圖6），分析15 年間象山港岸線的變化情況。可以看到，2002—2017 年間，郭巨—大嵩、大嵩—桐照末尾、鐵港、黃墩港、西滬港變化較大，其他地區變化較小。本文將變化較大區域用箭頭標明，并對照高清衛星圖對該區域變化原因進行分析。

圖6 2002—2017 年岸線變遷圖

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

大嵩—桐照末尾：2002—2010 年桐照西側岸線出現較大擴張（區域4），2010—2017 年桐照東側岸線出現較大擴張（區域3），2002—2010 年桐照西側區域進行了圍海造陸，而桐照東側出現岸線變化原因如區域1 相似，即將原有島嶼規劃進陸地（圖7（b））。

圖7 2002—2017 年岸線分析

4 結 論

本文以象山港為研究區域，結合實地考察進行遙感圖像岸線自動提取。考慮到不同方法對于不同岸線類型的適用度不同，故本文基于象山港各岸線類型，使用不同方法進行提取，研究結論如下。

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

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

（3）2002—2017 年，象山港區域各類自然岸線占比均呈現持續下降趨勢；反之，各類人工岸線占比均呈現持續上升趨勢。其中，對比2017 年與2002 年數據，自然岸線中的淤泥砂質岸線同比下降，人工岸線中的養殖區同比上升。

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

本文仍存在很多不足之處，由于未進行潮位校正，導致部分砂質及淤泥質提取岸線與真實岸線仍存在誤差，例如泥沙淤積較為嚴重的西滬港區域由于受季節及潮位的影響，2010 年的岸線與2002 年和2017 年的岸線有所差異（圖6 區域5），后續將開展海洋數值模型潮位模擬研究用于校正。