歐海沨,字城岱,滕興發(fā),關(guān)舒丹
(西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院,四川 成都 611756)
冰川受人類活動影響較小,受氣候變化影響的同時,也對其作出快速的響應(yīng)[1-2]。在全球氣候變暖的影響下,我國大部分冰川呈現(xiàn)出退縮減弱趨勢[3],導(dǎo)致由冰川變化引發(fā)的區(qū)域淹沒、冰湖潰決、滑坡和碎屑流等自然災(zāi)害的規(guī)模增大、數(shù)量增加。提取并分析冰川表面運動信息,不僅能揭示全球氣候變暖背景下冰川的動態(tài)變化信息,同時也能為冰川變化誘導(dǎo)的災(zāi)害預(yù)警提供重要依據(jù)[4-6]。由于冰川多分布在海拔較高、地勢嚴峻、條件惡劣的地區(qū),傳統(tǒng)的監(jiān)測技術(shù)受到諸多限制。而空間遙感技術(shù)體現(xiàn)出更多的優(yōu)勢,包括覆蓋范圍大、數(shù)據(jù)源廣、監(jiān)測周期長、無須現(xiàn)場勘測等。目前,在冰川表面流速的監(jiān)測中,普遍采用的遙感手段包括3種:合成孔徑雷達(SAR)影像的差分干涉[7-8]、SAR像素偏移追蹤技術(shù)(POT)[9-10]及光學(xué)影像互相關(guān)技術(shù)(OICC)[11-12]。
其中,光學(xué)影像具有更加豐富的紋理信息,更有利于監(jiān)測地表的幾何變化,同時,接近垂直的入射角也可避免陡坡之間的相互遮擋。雖然光學(xué)影像質(zhì)量易受云層的限制,但隨著高性能光學(xué)遙感衛(wèi)星數(shù)量的不斷增加,光學(xué)影像的存檔數(shù)量也在不斷積累,為及時、準確獲取冰川動態(tài)變化信息提供了可能。本文以加拉白壘峰北側(cè)的崗普冰川為研究對象,位于雅魯藏布江大拐彎地帶,該地區(qū)冰川多屬于海洋型冰川,運動速度較快;借鑒短基線集(SBAS)思想,利用OICC-SBAS技術(shù)提取冰川表面運動狀態(tài),并分析其時空變化特征。
加拉白壘峰位于念唐古拉山脈,屬于高原溫帶半濕潤季風氣候,具有雨量豐沛、干濕季節(jié)分明、晝夜溫差大等特點,受降雨和溫度的影響,該區(qū)域內(nèi)冰崩和雪崩等災(zāi)害十分頻繁。研究對象崗普冰川為加拉白壘峰周邊規(guī)模較大的冰川之一,地勢南高北低,東側(cè)發(fā)育著一條大型的支流(如圖1所示),區(qū)域內(nèi)最高海拔為4169 m,最低海拔為3438 m。
圖1 研究區(qū)地理位置及冰川地貌
Sentinel-2包含A、B兩顆衛(wèi)星,已收集了約7年的存檔數(shù)據(jù),且具有重訪周期短、中等分辨率和免費開放等特點。Sentinel-2傳感器共提供13個光譜波段,空間分辨率位于10~60 m之間,本文采用空間分辨率最高的近紅外波段。從歐洲航天局(European Space Agency,ESA)網(wǎng)站獲取了2015年12月6日至2022年3月9日的12景Sentinel-2遙感影像,詳細參數(shù)信息見表1。
表1 Sentinel-2影像參數(shù)
為獲取崗普冰川的時序表面運動速度及時空變化特征,本文采用OICC-SBAS時間序列反演算法,如圖2所示。該算法共包含5個模塊:影像配對策略、數(shù)據(jù)預(yù)處理、COSI-Corr批處理、誤差后處理及時序反演。
圖2 數(shù)據(jù)處理流程
不同影像之間存在著太陽光照條件、軌道位置、采集時間等因素的差異,這些因素與圖像相關(guān)性測量值之間存在著一定的關(guān)系,被稱為基線變量[13],包括太陽方位角差值(SAD)、太陽高度角差值(SED)、時間基線(TB)和空間基線(SB)。通過設(shè)置基線變量的閾值,可優(yōu)化影像對的選擇,具體見表2。在后續(xù)處理中手動去除質(zhì)量較差的影像對,最終獲取14組影像對。由于獲取的Sentinel-2數(shù)據(jù)為L1C級數(shù)據(jù),即經(jīng)過了粗略幾何校正和正射校正[13-14],可直接用于影像互相關(guān)處理。預(yù)處理只需提取覆蓋研究區(qū)域的范圍,即裁剪出包含相同地理覆蓋的子圖像,并按時間順序排列,形成圖像堆棧。
表2 基線閾值與批處理相關(guān)參數(shù)
以COSI-Corr軟件為數(shù)據(jù)處理平臺,經(jīng)多次試驗,具體參數(shù)設(shè)置為:初始搜索窗口和最終搜索窗口均設(shè)置為64×64,步長設(shè)置為1,其他參數(shù)(穩(wěn)健性迭代次數(shù)和掩膜閾值)均采用系統(tǒng)默認值,確保在降低測量不確定度的同時具有較高的計算效率。
初始相關(guān)結(jié)果中仍然存在大量的噪聲和偽影,包括失相關(guān)噪聲誤差、軌道誤差、條帶誤差、衛(wèi)星姿態(tài)角誤差,需要進一步后處理。針對異常值和失相關(guān)噪聲,將信噪比值低于0.96的像素和超過TB閾值估計的最大預(yù)期位移的像素進行掩膜。對于軌道誤差,采用一次多項式仿射變換模型解決;對于條帶誤差,通過均值相減法去除;對于衛(wèi)星姿態(tài)角誤差,采用改進的均值相減法[15]去除;最后再對相關(guān)圖進行非局部均值濾波,獲得進行了誤差處理的東西向形變圖和南北向形變圖。
基于互相關(guān)技術(shù)的位移時間序列反演借鑒了InSAR時序反演中的小基線時序分析技術(shù)(SBAS)[15]。假設(shè)有N+1幅光學(xué)影像,獲取時間序列為[t0,t1,…,tn],M為互相關(guān)影像對。以影像對在每個連續(xù)時間間隔上的平均速率v作為未知參數(shù),構(gòu)建觀測方程為
B·v=δd
(1)
式中,δd為影像對的位移序列;B為由相鄰影像的時間間隔組成系數(shù)矩陣。當系數(shù)矩陣B滿秩時,可采用最小二乘求解代求參數(shù)v,即
v=(BTB)-1BTδd
(2)
但當系數(shù)矩陣B秩虧時,需要采用SVD分解法求解參數(shù)v。為了更好地考慮去相關(guān)隨時間的變化,根據(jù)主影像和從影像之間的時間間隔,加入權(quán)重矩陣P,P=diag(w1,w2,…,wM),其中
(3)
Rl為原始測量值和反演位移之間差值,即Rl=δd-Bv。高殘差代表高偏差或低質(zhì)量的相關(guān)值,這些值會直接影響反演的精度。為了避免該問題,利用殘差Rl對反演進行第2次迭代,即
(4)
式中,R0為一個常數(shù),近似于“成功”相關(guān)性的噪聲標準差。因此,權(quán)重Wil限制了非常高的殘差異常值對反演結(jié)果的影響,但不會人為地高估殘差接近于零的相關(guān)值。獲取反演結(jié)果后,可以通過掩膜和濾波對結(jié)果進一步改進。首先,根據(jù)文獻[16]中的方法,過濾方向標準差較大的像素值,并去除運動方向與實際斜坡方向大于90°的值;然后,掩膜后對導(dǎo)出的位移時間序列進行方向濾波[17]。
為驗證SAD與SED對其相關(guān)性結(jié)果的影響程度,在研究區(qū)域內(nèi)選取3個穩(wěn)定區(qū)域(如圖1(b)中的Roi1、Roi2、Roi3),對3個區(qū)域東西向和南北向的形變均值進行精度評定,最終獲得的相互關(guān)系圖如圖3所示。通過統(tǒng)計發(fā)現(xiàn),當SAD與SED較大時,其均值較高;當SAD與SED較小時,其均值較低。這充分說明在適當?shù)臅r間基線下,選擇較小的SAD與SED,能夠有效提高結(jié)果精度。
圖3 太陽高度角差值和太陽方位角差值與穩(wěn)定區(qū)域形變均值的關(guān)系
由于缺乏冰川的實測數(shù)據(jù),在精度評定時,通常采用穩(wěn)定區(qū)域的均值和標準差等先驗信息進行輔助。理論上穩(wěn)定區(qū)域的速率為零,因此提取3個穩(wěn)定區(qū)域的殘余運動,并繪制相應(yīng)的殘余運動頻率分布圖(如圖4所示)。其中,Roi1、Roi2靠近冰川,地形起伏較大,且在一定程度上受積雪影響,殘余運動偏大;Roi3離冰川較遠,殘余運動明顯減小。誤差統(tǒng)計結(jié)果見表3。結(jié)果顯示,誤差值相對于冰川表面最大流速范圍為0.04%~0.7%,穩(wěn)定區(qū)域平均運動速度遠小于冰川實際運動速度,保證了冰川時序流速結(jié)果的可靠性與準確性。
表3 穩(wěn)定區(qū)域運動速度統(tǒng)計結(jié)果(m/a)
圖4 穩(wěn)定區(qū)域殘余運動頻率分布
本文通過2015—2022年獲取的12幅Sentinel-2光學(xué)影像,基于OICC-SBAS反演算法提取了崗普冰川的年平均表面流速(如圖5所示)。由于兩側(cè)山體擠壓,冰川中間的表面流速大于兩側(cè),符合冰川運動規(guī)律。同時,利用東西向和南北向偏移量獲取的偏移方向,與冰川實際運動方向相符(如圖6(a)所示)。冰川主流以南北向形變?yōu)橹?支流以東西向形變?yōu)橹?這主要是由坡度、坡向等地形因子聯(lián)合導(dǎo)致的。為了進一步分析冰川表面流速的空間分布特征,在崗普冰川上提取了3條剖面線(剖面線位置如圖5(b)中的AA′、BB′、CC′),并聯(lián)合高程、坡度信息進行分析,結(jié)果如圖6所示。通過提取冰川剖線與地形因子的關(guān)系,得到如下結(jié)論。
圖5 冰川年均表面流速
圖6 崗普冰川剖線表面流速與地形因子關(guān)系
(1)剖面線AA′:剖面線AA′為冰川主流中心線,以南北向形變?yōu)橹?。沿該剖面線的高程和坡度變化較平緩,流速波動較小,但也可以看出3個波動現(xiàn)象,分別位于剖面線的0.5、3.2、5.6 km處(分別為圖6(a)中的T1、T2、T3)。距T1處的加速主要是由于上游支流在此匯集,從而導(dǎo)致冰川加速;而在T2處,雖有東側(cè)支流的匯入,但冰川流速減弱,主要是由于該支流在匯入主冰川之前其流速已經(jīng)減小至1.24 m/a,因此對主冰川流速未表現(xiàn)出較大的推動作用;T3位于冰舌末端,由于冰川物質(zhì)積累和地形變化,也有一個加速現(xiàn)象。
(2)剖面線BB′:剖面線BB′位于冰川的東側(cè)支流,以東西向流速為主。沿剖面線BB′的地形變化較明顯,該支流流速隨高程的減小而逐漸減小;在剖線末端,由于匯入主流,導(dǎo)致南北向流速有一個明顯的增速現(xiàn)象。
(3)剖面線CC′:該剖面線橫跨上游東西兩條支流,地形為中間低、兩頭高。該支流流速在1.0、3.9 km處(分別為圖6(c)中的T1、T2)隨著冰磧物的積累和地形變換,其東西向流速有一段快速加速過程;后因逐漸匯入主流,開始以南北向流速為主導(dǎo)。冰川西側(cè)的上游支流為表面流速高值的主要聚集區(qū)域,該區(qū)域為冰川的物質(zhì)積累區(qū)域,形變量級較大且形變特征復(fù)雜,年最高流速可達82.5 m/a。
綜上所述,崗普冰川主體流速受支流影響,但整體無劇烈變化,維持在40~50 m/a之間。東側(cè)支流隨著高程減小,其冰川流速也逐漸減小,位于20~40 m/a之間。崗普冰川上游為整條冰川流速高值的主要聚集區(qū)域,形變量大且形變特征較為復(fù)雜。
通過OICC-SBAS技術(shù)處理得到的2015—2022年的冰川時序形變結(jié)果如圖7所示??梢悦黠@看出,2017年2月—2018年3月、2020年11月—2021年12月,這兩時段內(nèi)的累計形變明顯大于其他時期。一方面是由于這兩時段的時間間隔較長,另一方面則是受氣候因素的影響。
圖7 時序形變處理結(jié)果
為了研究崗普冰川的變化規(guī)律,通過美國國家海洋和大氣管理局(NOAA)網(wǎng)站獲取距離研究區(qū)域最近的林芝站點降水量的逐日統(tǒng)計結(jié)果。在崗普冰川上選取3個區(qū)域(分別為圖7中的A、B、C),分析冰川形變與氣候因素的關(guān)系(如圖8所示)??梢钥闯?研究時段內(nèi)年降水量具有明顯差異,2017年2月—2018年3月、2020年11月—2021年12月均具有較高的累計降水量。其中,2016年3月—2017年2月具有最高的累計降水量。該時段時間間隔較短,導(dǎo)致累計形變不是最大,但該時段形變?nèi)匀伙@著高于其他累計降水量較小的時段??傮w而言,冰川的形變量與降水變化較為一致,在累計降水量高且時間間隔較長的時段內(nèi),冰川累計形變較高,而在累計降水量較低的時段內(nèi),冰川的形變量有明顯減少。
圖8 崗普冰川累計形變與氣溫降水關(guān)系
本文利用2015—2022年獲取的12景Sentinel-2影像,采用OICC-SBAS反演算法,提取崗普冰川表面時序運動速度分布信息,并進行定性、定量分析。空間分布特征上,崗普冰川位移受地形和氣候因子的共同作用,流速高值主要集中在上游的物質(zhì)積累區(qū),該區(qū)形變特征較為復(fù)雜,最高流速超過80 m/a。冰川主體運動速度主要為40~50 m/a,受地形和支流匯集影響,表現(xiàn)出輕微波動,但整體較為穩(wěn)定。東側(cè)支流作為崗普冰川水資源的重要補給,其運動速度隨著高程減小而緩慢減小,主要維持在20~40 m/a。時序變化特征上,研究時段內(nèi)冰川流速無較大變化,年際之間雖然有輕微的差距,但整體呈現(xiàn)出較高的相似性。即2015—2022年,崗普冰川都維持在一個穩(wěn)定流動的狀態(tài)。且崗普冰川的流速變化與降水變化較為一致,能較好地反映出氣候變化的特征。
同時,本文利用東西向和南北向偏移量獲取的最終偏移方向與冰川實際流動方向相符,偏移量誤差僅占冰川最大流速的0.04%~0.7%,穩(wěn)定區(qū)域的殘余運動遠小于冰川的實際運動速度,有力地證實了冰川時序表面流速結(jié)果具有較高的可靠性和檢測精度。