基于多算法水邊線提取的潮灘DEM構建

李天祺，王建超，吳芳，趙政，張文凱

(中國自然資源航空物探遙感中心，北京 100083)

0 引言

淤泥質潮灘通常指平均大潮高潮線與平均大潮低潮線之間的地帶，也稱潮間帶[1-2]。在我國約有四分之一的海岸屬于淤泥質海岸[3]，獲取淤泥質潮灘地形對研究岸線變遷、海岸帶生態變化及沿海工程建設有著重要意義[4-5]。淤泥質潮灘受波浪、泥沙、沿岸流和地質等諸多因素的影響，具有面積寬廣、水淺灘平、變化頻繁的特點，使得該區域地面調查與地形測繪的難度較大、成本較高，且存在一定風險[6-7]。遙感技術具有高時效、大范圍、高頻率的特點，利用多期遙感水邊線復合潮位數據構建潮灘數字高程模型(digital elevation model,DEM)，成為獲取大范圍淤泥質潮灘地形信息的有效途徑。

遙感水邊線是衛星過境時獲取的瞬時水邊線，準確提取水邊線是潮灘DEM構建的關鍵。現有的水邊線提取方法主要有： 邊緣檢測法[8]、閾值分割法[9]、區域生長法[10]、主動輪廓模型法[11]、面向對象分類[12]等。Mason等[13]基于ERS SAR影像，通過紋理分割提取多期水邊線，應用水動力模型對水邊線高程進行賦值，構建了英國東岸Humber/Wash區域潮灘DEM； 沈芳等[14]基于Landsat TM影像，對比了不同波段的水邊線提取精度，利用閾值分割法提取水邊線，通過潮汐預測數據插值對水邊線賦值，構建了長江口九段沙DEM； 穆敬等[15]基于BJ-1影像，利用面向對象分類法提取水邊線，通過潮汐網格數據對水邊線賦值，構建了黃驊市潮灘DEM。目前國內外學者在構建潮灘DEM的研究中大多關注于準確模擬水邊線的瞬時潮位和提升單一方法的水邊線提取精度[8-15]。但是，瞬時水邊線受衛星過境時的潮情、天氣以及灘面等因素的影響，在影像中的光譜與紋理差異較大，基于單一方法提取水邊線，難以保證水邊線的提取精度，從而影響潮灘DEM反演結果。吳迪等[16]雖提出了基于多源多算法的水邊線提取模型，但未對水邊線進行系統劃分，且沒有將其應用于潮灘DEM構建研究。

針對上述問題，本文從邊緣類型與灘面噪聲2方面劃分水邊線，根據瞬時潮情對水邊線進行分類提取，通過潮汐數據構建多項式計算瞬時潮位，提出基于多算法水邊線提取的潮灘DEM構建方法，并通過野外實測數據對潮灘DEM反演精度進行驗證。

1 研究區概況及其數據源

1.1 研究區概況

研究區位于N37°21.5′～37°32.0′，E118°53.5′～118°58.5′，位于萊州灣西側，濰坊港驗潮站以北，如圖1所示。研究區主要為淤泥質潮灘，灘面廣闊，坡度平緩，南北長約20 km，寬度為2～4 km，近岸區域有人工堤壩和大量鹽田，灘面生長高耐鹽紅海灘植被，近海區域無植被生長。潮汐類型屬于不規則半日潮混合潮，平均潮差為125 cm，最大潮差為259 cm。

圖1 研究區地理位置Fig.1 Location map of the study area

1.2 研究數據源及預處理

本文選用GF-1 WFV數據，空間分辨率為16 m，重訪周期為2 d，幅寬為800 km。數據獲取時間為2017年2—7月，篩除受厚云、雪及海冰污染的數據，共選取了18期數據用于研究。數據預處理包括絕對輻射定標、幾何精糾正及研究區裁剪，其中幾何精糾正誤差控制在0.5個像元以內。

由于受野外測量條件限制，缺少研究區的實測潮位數據，本文使用濰坊港驗潮站的潮汐預測數據用于水邊線的分類與衛星過境瞬時的水邊線潮位計算。潮汐預測數據是在潮汐觀測數據的基礎上，利用調和分析法推算的預測潮位[14]。濰坊港驗潮站所在海域為不規則半日混合潮，潮汐預測數據包含每日整時的潮位值，以及每日的潮位極值與對應時刻，潮高基準面為平均海平面以下120 cm。

2 多算法水邊線提取

2.1 基于潮情的水邊線分類

水邊線作為遙感影像中一種重要的邊緣信息，是水陸區域的分界。不同遙感水邊線提取算法具有不同的適用范圍[17]，其提取效果主要取決于目標邊緣的類型與邊緣鄰域內的噪聲水平。根據邊緣類型與噪聲水平對水邊線分類，選擇適用算法進行提取，有利于準確地獲取水邊線信息。

水邊線兩側水陸區域的灰度變化梯度主要受漲落潮的影響： 漲潮時灘面干濕分明，水邊線兩側區域灰度差異顯著； 落潮時受灘面殘余水和水中泥沙影響，水邊線兩側灰度差異較小。因此按灰度變化梯度，邊緣類型可分為強邊緣與弱邊緣。水邊線的鄰域噪聲水平，主要由潮位高低決定： 由于灘面地物種類由近海向內陸逐漸增多，高潮位時水邊線在潮溝、植被、建筑物分布較多的近岸區域，灘面破碎復雜，鄰域內噪聲較高； 低潮位時水邊線近海，灘面相對均質單一。

綜上，本文通過分析潮汐數據，確定各期水邊線的瞬時潮情： 漲、落潮與高、低潮位，根據瞬時潮情將水邊線劃分為單一灘面弱邊緣、復雜灘面弱邊緣、單一灘面強邊緣和復雜灘面強邊緣4類，具體分類規則如圖2所示。

Test and optimization on heat transfer performance of electric heater under lower air pressure

圖2 分類規則

在本實驗中，影像中的薄云或霧明顯削弱灰度變化梯度，因此通過人工目視對影像質量進行判斷，將包含薄云或霧影像的水邊線分為弱邊緣水邊線； 通過分析潮汐數據對晴朗影像的水邊線進行劃分。若H(10),H(11)和H(12)3個時刻的潮位依次升高，則瞬時潮情(11—12時之間)為漲潮，灘面不存在殘余水的影響，則水邊線劃分為強邊緣水邊線，若3個時刻的潮位依次下降(落潮)或10時的潮位高于11時潮位，則存在灘面殘余水或水中泥沙的影響，劃分為弱邊緣水邊線。

獲得強、弱邊緣水邊線后，按潮位高低對水邊線所在灘面特征分類： 若水邊線獲取當日的11時與12時的平均潮位高于多期平均值H，則該期水邊線近岸，對應復雜灘面； 若低于多期平均值H則該期水邊線近海，對應單一灘面。最后得到4類水邊線分別為： 單一灘面弱邊緣水邊線、復雜灘面弱邊緣水邊線、單一灘面強邊緣水邊線與復雜灘面強邊緣水邊線。對本文所用的18期GF-1 WFV影像的水邊線進行分類，分類結果如表1所示。

表1 水邊線分類結果Tab.1 Waterline classification results

2.2 基于邊緣檢測水邊線提取

邊緣檢測是基于圖像灰度的一階或二階導數來提取圖像中灰度急劇變化的區域邊界[18]，具有模型參數少、運行速度快的優點，但對噪聲敏感，在復雜灘面提取的水邊線存在嚴重的毛刺或不連續。本文選擇定位準確度高、單邊響應良好的Canny邊緣檢測算法[19]，以20170429影像為例，對單一灘面強邊緣水邊線進行提取，結果如圖3所示。

圖3 基于Canny邊緣檢測的水邊線提取結果Fig.3 Waterline extraction resultbased on Canny edge detection

2.3 基于閾值分割水邊線提取

由于弱邊緣水邊線兩側灰度變化不明顯，邊緣檢測不能有效提取邊緣信息，閾值分割通過設定閾值，如歸一化水體指數(normalized different water index,NDWI)[20]，可以快速有效地區分水陸區域，但對于高潮位的復雜灘面，受潮溝或其他積水區域的影響，水邊線提取不準確，因此，本文通過閾值分割提取單一灘面弱邊緣水邊線。

以20170617影像為例，首先，計算NDWI，增強影像中的水體信息，計算公式為：

(1)

式中Green與NIR分別代表GF-1 WFV影像中的綠光波段與近紅外波段的反射率。NDWI的取值范圍為[-1, 1]； 然后，通過OTSU算法[21]計算閾值，對NDWI灰度影像進行二值化； 最后，對二值圖做開運算處理，對邊界進行平滑，得到水邊線提取結果，如圖4所示。

圖4 基于NDWI閾值分割的水邊線提取結果Fig.4 Waterline extraction result based onNDWI threshold segmentation

2.4 基于面向對象法水邊線提取

面向對象法是通過影像中的光譜信息、空間信息和紋理信息對影像進行分割，然后以分割后的對象為基本單元進行分類，從而提取水邊線的一種方法[22]。面向對象法具有較強的抗噪能力，有效抑制復雜灘面上其他地物的影響，可準確保留強邊緣水邊線的細節信息，但該方法對紋理變化過于敏感，提取的弱邊緣水邊線不平滑，因此，本文通過面向對象法提取復雜灘面強邊緣水邊線。

首先，通過多尺度分割將影像分割為多個對象，然后根據各個對象的光譜相似度進行合并，多次迭代后最終得到水邊線提取結果，以20170609影像為例，初始分割結果和水邊線提取結果如圖5所示。

2.5 基于改進分水嶺水邊線提取

圖6 改進分水嶺算法流程Fig.6 Flow chart of improved watershed algorithm

以20170310影像為例，對影像NDWI低通濾波去除高頻噪聲后，潮灘上的潮溝、河口處的沙洲及較小的人工建筑物均已去除； 二值化處理后通過數學形態學處理，將鹽田、河流等內陸水體完全去除，在邊界處做緩沖區處理，生成待分割區域的標記信息； 利用原始影像B4(R)，B3(G)，B2(B)假彩色圖像結合標記信息進行分水嶺分割，提取結果如圖7所示。

圖7 基于改進分水嶺的水邊線提取結果Fig.7 Waterline extraction result basedon improved watershed algorithm

3 潮灘DEM構建與精度驗證

3.1 潮灘DEM構建

利用濰坊港驗潮站的潮汐預測數據對影像獲取的瞬時潮位進行插值。本文所用GF-1 WFV數據的獲取時間均在每日上午11—12時之間，利用10—13時的4個整時潮位值構建三次多項式，若在10—13時之間存在潮位極值，則用潮位極值替代距影像獲取時間較遠的整點潮位值，通過多項式插值計算水邊線的瞬時潮位。在研究資料有限時，這是一種行之有效的方法[14]，水邊線對應瞬時潮位如表2所示。

表2 水邊線瞬時潮位Tab.2 Instantaneous tidal level of waterlines (cm)

受灘面上潮溝或河流入海口處泥沙的影響，提取的多期水邊線在局部存在交叉或重疊的現象。針對上述問題，本文按照單一灘面強邊緣、復雜灘面強邊緣、單一灘面弱邊緣和復雜灘面弱邊緣的優先順序選用水邊線，結果如圖8所示。

圖8 水邊線選用結果Fig.8 Selection result of waterlines

將篩選后的水邊線重采樣為16 m×16 m的水邊點，利用不規則三角網進行空間插值，并將插值結果從潮高基準面轉換為以WGS84為基準的高程數據，得到潮灘DEM反演結果，如圖9所示。

圖9 潮灘DEM反演結果Fig.9 Inversion result of tidal flat DEM

3.2 精度驗證

為驗證基于多算法水邊線提取的潮灘DEM構建方法有效性，本文利用野外實測高程數據對潮灘DEM反演結果進行精度評價。實測高程數據獲取時間為2017年12月，采用Trimble R5 GPS接收機測量，靜態觀測不少于20分鐘。實測高程數據通過精密星歷解算為WGS84坐標數據，18個驗證點分布如圖10所示。

圖10 驗證點分布Fig.10 Distribution of verification points

對實測高程與反演結果進行相關分析，實測高程與對應DEM反演結果的相關性較高，R2為0.864 9，如圖11所示，誤差分布在0.31～0.78 m之間，中誤差為0.173 4 m。

圖11 相關性分析Fig.11 Correlation analysis

4 結論

本文針對遙感水邊線法潮灘DEM構建中，應用單一方法提取多期水邊線精度不高的問題，提出了基于多算法水邊線提取的潮灘DEM構建方法，并以萊州灣西側潮灘為例驗證了方法的有效性。主要結論如下：

1)分析了不同潮情下水邊線的邊緣類型與灘面噪聲的差異，并基于潮汐數據對水邊線進行分類，將研究數據中的多期水邊線分為： 單一灘面弱邊緣水邊線、復雜灘面弱邊緣水邊線、單一灘面強邊緣水邊線與復雜灘面強邊緣水邊線4種類型，該結果可進一步用于水邊線提取與潮灘DEM構建。

2)基于多算法提取水邊線，通過潮汐預測數據插值獲得多期水邊線瞬時潮位，并根據水邊線類別進一步篩選，構建了潮灘DEM。研究區DEM反演結果精度較高，可準確反映潮灘的近似地形，方法可有效用于潮灘地形信息的獲取。

但本文所選研究區與研究數據的代表性有限，可進一步完善水邊線分類規則，以適用于其他區域。通過潮汐預測數據插值計算瞬時潮位與水邊線實際潮位存在一定偏差，使用實測潮位數據可進一步提高研究區潮灘DEM的反演精度。