岸外沙洲潮間帶地形的增強型遙感構建方法

周永，張東,2*，鄧慧麗，徐南，張慧銘，郝昕，沈永明

( 1. 南京師范大學 海洋科學與工程學院，江蘇 南京 210023；2. 江蘇省地理信息資源開發與利用協同創新中心，江蘇 南京210023)

1 引言

潮灘位于海陸交互區域，具有獨特的地形[1]、沉積物[2]和水文特征[3]。作為眾多生物的棲息場所，潮灘在保持生物多樣性、維持生態系統穩定等方面具有極高的生態價值[4]。隨著沿海經濟的不斷發展，潮間帶的經濟價值日益凸顯[5]。但人類開發活動如灘涂圍墾、圍海養殖、港口建設等在帶來大量土地和經濟效應的同時，也給沿海地區尤其是潮灘地區造成了巨大影響[6]。近年來，隨著自然影響和人類開發活動力度的加大，潮灘環境變化加快，沖淤變化顯著，潮灘資源正在面臨極大的危機。針對潮灘的動態變化監測與利用保護刻不容緩[7]。

江蘇岸外輻射沙洲作為江蘇淤泥質潮灘的重要部分，在資源開發利用和海岸演變研究等方面具有重要價值[8]。但受較遠的離岸距離、惡劣的勘測環境和周期性淹沒下灘面持續變動等因素影響，岸外輻射沙洲潮間帶區域長期缺少全面、準確的地形監測數據。基于LiDAR[9]、InSAR[10]技術的潮灘地形獲取方法因其在數據獲取性和經濟性上的不足，制約著在該區域的大面積監測應用。遙感水邊線法利用多時相衛星遙感影像數據構建潮灘數字高程模型（DEM），具有數據獲取方便、簡單實用、精度較高等優勢，成為潮間帶地形獲取的重要手段[8,11]。目前該方法已運用于德國瓦登海[12]、韓國戈姆索和咸平郡灣[13]等多個區域的潮灘地形反演，國內Wang等[14]、Kang等[15]和鄭宗生[16]也利用遙感水邊線方法完成了長江口、江蘇沿海等淤泥質潮灘區域的地形構建，并取得較好結果。

輻射沙洲區域水動力條件復雜，日益頻繁的人類活動加劇了水動力環境的調整。雖然在較短周期內，相似潮位下的沙體輪廓范圍大體一致，但沙體表面形態變化（如潮溝擺動）會引起水邊線的位置發生變化，從而導致在較窄的平面條帶范圍內集中的大量水邊線存在相互交叉現象[17]，造成生成的DEM產生地形破碎化和異常起伏。因此，選擇合適的水邊線是構建高精度潮間帶DEM的基礎。Sagar等[18]提出了以長時間序列影像為基礎的潮灘高程信息提取方法，通過將序列影像的時間域轉化為潮高域來構建合成圖層，在獲取平均地形變化特征的同時減少影像異常值的影響。該方法在一定程度上可以借鑒用于解決交叉水邊線的問題，但其應用于澳大利亞大陸尺度的潮灘地形構建中，只關注了單一潮位站點的潮位序時變化以及長周期下大量像元灰度序列數據的統計特征，而忽略了區域內的潮差空間分布差異以及長時間周期內可能的潮灘沖淤變化，因而難以直接應用到輻射沙洲等潮灘地形變化較快的區域。此外，潮間帶沙體的坡度相對平緩，高差不大，灘面分布有大量的樹枝狀和貫通狀潮溝，然而現有的DEM構建方法只關注沙體整體的形態，卻忽略了潮溝的微地貌表達[19]，這也直接影響了潮間帶DEM的構建精度和呈現效果。

因此，本文擬針對常規水邊線方法在應用于地形多變的岸外沙洲潮間帶DEM遙感構建中出現的水邊線交叉誤差和缺乏考慮潮溝微地形的問題，研究增強型潮灘DEM構建方法，實現岸外沙體DEM的快速構建，為利用衛星遙感數據監控潮間帶地形變化提供技術支撐。

2 研究區與數據

2.1 研究區概況

輻射沙洲位于江蘇中部沿海，北起射陽河口，南至蒿枝港，由70多條沙脊和潮流通道組成，沙體表面平坦開闊，在潮汐作用下被周期性淹沒。中部條子泥-蔣家沙-竹根沙區域作為輻射沙洲的核心區，紫菜養殖[20]、風電開發[21]等人類開發活動頻繁，沙體短周期內變化監測需求較大，因而被選作研究區域，進行潮間帶DEM構建試驗。研究區范圍如圖1所示，位于32°38′～32°59′ N，120°58′～121°25′ E之間，南北跨度約37 km，東西跨度超過40 km。潮汐類型以正規半日潮為主，淺海分潮顯著。潮波系統受東海前進潮波與南黃海旋轉潮波共同影響，二者在弶港岸外輻合，波能集中，區域內平均潮差為3.9 m，最大潮差可達9.28 m，且潮差以弶港為中心分別向南、北逐漸減小[22]。

2.2 遙感影像數據

綜合可覆蓋研究區范圍、清晰少云等要求，選擇了25景10 m空間分辨率的Sentinel 2 MSI影像和8景15 m分辨率的Landsat 8 OLI影像作為遙感數據源，進行潮間帶DEM構建，影像的成像時間分布如圖2所示，時間跨度為2017年10月至2019年6月。所有影像經過大氣校正、幾何精校正、圖像增強等預處理，幾何精校正精度優于1個像元。由于空間分辨率不同，影像均重采樣至30 m。

圖2 影像成像時間分布Fig. 2 Image acquisition time

2.3 潮位數據

由于地球曲率變化、淺海潮波變形等影響，輻射沙洲岸段潮位變化呈現出的“跳繩效應”[22]，條子泥-蔣家沙-竹根沙海域的潮差變化較大，無法用單一站點的潮位代表水邊線上不同位置的潮位空間差異。因此，在研究區北、西、南3個方向上分別選擇了大豐港、弶港和洋口港3個潮位站點，進行潮位數據收集，站點位置如圖1所示。根據潮汐調和分析原理，利用站點的短周期（2個月以上）逐時潮位資料，即可計算其潮汐調和常數，實現站點在影像成像時刻的潮位模擬并具有較好精度[23]。因此，收集了上述3個站點2018年6月至2019年1月間逐小時的潮汐表潮位過程數據（http://global-tide.nmdis.org.cn/），利用T_TIDE潮汐調和分析工具進行分潮模擬[24]，根據信噪比大小確定出顯著分潮，模擬得到研究區的潮位誤差小于20 cm。

2.4 DEM驗證數據

收集了2018年蔣家沙-竹根沙局部海域的水下聲吶實測地形數據，空間分辨率為250 m，高程基準為1985國家高程基準。在蔣家沙沙體上按照經緯度方向、在竹根沙沙體上按照平行和垂直于沙脊走向方向劃定4條驗證剖面，分別命名為JJS-H、JJS-Z、ZGS-H和ZGS-Z，位置如圖1所示。通過提取斷面上的實測點高程，驗證遙感方法構建的潮間帶DEM精度，評價DEM質量。

圖1 研究區概況Fig. 1 Sketch map of the study area

3 研究方法

常規遙感水邊線方法是基于潮間帶區域在潮汐作用下被周期性淹沒的特性，將從遙感影像上提取并賦予潮位信息的瞬時水邊線視為“等高線”，利用空間插值處理，得到潮間帶DEM。本文在遙感水邊線法的基礎上，提出面向岸外沙洲潮間帶地形的增強型遙感構建方法，實現以下兩方面的改進：（1）根據較短周期、較小潮差內淹沒概率相似的原理，引入影像序列合成方法，提取合成影像的海陸邊界線，解決常規水邊線方法中對于較窄區域內多時相水邊線交叉重疊造成的DEM高程異常和破碎化難題；（2）提出潮溝DEM構建方法，通過在常規水邊線方法得到的DEM上鑲嵌潮溝DEM，實現對潮灘表面微地形的表達，提高潮間帶DEM精度，改善可視化效果。具體技術流程如圖3所示，主要技術方法闡述如下。

圖3 基于遙感水邊線的增強型地形構建方法技術流程圖Fig. 3 Flowchart of enhanced terrain construction method based on waterlines

3.1 基于空間及潮位關系特征的影像排序

構建高精度的潮間帶DEM，關鍵是用于空間插值的瞬時水邊線在平面上具有位置分離性，在高度變化上具有良好趨勢性。因此，首先需要對從33景預處理后的遙感影像上提取得到瞬時水邊線按照先平面位置、后潮位關系進行排序。

由于潮位滯后效應影響，大豐港、弶港和洋口港3個潮位觀測站點的潮位序時變化不同，需要選擇合適的站點作為影像序列的潮高排序基準。基于此，分別以其中1個站點的潮高為基礎，對全部33景影像按高低排序，然后對3個站點的潮位變化情況進行線性擬合，通過對比趨勢線的擬合決定系數R2和斜率，選擇潮位排序的基準站點。圖4a至圖4c分別顯示了以大豐港、弶港和洋口港潮位站為基準進行影像成像時刻潮位排序后的潮位分布和線性擬合結果，表1列出了對應的擬合決定系數R2和斜率值。

表1 不同排序基準下各站點的潮位變化擬合決定系數與斜率Table 1 The coefficient of determination and slope for fitting tide level changes at different tide stations under different sequencing conditions

相較而言，如圖4c所示，以洋口港站為排序基準下，大豐港站和弶港站潮位變化趨勢線性擬合效果相對較好，大豐港站、弶港站、洋口港站的擬合決定系數R2分別為0.74、0.32和0.93，潮位點變化離散起伏較小，表明以洋口港站為排序基準，大豐港與弶港站的潮位變化趨勢更具有規律性。進一步對比趨勢線斜率，以洋口港站為排序基準時，大豐港站、弶港站、洋口港站潮位趨勢斜率分別為-0.11、-0.06和-0.13，3個站點的潮位變化斜率變化差異最小，潮位變化速率最為接近。因此選定洋口港站作為基準站點進行潮位排序，其中，12條具有正確空間與潮位關系特征的水邊線，后續直接參與DEM插值構建；其余21條水邊線根據潮高大小相似、平面位置相近、線段相交的原則，分成5組，各組潮差范圍在0.13～0.4 m之間，平均潮差間隔0.24 m。5組水邊線對應的影像組如表2所示，用于下一步的影像序列合成及海陸邊界線提取。

表2 各影像組中的影像與成像時刻潮位高度Table 2 Images and their correspondent tide levels for each image group

圖4 各潮位站點在不同排序條件下影像成像時刻的潮位分布和線性擬合Fig. 4 Tide level distribution and linear fitting of imaging time at various tide stations under different sequencing conditions

3.2 影像序列合成與海陸邊界線提取

影像序列合成方法通過統計一個特定周期內影像各像元的灰度值域分布，來指示該像元位置被潮水淹沒的頻率狀況，然后結合潮高變動范圍，將時間序列影像的灰度值域轉變為潮高域，在特定潮高范圍內獲得預設淹沒頻率下的海陸邊界線[18]。以其代替分布于較窄區域內的多時相水邊線，可去除水邊線交叉所帶來的高度異常起伏。

將用于序列合成的遙感影像按組裁剪為相同大小，采用改進的歸一化差值水體指數（MNDWI）進行圖像增強處理[25]，來加大影像中的海陸差異，并進行影像序列合成。MNDWI公式如下：

式中，Rg為 綠波段反射率，在Landsat 8 OLI和Sentinel 2 MSI中同為B3波段；Rswir為短波紅外波段反射率，在Landsat 8 OLI和Sentinel 2 MSI中分別對應B6波段和B11波段。在淤泥質潮灘，由于短波紅外波段對渾濁海水的光譜反射率較低，因此使用MNDWI處理后，海陸邊界色調對比度更強，有利于海陸邊界線的識別。

與以長時間序列影像為基礎，取MNDWI中位數作為合成影像像元值的方式不同[18]，由于各影像組內原始影像的時間跨度短、數量少，經過篩選和排序后各像元MNDWI序列中基本不存在離群值，利用中位數去除異常值影響的作用并不顯著，因此選取序列的平均值來完成新圖像的合成。這樣通過計算像元的MNDWI平均值作為合成圖像的像元值，將較短周期內各影像中呈現的沙體變化進行均值化處理，可以在弱化單一成像時刻下瞬時沙體形態的同時，獲得代表地形平均變化特征的海陸邊界位置信息。

圖5 顯示了以影像組HC-4為例的5景遙感影像的MNDWI序列合成與海陸邊界線提取過程。可以看到，取序列MNDWI均值形成的合成圖像（圖5c）中出露沙體與海水的邊界總體清晰，僅在部分區域如蔣家沙東部和西南部呈現海陸灰度值相似、邊界不明顯狀況。通過觀察該合成圖像的MNDWI分布直方圖，取0.42的谷值作為分割閾值，生成海陸二值圖像（圖5d），進一步利用Sobel邊緣檢測算子提取出海陸邊界線。以邊界線為基礎、300 m為距離構建緩沖區，疊加合成前5景原始影像提取的瞬時水邊線，除部分區域如蔣家沙西北側與陸岸相連區域以及東南側潮溝擺動頻繁區域以外，各瞬時水邊線均位于該緩沖區內，表明序列合成影像提取的海陸邊界線結果能夠較準確的還原該潮位區間內沙體的平均輪廓信息。圖5e顯示了5組合成影像提取的海陸邊界線結果。

圖5 影像組HC-4的MNDWI序列合成與海陸邊界線提取過程Fig. 5 MNDWI sequence synthesis and sea-land boundary extraction of image group HC-4

3.3 初始潮間帶DEM構建

由于提取的水邊線或合成的海陸邊界線不能被簡單視為等高線，因此將提取出的12條具有正確空間與潮位關系特征的水邊線和從合成影像組中提取的5條合成海陸邊界線，按照30 m的間隔進行空間離散，針對其環狀的形態特征，選用反距離平方加權方法進行空間潮位內插[26]。其中，12條水邊線的潮位插值依據為各潮位控制站點在影像成像時刻的潮位，5條合成海陸邊界線的潮位插值依據為各潮位控制站點在合成組內各影像成像時刻潮位的平均值。對所有的離散點通過潮位控制站點的潮位值進行反距離平方加權潮位插值及賦值處理后，利用不規則三角網（TIN）方法，構建得到研究區的初始潮間帶DEM。

3.4 具有潮溝微地形特征的潮間帶DEM構建

潮間帶表面潮溝大多寬度窄、深度淺、高差小，通過空間插值方法直接反演難度較大。因此參考平原河道地形嵌入思路，通過單獨構建潮溝區域地形DEM，然后將其鑲嵌至初始潮間帶DEM上，可以實現對潮間帶微地形的細致模擬[27]。

選取覆蓋研究區的2018年2月23日Sentinel 2 MSI衛星遙感影像進行潮溝邊界線和中線提取，該影像潮灘出露基本為最大范圍，潮溝形態清晰完整。首先參照水邊線提取流程，提取得到矢量格式的潮溝兩側邊界線，用于控制潮溝范圍，潮溝寬度從0 m至575.92 m；然后將兩側邊界線封閉形成面圖層，提取出潮溝中線，用于控制潮溝深度。對潮溝邊界線與中線按照30 m間隔離散，通過前述3個潮位站點的潮位值利用反距離平方加權進行潮位空間插值，得到離散點的高程。對于提取的潮溝邊界線與中線，3個潮位站點的潮位取值方法如下：潮溝邊界線直接取所提取的影像的成像時刻潮位；潮溝中線由于并非直接提取自影像，因此假定所提潮溝的底部在全年最低潮位下均被海水淹沒，選擇排序基準洋口港站在2018年最低潮位時刻對應的3個潮位站點的潮位，作為其潮位插值基準。然后同樣利用不規則三角網（TIN）方法，插值構建出潮溝DEM。利用ArcGIS的圖像鑲嵌工具，將構建后的潮溝DEM以像元為基礎嵌入原始潮間帶DEM，完成對潮溝位置處的地形替換和邊界平滑處理，即可得到最終具有潮溝微地形特征的潮間帶DEM。

4 結果

4.1 DEM構建結果

對比常規遙感水邊線方法與增強型地形構建方法構建的沙體地形，常規遙感水邊線方法以全部33景影像所提取的瞬時水邊線為基礎進行空間插值，構建的潮間帶地形如圖6a所示，高度范圍為-2.30～1.69 m；增強型潮灘構建方法構建的DEM如圖6b所示，高程范圍為-2.68～1.69 m。將二者所構建的DEM作差值運算，得到的差值分布呈現明顯正態分布特征，平均差值為-2.93 cm，差值范圍在±0.3 m以內的區域占64.64%，可見兩者構建的DEM結果整體上保持了相似的地形信息。但是由于增強型地形構建方法進行了潮溝鑲嵌處理，將原始潮灘上的部分平緩地形替換為V字型潮溝地形，造成灘面高程的最低值下探，因此高程變化范圍略大于前者。

圖6 地形構建結果對比Fig. 6 Comparison of terrain construction results

在地形呈現方面，常規遙感水邊線方法構建的DEM表面多呈點狀或線狀分布，形態較為破碎，起伏變化雜亂，且無明顯的潮溝形態。高程較高的沙脊區域分布于蔣家沙西部和高泥的西北部，雖然能看出呈塊狀趨勢，但總體較為破碎。而增強型潮灘構建方法減少了因為水邊線交叉冗余所造成的高度異常起伏與地形破碎，增強了潮間帶灘面的潮溝微地形特征，因此可以看到構建的潮間帶DEM表面起伏過渡平滑自然，地貌單元多呈塊狀分布，潮溝形態清晰，多以線狀或樹枝狀分布于沙體邊緣，尤其在蔣家沙南部較為集中。沙脊區域位置與前者一致，但空間形態更加完整。

4.2 與DEM實測地形對比分析

4.2.1地形精度分析

增強型地形構建方法獲得的潮間帶DEM空間分辨率為30 m，能夠較為精細的描述地形起伏細節，但實測剖面是由250 m間隔下的剖面點組成，不同分辨率基準無法直接進行精度對比，因此首先采用升尺度方法，將模擬DEM的空間分辨率插值為250 m；然后采用平均絕對誤差（Mean Absolute Error, MAE）、均方根誤差（Root Mean Square Error, RMSE）和相關系數（r）3個指標進行DEM高程精度驗證。

根據4條剖面的同名點高程數據對比，利用增強型地形構建方法得到的條子泥-蔣家沙-竹根沙DEM平均的MAE為0.43 m，RMSE為0.54 m，地形起伏變化的相關系數r為0.75。從各斷面來看，如圖7所示，誤差最大剖面為JJS-Z，其MAE與RMSE分別為0.53 m和0.65 m；誤差最小的為ZGS-Z，MAE為0.37 m，RMSE為0.34 m。地形趨勢變化中相關系數最高的斷面同樣為ZGS-Z，r為0.85，最低為JJS-H，r為0.65。4條剖面中，竹根沙沙體的誤差結果較蔣家沙沙體偏小，主要是由于蔣家沙、條子泥區域與陸岸的距離較近，受陸岸圍墾等人類開發活動影響較大，水動力環境復雜，沙體東側的東大港大型潮溝系統向東南方向擺動，潮溝拓寬浚深；西側的西大港潮溝向南貫通，造成蔣家沙沙體沖淤頻繁，灘面小型潮溝變化活躍，對水邊線的提取與DEM插值精度產生影響。而竹根沙區域距離陸岸較遠，水動力環境相對蔣家沙而言較為穩定，沙體在研究周期內變化較小，因此DEM精度相對較高。此外，增強型地形構建方法得到的地形DEM反映的是模擬時段內的平均地形信息，具有一定的地形概化與平均效果，而實測數據顯示的是地表測量時點的瞬時形態，因而DEM獲取方法的差異及數據所代表的時效性也是造成兩者地形之間存在差異的一個重要原因。

圖7 各剖面實測地形與模擬地形的精度和趨勢相關性對比Fig. 7 Comparison of the accuracy and trend correlation between the measured and simulated terrain at each section

4.2.2剖面形態分析

圖8 顯示了增強型地形構建方法獲得的沙體DEM與實測DEM的剖面形態對比。可以看出，JJSH、JJS-Z和ZGS-Z 3條剖面中實測與模擬地形的起伏位置總體基本吻合，但在地形起伏程度上，除了ZGS-H剖面的模擬地形起伏較實測地形稍大以外，其余剖面實測地形的起伏更為顯著。如圖8a所示，JJSH剖面中實測地形在距離起點3 000 m和9 000 m處有各有一個明顯的沙脊隆起，高程分別達到0.3 m和0.8 m；而在模擬地形中，相同位置上同樣有小的地形上升，但上升高度并不明顯。圖8b所示的JJS-Z剖面中，實測地形在距起點2 000 m和7 500 m處都存在明顯的潮溝下切，深度可達到約-2.5 m；在模擬剖面中兩處位置同樣存在潮溝，但下切深度僅達到-2 m左右，較實測地形存在一定差異。這主要是由于序列合成處理將沙體在不同時相的變化平均化，且潮位賦值也同樣以潮高變化的平均值為主，因此相比于實測地形，模擬DEM的高程變化幅度呈現被壓縮狀態。圖8d中ZGS-Z剖面的反演結果在地形起伏方面基本與實測剖面一致，其中在距起點約0～2 000 m的范圍內呈現出地形上升，而后2 000～6 500 m呈現出緩慢下降趨勢，表明竹根沙中沙體北側較南側更陡。

ZGS-H剖面（圖8c）與前面3個剖面呈現出相反特征，在距離起點2 000～8 000 m的范圍中，實測剖面相對較為平緩，但是模擬地形起伏較為明顯。對比影像數據可知，起伏位置都表現為潮溝分布。由于模擬DEM的空間分辨率高，結果能呈現更為細致的起伏結果，且升尺度后的同名位置點仍然能保持周期影像內的起伏信息，但實測地形在250 m的分辨率下，部分較小潮溝信息無法呈現，因而該段區域內的地形呈現出輕微起伏下的平緩狀況。

圖8 剖面形態對比Fig. 8 Morphology comparison of each section

5 結論

本文提出了一種岸外沙洲潮間帶增強型地形遙感構建方法，通過將常規遙感水邊線法與序列影像合成法相結合，完成小范圍和短周期下的潮灘地形反演；進一步引入潮溝DEM鑲嵌技術，實現了表現潮灘微地貌特征的岸外沙洲潮間帶DEM快速構建。研究結果表明：

（1）根據模擬的潮間帶DEM與實測地形之間的4條驗證剖面對比，平均的MAE為0.43 m，RMSE為0.54 m，相關系數為 0.75，兩者的地形起伏位置較為吻合。但是由于增強型地形構建方法得到的是模擬周期內的平均地形變化狀況，因此地形整體起伏偏小。

（2）對比增強型地形遙感構建方法與常規遙感水邊線法得到的潮灘DEM結果，兩者的平均差值為-2.93 cm，差值呈現明顯的正態分布特征且范圍在±0.3 m以內的區域占64.64%，可見增強型地形遙感構建方法能夠較好地保留潮灘平均高度變化信息。同時該方法減少了因灘面潮溝擺動導致水邊線交叉造成的高程誤差和地形破碎，構建的潮間帶DEM表面起伏過渡自然，潮溝形態清晰，因此適用于淤泥質潮間帶地形的快速構建，并保證較高的DEM模擬精度。

該方法在合成影像組的確定過程中，需要首先對所有影像的水邊線進行提取，而后在水邊線空間排序的基礎上再對比潮位高度的變化范圍，從而確定影像組對應的影像數據。但是從以上影像組提取的水邊線并未參與后續的潮間帶DEM構建，所以后續研究需要對影像中沙體出露范圍判斷和影像成像時刻潮位匹配流程進行改進優化，來提高本潮間帶增強型地形遙感構建方法的建模效率。

致謝：感謝歐洲航天局（ESA）和美國地質調查局（USGS）提供的免費衛星影像資料，以及國家海洋信息中心提供的潮位資料信息。