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

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

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

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

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

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

1 研究方法

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

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

1.1 湖泊水面邊界線提取

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

1.2 PCL光子高程點提取

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

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

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

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鄰域范圍內(nèi)的光子點。

根據(jù)鄰域內(nèi)的所有光子計算中心光子點m的光子密度如下

(4)

(5)

(6)

式中,Pm為中心光子點m的鄰域點集;N為Pm中光子點數(shù)量;Wn是采用高斯核計算的橢圓鄰域中每個光子點的權(quán)重[33];δ是高斯核的標(biāo)準(zhǔn)差;Dm為中心光子點m的點密度,由鄰域內(nèi)所有點的高斯權(quán)重求和得到。

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

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

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

在計算PCL的地面光子高程信息后,利用與之時間對應(yīng)的湖泊圖層,提取湖盆范圍內(nèi)、湖面邊界外的高程點數(shù)據(jù),并將其高程基準(zhǔn)面轉(zhuǎn)換為大地水準(zhǔn)面,作為湖盆內(nèi)無水區(qū)的高程點。對不同時相的所有高程點進行提取,最后合并為覆蓋湖盆的PCL光子高程點。

1.3 湖泊水面邊界線高程內(nèi)插

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

圖3 內(nèi)插湖泊邊界高程Fig.3 Altitude interpolation of lake boundaries

1.4 湖盆地形構(gòu)建

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

1.5 湖泊水量時序重建

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

(7)

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

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

2.1 研究區(qū)概況

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

2.2 數(shù)據(jù)源選擇

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

表1 選擇的ICESat-2和Sentinel-2數(shù)據(jù)對應(yīng)列表

3 試驗與結(jié)果分析

3.1 試驗過程

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

圖5 構(gòu)建臺特瑪湖DEMFig.5 Construction of Taitema lake DEM

在獲取精確的湖盆DEM后,結(jié)合2016—2022年臺特瑪湖的水面邊界,根據(jù)式(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 精度評價與誤差分析

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

3.2.1 檢查點法精度驗證

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

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

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

3.2.2 等高線回放法精度驗證

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

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

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

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

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

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

3.3 討 論

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

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

4 結(jié) 論

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

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

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

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