李天祺,王建超,吳芳,趙政,張文凱
(中國自然資源航空物探遙感中心,北京 100083)
淤泥質(zhì)潮灘通常指平均大潮高潮線與平均大潮低潮線之間的地帶,也稱潮間帶[1-2]。在我國約有四分之一的海岸屬于淤泥質(zhì)海岸[3],獲取淤泥質(zhì)潮灘地形對研究岸線變遷、海岸帶生態(tài)變化及沿海工程建設有著重要意義[4-5]。淤泥質(zhì)潮灘受波浪、泥沙、沿岸流和地質(zhì)等諸多因素的影響,具有面積寬廣、水淺灘平、變化頻繁的特點,使得該區(qū)域地面調(diào)查與地形測繪的難度較大、成本較高,且存在一定風險[6-7]。遙感技術具有高時效、大范圍、高頻率的特點,利用多期遙感水邊線復合潮位數(shù)據(jù)構(gòu)建潮灘數(shù)字高程模型(digital elevation model,DEM),成為獲取大范圍淤泥質(zhì)潮灘地形信息的有效途徑。
遙感水邊線是衛(wèi)星過境時獲取的瞬時水邊線,準確提取水邊線是潮灘DEM構(gòu)建的關鍵?,F(xiàn)有的水邊線提取方法主要有: 邊緣檢測法[8]、閾值分割法[9]、區(qū)域生長法[10]、主動輪廓模型法[11]、面向?qū)ο蠓诸怺12]等。Mason等[13]基于ERS SAR影像,通過紋理分割提取多期水邊線,應用水動力模型對水邊線高程進行賦值,構(gòu)建了英國東岸Humber/Wash區(qū)域潮灘DEM; 沈芳等[14]基于Landsat TM影像,對比了不同波段的水邊線提取精度,利用閾值分割法提取水邊線,通過潮汐預測數(shù)據(jù)插值對水邊線賦值,構(gòu)建了長江口九段沙DEM; 穆敬等[15]基于BJ-1影像,利用面向?qū)ο蠓诸惙ㄌ崛∷吘€,通過潮汐網(wǎng)格數(shù)據(jù)對水邊線賦值,構(gòu)建了黃驊市潮灘DEM。目前國內(nèi)外學者在構(gòu)建潮灘DEM的研究中大多關注于準確模擬水邊線的瞬時潮位和提升單一方法的水邊線提取精度[8-15]。但是,瞬時水邊線受衛(wèi)星過境時的潮情、天氣以及灘面等因素的影響,在影像中的光譜與紋理差異較大,基于單一方法提取水邊線,難以保證水邊線的提取精度,從而影響潮灘DEM反演結(jié)果。吳迪等[16]雖提出了基于多源多算法的水邊線提取模型,但未對水邊線進行系統(tǒng)劃分,且沒有將其應用于潮灘DEM構(gòu)建研究。
針對上述問題,本文從邊緣類型與灘面噪聲2方面劃分水邊線,根據(jù)瞬時潮情對水邊線進行分類提取,通過潮汐數(shù)據(jù)構(gòu)建多項式計算瞬時潮位,提出基于多算法水邊線提取的潮灘DEM構(gòu)建方法,并通過野外實測數(shù)據(jù)對潮灘DEM反演精度進行驗證。
研究區(qū)位于N37°21.5′~37°32.0′,E118°53.5′~118°58.5′,位于萊州灣西側(cè),濰坊港驗潮站以北,如圖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
本文選用GF-1 WFV數(shù)據(jù),空間分辨率為16 m,重訪周期為2 d,幅寬為800 km。數(shù)據(jù)獲取時間為2017年2—7月,篩除受厚云、雪及海冰污染的數(shù)據(jù),共選取了18期數(shù)據(jù)用于研究。數(shù)據(jù)預處理包括絕對輻射定標、幾何精糾正及研究區(qū)裁剪,其中幾何精糾正誤差控制在0.5個像元以內(nèi)。
由于受野外測量條件限制,缺少研究區(qū)的實測潮位數(shù)據(jù),本文使用濰坊港驗潮站的潮汐預測數(shù)據(jù)用于水邊線的分類與衛(wèi)星過境瞬時的水邊線潮位計算。潮汐預測數(shù)據(jù)是在潮汐觀測數(shù)據(jù)的基礎上,利用調(diào)和分析法推算的預測潮位[14]。濰坊港驗潮站所在海域為不規(guī)則半日混合潮,潮汐預測數(shù)據(jù)包含每日整時的潮位值,以及每日的潮位極值與對應時刻,潮高基準面為平均海平面以下120 cm。
水邊線作為遙感影像中一種重要的邊緣信息,是水陸區(qū)域的分界。不同遙感水邊線提取算法具有不同的適用范圍[17],其提取效果主要取決于目標邊緣的類型與邊緣鄰域內(nèi)的噪聲水平。根據(jù)邊緣類型與噪聲水平對水邊線分類,選擇適用算法進行提取,有利于準確地獲取水邊線信息。
水邊線兩側(cè)水陸區(qū)域的灰度變化梯度主要受漲落潮的影響: 漲潮時灘面干濕分明,水邊線兩側(cè)區(qū)域灰度差異顯著; 落潮時受灘面殘余水和水中泥沙影響,水邊線兩側(cè)灰度差異較小。因此按灰度變化梯度,邊緣類型可分為強邊緣與弱邊緣。水邊線的鄰域噪聲水平,主要由潮位高低決定: 由于灘面地物種類由近海向內(nèi)陸逐漸增多,高潮位時水邊線在潮溝、植被、建筑物分布較多的近岸區(qū)域,灘面破碎復雜,鄰域內(nèi)噪聲較高; 低潮位時水邊線近海,灘面相對均質(zhì)單一。
綜上,本文通過分析潮汐數(shù)據(jù),確定各期水邊線的瞬時潮情: 漲、落潮與高、低潮位,根據(jù)瞬時潮情將水邊線劃分為單一灘面弱邊緣、復雜灘面弱邊緣、單一灘面強邊緣和復雜灘面強邊緣4類,具體分類規(guī)則如圖2所示。
Test and optimization on heat transfer performance of electric heater under lower air pressure
圖2 分類規(guī)則
在本實驗中,影像中的薄云或霧明顯削弱灰度變化梯度,因此通過人工目視對影像質(zhì)量進行判斷,將包含薄云或霧影像的水邊線分為弱邊緣水邊線; 通過分析潮汐數(shù)據(jù)對晴朗影像的水邊線進行劃分。若H(10),H(11)和H(12)3個時刻的潮位依次升高,則瞬時潮情(11—12時之間)為漲潮,灘面不存在殘余水的影響,則水邊線劃分為強邊緣水邊線,若3個時刻的潮位依次下降(落潮)或10時的潮位高于11時潮位,則存在灘面殘余水或水中泥沙的影響,劃分為弱邊緣水邊線。
獲得強、弱邊緣水邊線后,按潮位高低對水邊線所在灘面特征分類: 若水邊線獲取當日的11時與12時的平均潮位高于多期平均值H,則該期水邊線近岸,對應復雜灘面; 若低于多期平均值H則該期水邊線近海,對應單一灘面。最后得到4類水邊線分別為: 單一灘面弱邊緣水邊線、復雜灘面弱邊緣水邊線、單一灘面強邊緣水邊線與復雜灘面強邊緣水邊線。對本文所用的18期GF-1 WFV影像的水邊線進行分類,分類結(jié)果如表1所示。
表1 水邊線分類結(jié)果Tab.1 Waterline classification results
邊緣檢測是基于圖像灰度的一階或二階導數(shù)來提取圖像中灰度急劇變化的區(qū)域邊界[18],具有模型參數(shù)少、運行速度快的優(yōu)點,但對噪聲敏感,在復雜灘面提取的水邊線存在嚴重的毛刺或不連續(xù)。本文選擇定位準確度高、單邊響應良好的Canny邊緣檢測算法[19],以20170429影像為例,對單一灘面強邊緣水邊線進行提取,結(jié)果如圖3所示。
圖3 基于Canny邊緣檢測的水邊線提取結(jié)果Fig.3 Waterline extraction resultbased on Canny edge detection
由于弱邊緣水邊線兩側(cè)灰度變化不明顯,邊緣檢測不能有效提取邊緣信息,閾值分割通過設定閾值,如歸一化水體指數(shù)(normalized different water index,NDWI)[20],可以快速有效地區(qū)分水陸區(qū)域,但對于高潮位的復雜灘面,受潮溝或其他積水區(qū)域的影響,水邊線提取不準確,因此,本文通過閾值分割提取單一灘面弱邊緣水邊線。
以20170617影像為例,首先,計算NDWI,增強影像中的水體信息,計算公式為:
(1)
式中Green與NIR分別代表GF-1 WFV影像中的綠光波段與近紅外波段的反射率。NDWI的取值范圍為[-1, 1]; 然后,通過OTSU算法[21]計算閾值,對NDWI灰度影像進行二值化; 最后,對二值圖做開運算處理,對邊界進行平滑,得到水邊線提取結(jié)果,如圖4所示。
圖4 基于NDWI閾值分割的水邊線提取結(jié)果Fig.4 Waterline extraction result based onNDWI threshold segmentation
面向?qū)ο蠓ㄊ峭ㄟ^影像中的光譜信息、空間信息和紋理信息對影像進行分割,然后以分割后的對象為基本單元進行分類,從而提取水邊線的一種方法[22]。面向?qū)ο蠓ň哂休^強的抗噪能力,有效抑制復雜灘面上其他地物的影響,可準確保留強邊緣水邊線的細節(jié)信息,但該方法對紋理變化過于敏感,提取的弱邊緣水邊線不平滑,因此,本文通過面向?qū)ο蠓ㄌ崛碗s灘面強邊緣水邊線。
首先,通過多尺度分割將影像分割為多個對象,然后根據(jù)各個對象的光譜相似度進行合并,多次迭代后最終得到水邊線提取結(jié)果,以20170609影像為例,初始分割結(jié)果和水邊線提取結(jié)果如圖5所示。
圖6 改進分水嶺算法流程Fig.6 Flow chart of improved watershed algorithm
以20170310影像為例,對影像NDWI低通濾波去除高頻噪聲后,潮灘上的潮溝、河口處的沙洲及較小的人工建筑物均已去除; 二值化處理后通過數(shù)學形態(tài)學處理,將鹽田、河流等內(nèi)陸水體完全去除,在邊界處做緩沖區(qū)處理,生成待分割區(qū)域的標記信息; 利用原始影像B4(R),B3(G),B2(B)假彩色圖像結(jié)合標記信息進行分水嶺分割,提取結(jié)果如圖7所示。
圖7 基于改進分水嶺的水邊線提取結(jié)果Fig.7 Waterline extraction result basedon improved watershed algorithm
利用濰坊港驗潮站的潮汐預測數(shù)據(jù)對影像獲取的瞬時潮位進行插值。本文所用GF-1 WFV數(shù)據(jù)的獲取時間均在每日上午11—12時之間,利用10—13時的4個整時潮位值構(gòu)建三次多項式,若在10—13時之間存在潮位極值,則用潮位極值替代距影像獲取時間較遠的整點潮位值,通過多項式插值計算水邊線的瞬時潮位。在研究資料有限時,這是一種行之有效的方法[14],水邊線對應瞬時潮位如表2所示。
表2 水邊線瞬時潮位Tab.2 Instantaneous tidal level of waterlines (cm)
受灘面上潮溝或河流入??谔幠嗌车挠绊懀崛〉亩嗥谒吘€在局部存在交叉或重疊的現(xiàn)象。針對上述問題,本文按照單一灘面強邊緣、復雜灘面強邊緣、單一灘面弱邊緣和復雜灘面弱邊緣的優(yōu)先順序選用水邊線,結(jié)果如圖8所示。
圖8 水邊線選用結(jié)果Fig.8 Selection result of waterlines
將篩選后的水邊線重采樣為16 m×16 m的水邊點,利用不規(guī)則三角網(wǎng)進行空間插值,并將插值結(jié)果從潮高基準面轉(zhuǎn)換為以WGS84為基準的高程數(shù)據(jù),得到潮灘DEM反演結(jié)果,如圖9所示。
圖9 潮灘DEM反演結(jié)果Fig.9 Inversion result of tidal flat DEM
為驗證基于多算法水邊線提取的潮灘DEM構(gòu)建方法有效性,本文利用野外實測高程數(shù)據(jù)對潮灘DEM反演結(jié)果進行精度評價。實測高程數(shù)據(jù)獲取時間為2017年12月,采用Trimble R5 GPS接收機測量,靜態(tài)觀測不少于20分鐘。實測高程數(shù)據(jù)通過精密星歷解算為WGS84坐標數(shù)據(jù),18個驗證點分布如圖10所示。
圖10 驗證點分布Fig.10 Distribution of verification points
對實測高程與反演結(jié)果進行相關分析,實測高程與對應DEM反演結(jié)果的相關性較高,R2為0.864 9,如圖11所示,誤差分布在0.31~0.78 m之間,中誤差為0.173 4 m。
圖11 相關性分析Fig.11 Correlation analysis
本文針對遙感水邊線法潮灘DEM構(gòu)建中,應用單一方法提取多期水邊線精度不高的問題,提出了基于多算法水邊線提取的潮灘DEM構(gòu)建方法,并以萊州灣西側(cè)潮灘為例驗證了方法的有效性。主要結(jié)論如下:
1)分析了不同潮情下水邊線的邊緣類型與灘面噪聲的差異,并基于潮汐數(shù)據(jù)對水邊線進行分類,將研究數(shù)據(jù)中的多期水邊線分為: 單一灘面弱邊緣水邊線、復雜灘面弱邊緣水邊線、單一灘面強邊緣水邊線與復雜灘面強邊緣水邊線4種類型,該結(jié)果可進一步用于水邊線提取與潮灘DEM構(gòu)建。
2)基于多算法提取水邊線,通過潮汐預測數(shù)據(jù)插值獲得多期水邊線瞬時潮位,并根據(jù)水邊線類別進一步篩選,構(gòu)建了潮灘DEM。研究區(qū)DEM反演結(jié)果精度較高,可準確反映潮灘的近似地形,方法可有效用于潮灘地形信息的獲取。
但本文所選研究區(qū)與研究數(shù)據(jù)的代表性有限,可進一步完善水邊線分類規(guī)則,以適用于其他區(qū)域。通過潮汐預測數(shù)據(jù)插值計算瞬時潮位與水邊線實際潮位存在一定偏差,使用實測潮位數(shù)據(jù)可進一步提高研究區(qū)潮灘DEM的反演精度。