利用OICC-SBAS技術提取崗普冰川表面流速時空變化特征

歐海沨,字城岱,滕興發,關舒丹

(西南交通大學地球科學與環境工程學院,四川 成都 611756)

冰川受人類活動影響較小,受氣候變化影響的同時,也對其作出快速的響應[1-2]。在全球氣候變暖的影響下,我國大部分冰川呈現出退縮減弱趨勢[3],導致由冰川變化引發的區域淹沒、冰湖潰決、滑坡和碎屑流等自然災害的規模增大、數量增加。提取并分析冰川表面運動信息,不僅能揭示全球氣候變暖背景下冰川的動態變化信息,同時也能為冰川變化誘導的災害預警提供重要依據[4-6]。由于冰川多分布在海拔較高、地勢嚴峻、條件惡劣的地區,傳統的監測技術受到諸多限制。而空間遙感技術體現出更多的優勢,包括覆蓋范圍大、數據源廣、監測周期長、無須現場勘測等。目前,在冰川表面流速的監測中,普遍采用的遙感手段包括3種:合成孔徑雷達(SAR)影像的差分干涉[7-8]、SAR像素偏移追蹤技術(POT)[9-10]及光學影像互相關技術(OICC)[11-12]。

其中,光學影像具有更加豐富的紋理信息,更有利于監測地表的幾何變化,同時,接近垂直的入射角也可避免陡坡之間的相互遮擋。雖然光學影像質量易受云層的限制,但隨著高性能光學遙感衛星數量的不斷增加,光學影像的存檔數量也在不斷積累,為及時、準確獲取冰川動態變化信息提供了可能。本文以加拉白壘峰北側的崗普冰川為研究對象,位于雅魯藏布江大拐彎地帶,該地區冰川多屬于海洋型冰川,運動速度較快;借鑒短基線集(SBAS)思想,利用OICC-SBAS技術提取冰川表面運動狀態,并分析其時空變化特征。

1 研究區與數據

1.1 研究區概況

加拉白壘峰位于念唐古拉山脈,屬于高原溫帶半濕潤季風氣候,具有雨量豐沛、干濕季節分明、晝夜溫差大等特點,受降雨和溫度的影響,該區域內冰崩和雪崩等災害十分頻繁。研究對象崗普冰川為加拉白壘峰周邊規模較大的冰川之一,地勢南高北低,東側發育著一條大型的支流(如圖1所示),區域內最高海拔為4169 m,最低海拔為3438 m。

圖1 研究區地理位置及冰川地貌

1.2 研究數據

Sentinel-2包含A、B兩顆衛星,已收集了約7年的存檔數據,且具有重訪周期短、中等分辨率和免費開放等特點。Sentinel-2傳感器共提供13個光譜波段,空間分辨率位于10～60 m之間,本文采用空間分辨率最高的近紅外波段。從歐洲航天局(European Space Agency,ESA)網站獲取了2015年12月6日至2022年3月9日的12景Sentinel-2遙感影像,詳細參數信息見表1。

表1 Sentinel-2影像參數

2 研究方法與數據處理

為獲取崗普冰川的時序表面運動速度及時空變化特征,本文采用OICC-SBAS時間序列反演算法,如圖2所示。該算法共包含5個模塊:影像配對策略、數據預處理、COSI-Corr批處理、誤差后處理及時序反演。

圖2 數據處理流程

2.1 影像配對策略與預處理

不同影像之間存在著太陽光照條件、軌道位置、采集時間等因素的差異,這些因素與圖像相關性測量值之間存在著一定的關系,被稱為基線變量[13],包括太陽方位角差值(SAD)、太陽高度角差值(SED)、時間基線(TB)和空間基線(SB)。通過設置基線變量的閾值,可優化影像對的選擇,具體見表2。在后續處理中手動去除質量較差的影像對,最終獲取14組影像對。由于獲取的Sentinel-2數據為L1C級數據,即經過了粗略幾何校正和正射校正[13-14],可直接用于影像互相關處理。預處理只需提取覆蓋研究區域的范圍,即裁剪出包含相同地理覆蓋的子圖像,并按時間順序排列,形成圖像堆棧。

表2 基線閾值與批處理相關參數

2.2 偏移量估計與誤差后處理

以COSI-Corr軟件為數據處理平臺,經多次試驗,具體參數設置為:初始搜索窗口和最終搜索窗口均設置為64×64,步長設置為1,其他參數(穩健性迭代次數和掩膜閾值)均采用系統默認值,確保在降低測量不確定度的同時具有較高的計算效率。

初始相關結果中仍然存在大量的噪聲和偽影,包括失相關噪聲誤差、軌道誤差、條帶誤差、衛星姿態角誤差,需要進一步后處理。針對異常值和失相關噪聲,將信噪比值低于0.96的像素和超過TB閾值估計的最大預期位移的像素進行掩膜。對于軌道誤差,采用一次多項式仿射變換模型解決;對于條帶誤差,通過均值相減法去除;對于衛星姿態角誤差,采用改進的均值相減法[15]去除;最后再對相關圖進行非局部均值濾波,獲得進行了誤差處理的東西向形變圖和南北向形變圖。

2.3 時序偏移量計算

基于互相關技術的位移時間序列反演借鑒了InSAR時序反演中的小基線時序分析技術(SBAS)[15]。假設有N+1幅光學影像,獲取時間序列為[t0,t1,…,tn],M為互相關影像對。以影像對在每個連續時間間隔上的平均速率v作為未知參數,構建觀測方程為

B·v=δd

(1)

式中,δd為影像對的位移序列;B為由相鄰影像的時間間隔組成系數矩陣。當系數矩陣B滿秩時,可采用最小二乘求解代求參數v,即

v=(BTB)-1BTδd

(2)

但當系數矩陣B秩虧時,需要采用SVD分解法求解參數v。為了更好地考慮去相關隨時間的變化,根據主影像和從影像之間的時間間隔,加入權重矩陣P,P=diag(w1,w2,…,wM),其中

(3)

Rl為原始測量值和反演位移之間差值,即Rl=δd-Bv。高殘差代表高偏差或低質量的相關值,這些值會直接影響反演的精度。為了避免該問題,利用殘差Rl對反演進行第2次迭代,即

(4)

式中,R0為一個常數,近似于“成功”相關性的噪聲標準差。因此,權重Wil限制了非常高的殘差異常值對反演結果的影響,但不會人為地高估殘差接近于零的相關值。獲取反演結果后,可以通過掩膜和濾波對結果進一步改進。首先,根據文獻[16]中的方法,過濾方向標準差較大的像素值,并去除運動方向與實際斜坡方向大于90°的值;然后,掩膜后對導出的位移時間序列進行方向濾波[17]。

3 結果與分析

3.1 冰川運動結果可靠性分析

為驗證SAD與SED對其相關性結果的影響程度,在研究區域內選取3個穩定區域(如圖1(b)中的Roi1、Roi2、Roi3),對3個區域東西向和南北向的形變均值進行精度評定,最終獲得的相互關系圖如圖3所示。通過統計發現,當SAD與SED較大時,其均值較高;當SAD與SED較小時,其均值較低。這充分說明在適當的時間基線下,選擇較小的SAD與SED,能夠有效提高結果精度。

圖3 太陽高度角差值和太陽方位角差值與穩定區域形變均值的關系

由于缺乏冰川的實測數據,在精度評定時,通常采用穩定區域的均值和標準差等先驗信息進行輔助。理論上穩定區域的速率為零,因此提取3個穩定區域的殘余運動,并繪制相應的殘余運動頻率分布圖(如圖4所示)。其中,Roi1、Roi2靠近冰川,地形起伏較大,且在一定程度上受積雪影響,殘余運動偏大;Roi3離冰川較遠,殘余運動明顯減小。誤差統計結果見表3。結果顯示,誤差值相對于冰川表面最大流速范圍為0.04%～0.7%,穩定區域平均運動速度遠小于冰川實際運動速度,保證了冰川時序流速結果的可靠性與準確性。

表3 穩定區域運動速度統計結果(m/a)

圖4 穩定區域殘余運動頻率分布

3.2 冰川表面流速空間分布特征

本文通過2015—2022年獲取的12幅Sentinel-2光學影像,基于OICC-SBAS反演算法提取了崗普冰川的年平均表面流速(如圖5所示)。由于兩側山體擠壓,冰川中間的表面流速大于兩側,符合冰川運動規律。同時,利用東西向和南北向偏移量獲取的偏移方向,與冰川實際運動方向相符(如圖6(a)所示)。冰川主流以南北向形變為主,支流以東西向形變為主,這主要是由坡度、坡向等地形因子聯合導致的。為了進一步分析冰川表面流速的空間分布特征,在崗普冰川上提取了3條剖面線(剖面線位置如圖5(b)中的AA′、BB′、CC′),并聯合高程、坡度信息進行分析,結果如圖6所示。通過提取冰川剖線與地形因子的關系,得到如下結論。

圖5 冰川年均表面流速

圖6 崗普冰川剖線表面流速與地形因子關系

(1)剖面線AA′:剖面線AA′為冰川主流中心線,以南北向形變為主。沿該剖面線的高程和坡度變化較平緩,流速波動較小,但也可以看出3個波動現象,分別位于剖面線的0.5、3.2、5.6 km處(分別為圖6(a)中的T1、T2、T3)。距T1處的加速主要是由于上游支流在此匯集,從而導致冰川加速;而在T2處,雖有東側支流的匯入,但冰川流速減弱,主要是由于該支流在匯入主冰川之前其流速已經減小至1.24 m/a,因此對主冰川流速未表現出較大的推動作用;T3位于冰舌末端,由于冰川物質積累和地形變化,也有一個加速現象。

(2)剖面線BB′:剖面線BB′位于冰川的東側支流,以東西向流速為主。沿剖面線BB′的地形變化較明顯,該支流流速隨高程的減小而逐漸減小;在剖線末端,由于匯入主流,導致南北向流速有一個明顯的增速現象。

(3)剖面線CC′:該剖面線橫跨上游東西兩條支流,地形為中間低、兩頭高。該支流流速在1.0、3.9 km處(分別為圖6(c)中的T1、T2)隨著冰磧物的積累和地形變換,其東西向流速有一段快速加速過程;后因逐漸匯入主流,開始以南北向流速為主導。冰川西側的上游支流為表面流速高值的主要聚集區域,該區域為冰川的物質積累區域,形變量級較大且形變特征復雜,年最高流速可達82.5 m/a。

綜上所述,崗普冰川主體流速受支流影響,但整體無劇烈變化,維持在40～50 m/a之間。東側支流隨著高程減小,其冰川流速也逐漸減小,位于20～40 m/a之間。崗普冰川上游為整條冰川流速高值的主要聚集區域,形變量大且形變特征較為復雜。

3.3 冰川表面流速時序變化特征

通過OICC-SBAS技術處理得到的2015—2022年的冰川時序形變結果如圖7所示。可以明顯看出,2017年2月—2018年3月、2020年11月—2021年12月,這兩時段內的累計形變明顯大于其他時期。一方面是由于這兩時段的時間間隔較長,另一方面則是受氣候因素的影響。

圖7 時序形變處理結果

為了研究崗普冰川的變化規律,通過美國國家海洋和大氣管理局(NOAA)網站獲取距離研究區域最近的林芝站點降水量的逐日統計結果。在崗普冰川上選取3個區域(分別為圖7中的A、B、C),分析冰川形變與氣候因素的關系(如圖8所示)。可以看出,研究時段內年降水量具有明顯差異,2017年2月—2018年3月、2020年11月—2021年12月均具有較高的累計降水量。其中,2016年3月—2017年2月具有最高的累計降水量。該時段時間間隔較短,導致累計形變不是最大,但該時段形變仍然顯著高于其他累計降水量較小的時段。總體而言,冰川的形變量與降水變化較為一致,在累計降水量高且時間間隔較長的時段內,冰川累計形變較高,而在累計降水量較低的時段內,冰川的形變量有明顯減少。

圖8 崗普冰川累計形變與氣溫降水關系

4 結 語

本文利用2015—2022年獲取的12景Sentinel-2影像,采用OICC-SBAS反演算法,提取崗普冰川表面時序運動速度分布信息,并進行定性、定量分析。空間分布特征上,崗普冰川位移受地形和氣候因子的共同作用,流速高值主要集中在上游的物質積累區,該區形變特征較為復雜,最高流速超過80 m/a。冰川主體運動速度主要為40～50 m/a,受地形和支流匯集影響,表現出輕微波動,但整體較為穩定。東側支流作為崗普冰川水資源的重要補給,其運動速度隨著高程減小而緩慢減小,主要維持在20～40 m/a。時序變化特征上,研究時段內冰川流速無較大變化,年際之間雖然有輕微的差距,但整體呈現出較高的相似性。即2015—2022年,崗普冰川都維持在一個穩定流動的狀態。且崗普冰川的流速變化與降水變化較為一致,能較好地反映出氣候變化的特征。

同時,本文利用東西向和南北向偏移量獲取的最終偏移方向與冰川實際流動方向相符,偏移量誤差僅占冰川最大流速的0.04%～0.7%,穩定區域的殘余運動遠小于冰川的實際運動速度,有力地證實了冰川時序表面流速結果具有較高的可靠性和檢測精度。