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

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

0 引言

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

遙感水邊線是衛(wèi)星過境時(shí)獲取的瞬時(shí)水邊線，準(zhǔn)確提取水邊線是潮灘DEM構(gòu)建的關(guān)鍵?，F(xiàn)有的水邊線提取方法主要有： 邊緣檢測(cè)法[8]、閾值分割法[9]、區(qū)域生長法[10]、主動(dòng)輪廓模型法[11]、面向?qū)ο蠓诸怺12]等。Mason等[13]基于ERS SAR影像，通過紋理分割提取多期水邊線，應(yīng)用水動(dòng)力模型對(duì)水邊線高程進(jìn)行賦值，構(gòu)建了英國東岸Humber/Wash區(qū)域潮灘DEM； 沈芳等[14]基于Landsat TM影像，對(duì)比了不同波段的水邊線提取精度，利用閾值分割法提取水邊線，通過潮汐預(yù)測(cè)數(shù)據(jù)插值對(duì)水邊線賦值，構(gòu)建了長江口九段沙DEM； 穆敬等[15]基于BJ-1影像，利用面向?qū)ο蠓诸惙ㄌ崛∷吘€，通過潮汐網(wǎng)格數(shù)據(jù)對(duì)水邊線賦值，構(gòu)建了黃驊市潮灘DEM。目前國內(nèi)外學(xué)者在構(gòu)建潮灘DEM的研究中大多關(guān)注于準(zhǔn)確模擬水邊線的瞬時(shí)潮位和提升單一方法的水邊線提取精度[8-15]。但是，瞬時(shí)水邊線受衛(wèi)星過境時(shí)的潮情、天氣以及灘面等因素的影響，在影像中的光譜與紋理差異較大，基于單一方法提取水邊線，難以保證水邊線的提取精度，從而影響潮灘DEM反演結(jié)果。吳迪等[16]雖提出了基于多源多算法的水邊線提取模型，但未對(duì)水邊線進(jìn)行系統(tǒng)劃分，且沒有將其應(yīng)用于潮灘DEM構(gòu)建研究。

針對(duì)上述問題，本文從邊緣類型與灘面噪聲2方面劃分水邊線，根據(jù)瞬時(shí)潮情對(duì)水邊線進(jìn)行分類提取，通過潮汐數(shù)據(jù)構(gòu)建多項(xiàng)式計(jì)算瞬時(shí)潮位，提出基于多算法水邊線提取的潮灘DEM構(gòu)建方法，并通過野外實(shí)測(cè)數(shù)據(jù)對(duì)潮灘DEM反演精度進(jìn)行驗(yàn)證。

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

1.1 研究區(qū)概況

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

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

1.2 研究數(shù)據(jù)源及預(yù)處理

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

由于受野外測(cè)量條件限制，缺少研究區(qū)的實(shí)測(cè)潮位數(shù)據(jù)，本文使用濰坊港驗(yàn)潮站的潮汐預(yù)測(cè)數(shù)據(jù)用于水邊線的分類與衛(wèi)星過境瞬時(shí)的水邊線潮位計(jì)算。潮汐預(yù)測(cè)數(shù)據(jù)是在潮汐觀測(cè)數(shù)據(jù)的基礎(chǔ)上，利用調(diào)和分析法推算的預(yù)測(cè)潮位[14]。濰坊港驗(yàn)潮站所在海域?yàn)椴灰?guī)則半日混合潮，潮汐預(yù)測(cè)數(shù)據(jù)包含每日整時(shí)的潮位值，以及每日的潮位極值與對(duì)應(yīng)時(shí)刻，潮高基準(zhǔn)面為平均海平面以下120 cm。

2 多算法水邊線提取

2.1 基于潮情的水邊線分類

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

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

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

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

圖2 分類規(guī)則

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

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

表1 水邊線分類結(jié)果Tab.1 Waterline classification results

2.2 基于邊緣檢測(cè)水邊線提取

邊緣檢測(cè)是基于圖像灰度的一階或二階導(dǎo)數(shù)來提取圖像中灰度急劇變化的區(qū)域邊界[18]，具有模型參數(shù)少、運(yùn)行速度快的優(yōu)點(diǎn)，但對(duì)噪聲敏感，在復(fù)雜灘面提取的水邊線存在嚴(yán)重的毛刺或不連續(xù)。本文選擇定位準(zhǔn)確度高、單邊響應(yīng)良好的Canny邊緣檢測(cè)算法[19]，以20170429影像為例，對(duì)單一灘面強(qiáng)邊緣水邊線進(jìn)行提取，結(jié)果如圖3所示。

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

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

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

以20170617影像為例，首先，計(jì)算NDWI，增強(qiáng)影像中的水體信息，計(jì)算公式為：

(1)

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

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

2.4 基于面向?qū)ο蠓ㄋ吘€提取

面向?qū)ο蠓ㄊ峭ㄟ^影像中的光譜信息、空間信息和紋理信息對(duì)影像進(jìn)行分割，然后以分割后的對(duì)象為基本單元進(jìn)行分類，從而提取水邊線的一種方法[22]。面向?qū)ο蠓ň哂休^強(qiáng)的抗噪能力，有效抑制復(fù)雜灘面上其他地物的影響，可準(zhǔn)確保留強(qiáng)邊緣水邊線的細(xì)節(jié)信息，但該方法對(duì)紋理變化過于敏感，提取的弱邊緣水邊線不平滑，因此，本文通過面向?qū)ο蠓ㄌ崛?fù)雜灘面強(qiáng)邊緣水邊線。

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

2.5 基于改進(jìn)分水嶺水邊線提取

圖6 改進(jìn)分水嶺算法流程Fig.6 Flow chart of improved watershed algorithm

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

圖7 基于改進(jìn)分水嶺的水邊線提取結(jié)果Fig.7 Waterline extraction result basedon improved watershed algorithm

3 潮灘DEM構(gòu)建與精度驗(yàn)證

3.1 潮灘DEM構(gòu)建

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

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

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

圖8 水邊線選用結(jié)果Fig.8 Selection result of waterlines

將篩選后的水邊線重采樣為16 m×16 m的水邊點(diǎn)，利用不規(guī)則三角網(wǎng)進(jìn)行空間插值，并將插值結(jié)果從潮高基準(zhǔn)面轉(zhuǎn)換為以WGS84為基準(zhǔn)的高程數(shù)據(jù)，得到潮灘DEM反演結(jié)果，如圖9所示。

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

3.2 精度驗(yàn)證

為驗(yàn)證基于多算法水邊線提取的潮灘DEM構(gòu)建方法有效性，本文利用野外實(shí)測(cè)高程數(shù)據(jù)對(duì)潮灘DEM反演結(jié)果進(jìn)行精度評(píng)價(jià)。實(shí)測(cè)高程數(shù)據(jù)獲取時(shí)間為2017年12月，采用Trimble R5 GPS接收機(jī)測(cè)量，靜態(tài)觀測(cè)不少于20分鐘。實(shí)測(cè)高程數(shù)據(jù)通過精密星歷解算為WGS84坐標(biāo)數(shù)據(jù)，18個(gè)驗(yàn)證點(diǎn)分布如圖10所示。

圖10 驗(yàn)證點(diǎn)分布Fig.10 Distribution of verification points

對(duì)實(shí)測(cè)高程與反演結(jié)果進(jìn)行相關(guān)分析，實(shí)測(cè)高程與對(duì)應(yīng)DEM反演結(jié)果的相關(guān)性較高，R2為0.864 9，如圖11所示，誤差分布在0.31～0.78 m之間，中誤差為0.173 4 m。

圖11 相關(guān)性分析Fig.11 Correlation analysis

4 結(jié)論

本文針對(duì)遙感水邊線法潮灘DEM構(gòu)建中，應(yīng)用單一方法提取多期水邊線精度不高的問題，提出了基于多算法水邊線提取的潮灘DEM構(gòu)建方法，并以萊州灣西側(cè)潮灘為例驗(yàn)證了方法的有效性。主要結(jié)論如下：

1)分析了不同潮情下水邊線的邊緣類型與灘面噪聲的差異，并基于潮汐數(shù)據(jù)對(duì)水邊線進(jìn)行分類，將研究數(shù)據(jù)中的多期水邊線分為： 單一灘面弱邊緣水邊線、復(fù)雜灘面弱邊緣水邊線、單一灘面強(qiáng)邊緣水邊線與復(fù)雜灘面強(qiáng)邊緣水邊線4種類型，該結(jié)果可進(jìn)一步用于水邊線提取與潮灘DEM構(gòu)建。

2)基于多算法提取水邊線，通過潮汐預(yù)測(cè)數(shù)據(jù)插值獲得多期水邊線瞬時(shí)潮位，并根據(jù)水邊線類別進(jìn)一步篩選，構(gòu)建了潮灘DEM。研究區(qū)DEM反演結(jié)果精度較高，可準(zhǔn)確反映潮灘的近似地形，方法可有效用于潮灘地形信息的獲取。

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