Sentinel-2和ICESat-2密集時序數據反演淺平型內陸湖泊水量變化

吳浩儒,李均力,包安明,張久丹,馬英蓮

1. 中國科學院新疆生態與地理研究所荒漠與綠洲生態國家重點實驗室,新疆 烏魯木齊 830011; 2. 中國科學院大學,北京 100049; 3. 新疆遙感與地理信息系統應用重點實驗室,新疆 烏魯木齊 830011; 4. 新疆維吾爾自治區測繪科學研究院,新疆 烏魯木齊 830002

湖泊水量是湖泊水情的關鍵因素,相對于面積和水位,水量更能客觀地指示湖泊的時空變化[1]。干旱區內陸尾閭湖是封閉型流域水資源的最終匯集點,其水量變化能夠客觀反映流域的水量平衡過程[2]。然而,大部分干旱區內陸湖泊地處偏遠,水文基礎資料稀缺,尚未有湖盆地形及水量的公開信息[3],重建無資料地區內陸湖泊準確的水量變化可為分析內陸湖泊變化的機理、科學配置干旱區水資源及制定河湖生態治理策略提供科學依據[4]。

無資料或缺資料地區湖泊水量變化提取方法目前主要有基于地形測量和基于遙感的方法。其中,湖盆地形測量多采用船載回聲測深法[5]和RTK地形測量法[6],這種方法測量精度高,但作業成本高、測量效率低,對于水深變動大、季節性變化顯著的尾閭湖,開展相關測量工作的難度更大。遙感技術具有大尺度、周期性及全天候的觀測特性,目前已成為獲取內陸水體水文信息的重要手段[7],相關文獻對光譜特征法、重力衛星法及水位-面積關系法等基于遙感的水量反演方法進行了歸納和總結[8-9]。其中,光譜特征法通過建立湖泊水深與水體光譜特征的經驗關系進行水深反演[10-12],適用于水質良好的淺水湖;重力衛星法通過測量研究區域的重力變化來估算水量變化[13-15],適用于大尺度、大區域內陸流域的水量變化,如塔里木河流域或咸海流域;湖泊水位-面積關系法主要利用測高衛星和光學衛星分別獲取水位、面積信息,建立水位、面積與湖泊水量變化之間的關系,并根據面積或水位推算湖泊水量或水量變化[16-19],適用于有時間一致的水位和水面積監測數據的湖泊。

內陸干旱區尾閭湖泊是河流流至廣袤剝蝕盆地中心潴積而成,在長期湖泊沉積和強烈的風生湖流作用下,通常湖水較淺且湖底地形平坦,湖盆的最大深度往往只有4～5 m,湖岸線極其復雜[20]。這種湖泊的變化主要表現為水面變化迅速而水位變化微小,現有方法很難根據水位-水面關系來確定水量變化;而且在多風條件下,湖水中泥沙和懸浮物較多導致水體光譜特征不穩定,光譜特征法反演水深的精度變動較大[21]。近年來,文獻[22]利用湖泊水面的等高特性,獲取不同時相湖泊水面的激光雷達測高信息,形成若干等高線,并構建湖盆地形,在邊界規則、地形漸進變化的自然湖泊上取得了較好的效果[23-24]。

隨著遙感技術的不斷發展,衛星遙感分辨率、覆蓋周期和測高精度都有了長足的進步,其中Sentinel-2光學遙感衛星可實現5 d的重訪周期[25],ICESat-2的激光測高儀發射的激光脈沖可獲取沿軌道方向間隔0.7 m、光斑直徑為17 m的光子信號,測高精度達到了2～3 cm[26-27],可以滿足微小地形變化的測量要求。本文結合以上兩種衛星的優勢,以塔里木河干流尾閭湖臺特瑪湖為例,提出一種基于湖盆地形重建的水量變化估算方法。該方法通過湖盆內多時相ICESat-2非水面高程點和多時相Sentinel-2水面線等高點構建高精度湖盆地形信息,并在此基礎上利用時序面積信息重建水量變化時序。

1 研究方法

湖盆地形重建主要技術流程如圖1所示,主要包括多時序水面邊界線的提取、多時序湖盆激光雷達高程點提取,水面邊界線的高程內插、地形構建和水量時序重建等。其中構建高精度的湖盆地形是獲取淺平型內陸湖泊水文參數的關鍵。單光子計數激光雷達(photon-counting LiDAR,PCL)的光子高程點可獲取平坦地面高精度的高程信息,但由于軌道間距較大,無法全面覆蓋湖盆,單靠湖盆內的激光雷達高程點無法構建完整的湖盆地形信息。淺平型內陸湖泊在上游河流來水時面積迅速擴張,河流斷流時水面迅速消退乃至干涸。特別是在水面消退期間,湖泊分裂成眾多小水面,且水體流動少、水面相對靜止,可以近似將每個小水面看作一個等高面。為此,本文提出改進的湖盆反演方法,利用湖盆內無水區域的激光光子高程點對處于退水期的多時相水面邊界線進行高程插值,不同水面邊界線的高程點加密了湖盆內的高程點。在此基礎上,利用加密的高程點構建TIN和數字高程模型,并采用精確的湖盆地形和多時相水面邊界重建湖泊水量變化。

圖1 淺平型湖盆的地形構建技術流程Fig.1 Flowchart of basin topographic retrieval for shallow-flat lakes

1.1 湖泊水面邊界線提取

干旱區內陸湖泊的水陸邊界明顯,水體指數(normal different water index,NDWI)[28]是區分陸地和水面的一種主要方法。為此,本文利用光學遙感時序影像生產的NDWI圖,采用“全局-局部”自適應閾值的方法[29]提取湖泊邊界的時序信息。該方法提取的水面時序信息的生產者精度在96.42%～99.73%之間。

1.2 PCL光子高程點提取

PCL激光脈沖信號從發射到回波接收的過程中會受到太陽輻射、大氣散射等影響,形成較多的噪聲,如何去除噪聲光子會直接影響高程信息的反演精度[30]。為此,本文提出一種改進的密度自適應的信號光子檢測算法來區分信號與噪聲(圖2)。其主要原理是激光光子到達反射面并回波時,信號光子沿反射面的分布密度遠高于噪聲光子[31]。本文的信號與噪聲光子檢測方法如下。

圖2 PCL光子高程信息提取Fig.2 PCL photons elevation information extraction

(1) 建立橢圓鄰域計算光子密度。對每個激光光子m,根據光子信號與噪聲的分布規律[32],采用長軸a=10 m、短軸b=1 m、傾角θ=0°構建橢圓鄰域范圍(圖2(b)),根據以下公式判斷光子點n是否在中心光子點m的橢圓鄰域內

dx=(xm-xn)cosθ+(ym-yn)sinθ

(1)

dy=(ym-yn)cosθ-(xm-xn)sinθ

(2)

(3)

式中,a、b和θ分別表示橢圓鄰域的長軸、短軸和傾角;xm和ym分別為中心光子點m的緯度和海拔;xn和yn分別為光子點n的緯度和海拔。如果dist(m,n)<1則說明光子n為中心光子點m鄰域范圍內的光子點。

根據鄰域內的所有光子計算中心光子點m的光子密度如下

(4)

(5)

(6)

式中,Pm為中心光子點m的鄰域點集;N為Pm中光子點數量;Wn是采用高斯核計算的橢圓鄰域中每個光子點的權重[33];δ是高斯核的標準差;Dm為中心光子點m的點密度,由鄰域內所有點的高斯權重求和得到。

(2) 計算光子最大點密度。順時針每旋轉15°計算一次中心光子密度,最后取得到的最大密度作為中心光子的點密度。通過自適應調節橢圓鄰域的傾角θ,可以使信號光子的橢圓鄰域減小地形變化的影響,捕捉到更多密集光子,更容易與噪聲光子進行區分(圖2(a)和圖2(c))。

(3) 確定信號與噪聲光子的閾值。當所有光子的最大點密度計算完畢,采用OTSU閾值法[34]計算信號光子和噪聲光子密度之間的閾值,從而提取信號光子。

(4) 基于信號點生成地形剖面線。根據提取單激光束所有的信號光子點,采用最小二乘回歸方法擬合出該軌跡下的地形剖面線,對信號光子較少的區域進行加密。

在計算PCL的地面光子高程信息后,利用與之時間對應的湖泊圖層,提取湖盆范圍內、湖面邊界外的高程點數據,并將其高程基準面轉換為大地水準面,作為湖盆內無水區的高程點。對不同時相的所有高程點進行提取,最后合并為覆蓋湖盆的PCL光子高程點。

1.3 湖泊水面邊界線高程內插

在獲取湖盆內無水區的PCL光子高程信息點后,還需要確定不同時相湖面邊界線的高程值,從而加密湖盆內的PCL光子高程信息點。本文采用反距離加權(inverse distance weighted,IDW)[35]的內插方法,用PCL光子點的高程值內插湖面邊界線上點的高程值,進而獲取湖邊界線高程值,具體方法如圖3所示。對于每條湖面邊界線都會有若干條PCL光子軌跡與其相交(圖3(a)),考慮到臺特瑪湖湖盆較為平坦,且PCL光子軌跡上點較為密集,可以在軌跡上選取交點兩側各4個距離最近的湖盆高程點進行反距離加權插值得到交點高程(圖3(b))。對于每條湖面邊界線,可以求出邊界交點的高程值,選定遠離平均值3倍標準差的值為異常值并將其剔除[36],然后取平均作為水面邊界的高程。

圖3 內插湖泊邊界高程Fig.3 Altitude interpolation of lake boundaries

1.4 湖盆地形構建

多時相的湖泊水面邊界高程點既能加密湖盆內的PCL光子高程點,水面等深線也能更好模擬湖盆的地形地貌特征。為此,利用加密后的湖盆高程點構建不規則三角網,并在此基礎上利用Kriging插值生成湖盆的數字高程模型(DEM),DEM的格網大小與提取水面邊界線的遙感影像分辨率一致。

1.5 湖泊水量時序重建

淺平型湖泊的季節性變化顯著,湖泊往往分裂為眾多的小湖,且每個小湖泊的水面高程存在一定差異,常用的“水位-面積-水量”方法不能適用這種情況,需要分別對每個小湖分別計算水量,并累加為湖泊的總水量。為此,每個時相的水面對應的湖泊水量計算方法為

(7)

式中,Nt為時相t湖泊的小水面數量;Ai為第i個小水面的像素數量;Hi為第i個小水面的高程;Hi,j為第i個小水面、第j個像素的高程;Si,j為第i個小水面、第j個像素的面積。對每個時相的湖泊水面線,根據式(7)分別計算湖泊的水量,最后形成湖泊的水量時序信息。

2 研究區與數據

2.1 研究區概況

本文選擇塔里木河干流下游的尾閭湖——臺特瑪湖作為試驗區構建湖盆并重建湖泊水量,并且以最大湖面為基準,根據地勢起伏程度向外擴張3～5 km作為湖盆范圍,如圖4所示。臺特瑪湖位于新疆維吾爾自治區若羌縣的北部,是塔里木河干流和車爾臣河交匯形成的尾閭湖。臺特瑪湖的西部是車爾臣河下游的康拉克湖群,南部的若羌河與米蘭河是間歇河流,僅在豐水期有地表徑流匯入湖區。湖泊東部為庫木塔格沙漠,歷曾有河道與羅布泊相連,后因河水量減少及風沙阻塞河道,臺特瑪湖開始在車爾臣河與塔里木河交匯處潴水成湖,湖泊周邊地形平坦,地勢呈西高東低的特征。

2.2 數據源選擇

湖泊水面線信息提取采用2016—2022年共132期無云和質量較好的Sentinel-2衛星遙感數據,覆蓋了不同季節從湖盆無水到最大水面的各種情況。高程信息數據從冰、云和陸地高程衛星2號(the ice,cloud,and land elevation satellite-2,ICESat-2)的地形激光測高儀系統(advanced topographic laser altimeter system,ATLAS)數據提取,數據源為美國冰雪數據中心發布的ICESat-2/ATLAS全球地理定位光子數據產品ATL03[37]。ATL03產品是獲取陸地表面高程的基礎數據,記錄了包括大量噪聲光子在內的傳感器接收的所有光子經緯度及橢球體高度信息[38]。本文挑選了28對時相相近的ICESat-2測高點數據和Sentinel-2影像用于湖盆高程點提取,詳細信息見表1。

表1 選擇的ICESat-2和Sentinel-2數據對應列表

3 試驗與結果分析

3.1 試驗過程

本文基于ICESat-2測高衛星數據和Sentinel-2遙感影像數據,提取了密集的湖盆高程點,構建了較為精確的湖盆DEM,具體過程如下:①分別提取多時相Sentinel-2水面信息和ICESat-2高程點信息,28個時相的激光光子經過去噪后的有效光子數見表1,共10 4653個光子高程點;②根據與ICESat-2時相對應湖泊水面,剔除湖泊水面內的ICESat-2高程點,從而得到不同時相湖盆內無水區的高程點共82 050個;③對所有不同時相Sentinel-2提取的水面邊界線,根據水面等高特性[39],采用ICESat-2高程點內插出每個單獨水面線的高程,得到967 689個高程點;④得到湖盆內的高程點共計1 049 739個高程點。經過水面線高程點的加密,不僅湖盆內的高程點數增加了11.8倍,而且這些水面線還可以作為地形構建的高程約束線。經過加密后的湖盆高程點構建的湖盆DEM如圖5(b)所示。

圖5 構建臺特瑪湖DEMFig.5 Construction of Taitema lake DEM

在獲取精確的湖盆DEM后,結合2016—2022年臺特瑪湖的水面邊界,根據式(7)可以計算臺特瑪湖的水量時序,如圖6所示。由圖6可知,湖泊的水量變化與面積變化是高度一致的,湖泊在的面積在2017—2020年處于面積和水量較大的階段,這個階段湖泊在1—2月水量達到最大,其中2019年1月19號水量達到最大值,為0.894 km3。而在2016年和2021年以后,水量在3—4月達到最大,最大面積在0.189～0.234 km3之間,僅為水量高峰期的1/4。湖泊水量最小的時段均在7—8月,在2016—2022年幾乎沒有變化。

圖6 臺特瑪湖水量時序重建Fig.6 Time-series reconstruction of water volumes of the Taitema lake

3.2 精度評價與誤差分析

由于臺特瑪湖尚未有公開的湖盆地形或水量數據,很難直接驗證湖泊水量時序的精度。由于湖泊在不同水深下面積變化極大,且分裂成數量眾多的小水面,同時由圖5可知,臺特瑪湖的湖盆地形較為復雜, 這種情況下湖泊水量的計算很難用面積-水位關系方法計算,而湖泊水量的精度主要取決于湖盆地形的精度。為此,本文采用較為常用的檢查點法和等高線回放方法[40]評估湖盆DEM的精度。

3.2.1 檢查點法精度驗證

檢查點法是在DEM上設置一些檢查點作為高程真值,并從DEM上提取檢查點的內插高程,比較各個檢查點的誤差,然后算出中誤差用貝塞爾公式[41]計算高程中誤差,如果中誤差的離散度小,則表明DEM的精度高。這里選擇Sentinel-2影像2018年5月26日提取的湖泊水面線,分別從湖盆DEM上內插得到每條水面線的高程點信息,共得到185 493個檢查點。由于湖泊的每個小水面可以看作是等高的,分別計算每個水面線的平均高程作為真值,進而計算中誤差。由于中誤差只有在服從正態分布時才能反映數據的精度,本文先采用正態QQ圖對誤差分布進行了正態檢驗(圖7(a))。

圖7 檢查點的誤差分布Fig.7 Error distribution of the check points

由圖7可知,正態QQ圖圍繞對角線呈近似直線分布(R2=0.984),檢查點誤差服從正態分布統計特征[42]。計算這些水面檢查點高程的平均中誤差RMSE=0.103 m,從誤差分布圖來看,99.82%的誤差在3倍中誤差內。這表明同一水面線高程點的中誤差離散度很小,殘差的實際值與理論值越接近,模型的擬合效果越好。因此,從檢查點法定精度驗證表明DEM的誤差滿足極限誤差的要求。

3.2.2 等高線回放法精度驗證

等高線回放法[43]指的是對一條湖泊水面線,計算水面線點的平均高程,并根據該高程從湖盆DEM內插出等高線,并與實際的湖泊水面線的點一一對應,通過計算兩者的距離差來反映DEM的誤差。本文以2018年5月26日Sentinel-2影像提取的水面邊界線為例,首先計算每條水面線的平均高程值,并以此值內插出等高線,這里選取6個大小不一的水面及對應的等高線(圖8),分別計算水面與等高面的面積、對應點的最大和平均距離,來反映湖盆DEM與水面線的相似度(表2)。

表2 湖盆等高線與水面邊界的誤差對比

圖8 用于等高線回放法驗證的水面Fig.8 The water samples for accuracy assessment by the contour playback method

結果顯示圖8中的6個湖面線與生成的等高線空間吻合度很高。等高面與水面的面積誤差最大僅為1.564%,模擬的等高面的誤差已經低于湖泊水面提取的誤差,這表明湖盆DEM精度能很好地描述湖泊水面的精細變化。而等高線與水面邊界線的最大距離和最大平均間距分別為39.525及9.794 m,這表明每條模擬的等高線的平均間距在1個像素(10 m)內,最大誤差不超過4個像素。結合圖9可知,等高線與水面邊界線距離小于10、20及30 m的像素點比例分別大于71.5%、94.4%及99.8%,這表明利用DEM生成都等高線與水面線具有非常高的一致性,本文方法構建的DEM所提取的等高線能夠很好地吻合實際的湖面邊界線。

圖9 等高線與湖面邊界間距的累積百分比Fig.9 Cumulative frequency of distance between contour and water boundary

綜上所述,本文構建的湖盆地形中誤差僅為0.103 m,這表明湖盆的單點高程精度很高,經過水面線加密后的ICESat-2高程點能更好地體現。而且湖盆DEM提取的等高線能與實際的湖面線吻合度高,說明湖盆DEM已經能模擬淺平型湖盆地形的精細變化。

3.3 討 論

干旱區內陸尾閭湖泊經過長時間的沉積和演化,湖盆地形平坦,湖泊水深較淺,雖湖泊總水量較少,但常年保有最小生態水面是維持湖泊周邊生態系統健康的關鍵所在。本文就是為了精確計算缺資料淺平型內陸湖泊的水量變化,進而為湖泊生態補水提供科學依據,保障水資源的高效利用。由圖5的臺特瑪湖湖盆地形可以看出,臺特瑪湖湖盆地形整體上西高東低,整個湖盆地區東西橫跨約96 km,而海拔僅變化了約15 m,地形起伏非常細微。這類淺平型的湖泊的水面積變化通常較為明顯,而水量和水位變化相對非常微小。由圖6構建的臺特瑪湖水量時序可以看出,臺特瑪湖在2016—2019年間水面積和水量均呈動態增加的趨勢,年均變化率分別為144.64 km2/a和0.28 km3/a,在2019—2022年間水面積和水量均呈動態減少的趨勢,年均變化率分別為150.48 km2/a和0.29 km3/a。

傳統的基于遙感的“面積-高程-庫容”方法在應用于這類湖盆水量反演時,由于SRTM(shuttle radar topography mission)高程精度在2～10 m范圍內[44-45],難以獲得精確的水面-水位關系。利用ICESat-2高精度的激光光子高程點構建缺資料地區的湖盆地形,并利用水面等高的特性作為地形約束線加密湖盆內的高程點,湖盆地形精度可達0.1 m左右,能夠描述湖盆內地形的微小起伏,對于水面周期性、季節性變化顯著的湖泊,可大大提高水量時序反演的精度。然而,該方法的前提條件是,ICESat-2過境的激光光子束能夠探測到湖泊從水面最大到干涸的全過程,對于橫向距離小于3 km(ICESat-2軌道間距),以及水面變化不大的湖泊該方法難以適用。

4 結 論

針對干旱區淺平型尾閭湖水量難以估算的問題,本文基于湖泊水面邊界等高的特性,提出利用湖泊水面線高程加密ICESat-2地面高程點,構建的湖盆地形結合水面時序數據模擬湖盆水量變化時序,主要思路是通過密集的湖盆內無水區高程點構建湖盆地形,并結合時序面積信息重建湖泊水量時序,重點解決缺資料地區湖泊水文參數的信息反演。以臺特瑪湖為主要試驗區,重建2016—2022年臺特瑪湖湖泊水量變化,驗證了湖盆地形的反演精度。研究得出以下結論。

(1) 本文提出了基于改進湖盆反演技術的湖泊水量時序重建方法,能夠提取精度要求更高的淺平型湖泊地形構建,適宜于水位變化不敏感的湖泊水量估算。

(2) 改進湖盆反演技術構建的湖盆地形信息的高程中誤差僅為0.103 m,而且DEM提取的等高線與實際水面線高度吻合,這表明改進湖盆反演技術滿足淺平尾閭湖的高精度水量信息的反演。

(3) 基于改進湖盆反演技術的水量時序重建方法能夠適用于干旱區淺平型的內陸尾閭湖,也可以推廣到全球范圍內水面變化顯著且有衛星測高數據過境的湖泊,為無資料淺平型湖泊水量信息的時序重建提供一種思路。