王子赫,康彥彥,王敏京

(1. 河海大學海洋學院,江蘇 南京 210098; 2. 自然資源部華東海岸帶野外科學觀測研究站,江蘇 南京 210007)

潮灘是由潮間帶淤泥帶部分與海底淺灘共同形成的沿海特殊地形,是海岸帶的重要組成部分,受潮汐的作用,高潮時被水淹沒,低潮時露出水面或變淺[1]。作為海陸相互作用的前沿地帶,近岸潮灘在海岸災害防護、土地資源利用及保護海洋生態(tài)環(huán)境等方面發(fā)揮著重要作用;但受風浪、潮流、風暴潮等因素的影響,其位置、形狀等不斷變化,這對沿海地形的監(jiān)測與更新提出了較高要求,因此對高精度潮灘地形反演方法的研究具有重要意義[2]。

考慮到潮灘普遍惡劣的自然環(huán)境,大規(guī)模的人工跑灘現(xiàn)場測量存在極大危險性和困難性,需要消耗大量人力物力[3]。而遙感測量技術憑借其時效性強、覆蓋范圍廣、經(jīng)濟耗費少等優(yōu)勢,受到大量學者的關注。InSAR與LiDAR技術的發(fā)展對于提取潮灘敏感地形信息、描述潮灘的形態(tài)特征具有很大潛力[4-5]。但受潮灘沉積物和殘留水體的影響,輻射信號受阻減弱,精度受限。同時InSAR技術的研究機理復雜、數(shù)據(jù)處理難度較大,很少被應用于實際的潮灘反演處理中[6]。雖然機載LiDAR技術的效率和精度都可以保證,但其受自然因素影響會產(chǎn)生數(shù)據(jù)缺失,同時獲取數(shù)據(jù)成本過高,對于中等精度高頻次大面積、連續(xù)、長期的更新監(jiān)測并不經(jīng)濟可行。

在眾多遙感監(jiān)測方法中,水邊線法以其適應性強、成本低及支持歷史地形構建等優(yōu)勢,成為目前潮灘地形監(jiān)測最具操作性的方法[7-8]。文獻[9]將水邊線視作等高線,利用多潮位賦值生成潮灘DEM模型。文獻[10]通過潮位測站基線計算水邊線高度值,構建了江蘇沿海輻射沙脊群DEM。文獻[11]利用二維水動力潮汐數(shù)值模型,模擬水邊線瞬時潮位,將高度值賦值于水線對應點,建立了江蘇沿海大面積潮灘歷史地形反演。水邊線法依托于潮位資料,賦值精度難以保證,同時在提取水邊線時存在必然的時間跨度,因此水邊線法反演地形存在局部區(qū)域差異。

2018年9月ICESat-2(ice,cloud and land elevation Satellite-2)激光衛(wèi)星成功發(fā)射,為潮灘地形構建帶來新的技術手段。該衛(wèi)星配備了高分辨率高級地形激光高度計系統(tǒng)(ATLAS),用于測量地球表面高程[12]。文獻[13]利用ICESat-2點云數(shù)據(jù)提高了航天飛機雷達地形任務(SRTM)高程產(chǎn)品衍生DEM的精度。文獻[14]通過構建ICESat-2高度數(shù)據(jù)與研究區(qū)域的淹沒頻率關系模型,獲得潮灘地形。但同一淹沒頻率下的實際高度同樣存在誤差,且使用遙感影像的時間跨度會更長,僅根據(jù)淹沒頻率很難反演真實的潮灘地形情況。

本文通過提取ICESat-2激光點云剖面數(shù)據(jù)高度值,結合多源光學影像與潮位進行反距離權重插值,提高水邊線賦值精度,得到一種新的潮灘地形反演方法,并與傳統(tǒng)水邊法DEM進行對比分析。

1 研究區(qū)與數(shù)據(jù)

1.1 研究區(qū)概況

選擇江蘇岸外輻射沙脊群南側水域的兩處大型沙洲——條子泥與高泥為研究對象。條子泥位于輻射沙脊群的中心地帶,東起東大港,西側與高泥相鄰,寬度約100 km,面積分別約為504.9、350 km2(如圖1所示)。受西大港、東大港、條魚港潮波系統(tǒng)及黃沙洋水域的共同作用,該區(qū)域泥沙交換活躍,潮灘地貌復雜多變,最大潮差超過9 m,具有輻射沙脊群水域最寬闊的潮面,灘面寬約14 km[15]。潮灘系統(tǒng)十分脆弱,其形態(tài)和動力地貌過程極易受人類活動的影響[16]。條子泥墾區(qū)一期建設總面積為6 095.62 hm2,大規(guī)模圍填海活動后,條子泥及其毗鄰近岸沉積動力環(huán)境的均衡態(tài)被打破,新的地貌格局逐漸形成[17]。

1.2 數(shù)據(jù)及預處理

研究初期收集2021年下半年的42組衛(wèi)星數(shù)據(jù),主要包括Landsat 8、GF-1、GF-6、HJ-1A衛(wèi)星,經(jīng)篩選獲取了10景清晰且含云量低于40%的光學圖像,其空間分辨率、幅寬等參數(shù)為適合水邊線法的應用,見表1。

本文使用的激光點云剖面數(shù)據(jù)是由ICESat-2衛(wèi)星搭載的ATLAS測高系統(tǒng)探測所得的。ATLAS是一套非掃描推帚式、高重頻、微脈沖、光子計數(shù)式激光測高系統(tǒng)。其沿軌道有3束激光束,相鄰光束之間的交叉軌道距離為3.3 km。每束激光被分為2個子束(兩子束激光距離較近,為90 m),其中一子束激光脈沖能量較強,約是較弱子束的4倍。地面探測波束結構如圖2所示,6束激光沿飛行方向同時向地面進行推掃[18]。

圖2 ICESat-2衛(wèi)星工作示意

激光點云高度數(shù)據(jù)采用NSIDC提供的ICESat-2二級產(chǎn)品——全球地理定位光子數(shù)據(jù)(ATL03),選取2021-07-09與2021-12-18兩組共12條激光軌道數(shù)據(jù)(如圖1所示)。其中黑色實線為2021-07-09軌道,由左至右分別為gt3r、gt3l、gt2r、gt2l、gt1r、gt1l的軌道;白色實線為2021-12-18軌道,由于ICESat-2在航行一段時間后會顛倒航線方向,因此該航道順序相反。原則上,在對數(shù)據(jù)進行篩選構建模型時,應選取軌道覆蓋潮灘面更為廣泛、能夠與水邊線產(chǎn)生更多交叉點,且激光脈沖能量較強的數(shù)據(jù)。因此選取2021-07-09中gt1l與2021-12-18中gt3r的高度值數(shù)據(jù)用于驗證,其余軌道數(shù)據(jù)用于構建DEM模型。

圖1中分別標注了東大港(DDG)、新條魚港(TYG)、洋口港(YKG)、大豐港(DFG)、東沙港(DSG)5處圍繞整個江蘇近海范圍的潮汐觀測站,能夠提供2021全年潮位數(shù)據(jù)?？堪兑粋葹镚PS RTK測量儀探測到的2021年10月的2組斷面實測高度數(shù)據(jù),每組數(shù)據(jù)按照約20 m的漸進步長獲取高度值,用于DEM模型的驗證。

2 研究方法

2.1 ICESat-2斷面數(shù)據(jù)處理

根據(jù)強光波束對不同表面散射回光束的能量參數(shù)值的不同,ATL03產(chǎn)品將目標地物分類為陸冰、海冰、海洋、陸地、內(nèi)陸水。產(chǎn)品根據(jù)每一光子事件的置信度值,將光子分為可能的信號光子點或背景光子點,通過信噪比與閾值的關系對光子事件進行歸類,以0、1、2、3、4為編號,分別對應噪聲、占位,以及低、中、高的置信度信號[12]。在進一步去噪前,選取高置信度信號(編號4)中被分類為陸地的信號光子事件,剔除其余光子事件,完成初步降噪,對應圖3中的降噪過程。

圖3 激光點云高度值提取過程

采用低通濾波窗口對其余光子點云剖面進行迭代運算。濾波降噪的主要思想是將信號中特定波段頻率濾除,是抑制和防止干擾的一項重要措施。本文采用的低通濾波能夠通過動態(tài)變化的方式改變頻率數(shù)值,且能夠適應不斷變化的潮灘地形高度,從而對激光點云數(shù)據(jù)進行降噪處理。其平滑與降噪的效果主要取決于截止頻率的設定,即

Wc=2·fc/fs

(1)

式中,Wc為3 dB截止頻率;fc為邊界頻率;fs為采樣頻率。去除背景噪聲后,激光點云剖面數(shù)據(jù)獲得最大限度的平滑,高度剖面線起伏變化位置周邊的噪聲得到削弱。對應圖3中的第4條降噪結果。

獲得最終潮灘高度數(shù)據(jù)前,需注意兩點:

(1)考慮到選取的光學影像分辨率為30 m,并以此基準提取水邊線數(shù)據(jù),為簡化數(shù)值計算過程、方便處理數(shù)據(jù);同時,為反映真實的地面高度情況,將濾波后的激光點云數(shù)據(jù)根據(jù)衛(wèi)星航行的距離(即光子事件間隔),選取每隔30 m范圍內(nèi)的光子事件為一組,提取每組數(shù)據(jù)中的高度中值,作為該區(qū)域的地形高度值。

(2)ATLAS后向散射的綠色激光束對清澈的Ⅰ類水體具有良好的穿透性,但對于輻射沙脊群水域渾濁的Ⅱ類水體,透射率極低,因此激光衛(wèi)星僅能獲取條子泥地區(qū)陸地與海表面的高度信息。而條子泥水域受潮汐作用的影響,海表面信息在不斷變化,本文選取的兩組ICESat-2衛(wèi)星過境軌道數(shù)據(jù)獲取時刻的潮位值均略高于該時間范圍的最低水位。因此去除多余的水平面數(shù)據(jù),獲得潮灘斷面的高度值。如圖3中最后獲得的潮灘地形高度數(shù)據(jù)。

2.2 批量獲取光學影像水邊線圖像數(shù)據(jù)

采用最大似然監(jiān)督分類方法對光學影像進行批量提取。提取前首先將光學遙感影像進行必要的預處理,包括輻射定標、大氣校正、幾何校正等。由于創(chuàng)建樣本時,同一類地物選擇過多將導致分類精度下降,因此應盡量選擇特征明顯的區(qū)域,如圖4(a)所示,將影像分為潮灘、水體兩類地物且作為訓練數(shù)據(jù)。使用最大似然法監(jiān)督分類選用多個波段對圖像進行分類,再經(jīng)濾波、聚類分析等分類后處理操作,得到圖4(b)。圖中高濁度懸沙或淺灘表層殘留水體等復雜地物,需結合潮汐、泥沙運動等因素,目視解譯對影像進行修改,提取影像中的水邊線數(shù)據(jù),結果如圖4(c)所示。白色邊緣線即為水邊線,線條連續(xù)清晰,能夠較好地展現(xiàn)影像拍攝時刻下潮灘的整體形態(tài)。采用提取的潮灘地形高度數(shù)據(jù)與序列水邊線在二維平面上進行疊加,產(chǎn)生多處交叉點,如圖4(d)所示。

圖4 最大似然法監(jiān)督分類批量水邊線提取

2.3 構建潮灘DEM模型

考慮到水邊線與激光剖面軌道的交點集中于激光衛(wèi)星軌道一側,其位置和數(shù)值均與距離較遠的潮位值存在一定差距。因此選用反距離權重(IDW)的方式構建海面并高度對水邊線進行賦值。該方法對于距預測位置最近的測量值對預測值的影響將會更大。假定每個測量點都存在局部影響,而這一影響會隨著距離的增大而減小[19]。

將圖像采集時刻對應的潮位與高度交點數(shù)據(jù)進行IDW計算,得到覆蓋整個研究區(qū)域范圍的柵格高度插值數(shù)據(jù),將結果賦值于對應水邊線中。為進一步提高模型精度,選取地形高度值中較為顯著的極值點與水邊線共同參與插值,這些點未能與水邊線相交,但能夠展現(xiàn)潮灘地形的部分細節(jié),作為補充修正。利用自然鄰域插值法獲得條子泥和高泥 2021年下半年的DEM。

本文使用的ICESat-2數(shù)據(jù),捕捉灘面的高度值范圍為-3～2 m。受潮汐水位的制約,選取最低潮位以上部分為掩膜興趣區(qū),提取該部分DEM后進行分析。

3 結果分析與討論

將預留的兩條軌道用于驗證模型的效果,RMSE分別為0.27、0.29 m。由此可知,IDW算法將潮灘地形高度值與潮汐水位值進行插值,獲得較準確的柵格高度值,模型與ICESat-2地形高度數(shù)據(jù)結合度較好,軌道處的細節(jié)得到了很好的展現(xiàn)。

圖5分別為傳統(tǒng)的潮位賦值水邊線法與本文基于ICESat-2高度值優(yōu)化的DEM模型。對比兩種DEM模型可知,優(yōu)化后的DEM模型展示了條子泥內(nèi)部更多的地形變化,深色部分主要集中于靠近江蘇近岸一側和沙洲系統(tǒng)的中心位置,淺色部分集中于沙洲系統(tǒng)邊緣的較低水位區(qū)域;經(jīng)過地形高度數(shù)值優(yōu)化后的模型能夠較好地識別出復雜多變的潮流、潮溝結構。

圖5 2021年下半年兩種方法獲得的DEM

為評估ICESat-2高度值對模型的優(yōu)化效果,將兩種方式構建的DEM模型結果與兩組實測斷面地形進行對比,如圖6所示。本文方法與傳統(tǒng)潮位賦值水邊線法獲得的DEM模型相比,R2由0.86提高至0.89,RMSE由0.45 m提高至0.34 m。兩種模型結果均可觀測到清晰的地形起伏變化,并與實測數(shù)據(jù)大致吻合。但在潮位較高的區(qū)域,僅依靠潮位的賦值無法準確反映潮灘上實際的地形。ICESat-2數(shù)據(jù)的輸入提高了賦值水邊線的精度,使得地形變化的細節(jié)被捕捉,進而顯現(xiàn)在模型中。模型與實測斷面存在的微小差異可能是由數(shù)據(jù)不同的采集日期導致的。因此,受研究區(qū)域廣闊、復雜且多變的地形,以及潮灘地區(qū)本就存在的時空動態(tài)變化等因素的影響,構建模型的數(shù)據(jù)應該盡可能控制在更短的時間跨度內(nèi)。

圖6 DEM結果誤差對比

4 結 語

本文通過對ICESat-2激光點云剖面數(shù)據(jù)進行低通濾波的降噪處理,提取了潮灘地形高度數(shù)據(jù),并與潮汐水位數(shù)據(jù)相結合,使用反距離權重的插值,進而提高水邊線賦值的精度,構建了新的DEM模型?；谲壍榔拭娣秶母叨戎祵φw潮灘地形進行優(yōu)化,實現(xiàn)了從二維到三維的轉換,結果符合潮灘地形基本規(guī)律,進一步驗證了ICESat-2數(shù)據(jù)應用于潮灘地形反演的可行性,為構建人類活動影響下的潮灘演化理論奠定了數(shù)據(jù)基礎,并為濕地生態(tài)環(huán)境保護和海岸防災減災提供了科學支撐。然而,目前的研究仍存在一定的局限性,如ICESat-2發(fā)射的綠色激光對于輻射沙脊群水域渾濁水體透射率極低,因此激光衛(wèi)星無法獲得水面以下的水深信息,潮位情況不可避免地阻礙衛(wèi)星獲取地形數(shù)據(jù)。但對于清澈水體綠色激光具有良好的穿透性,展現(xiàn)了ICESat-2在潮灘地形反演、水深勘探等領域的巨大潛力。