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

歐海沨,字城岱,滕興發(fā),關(guān)舒丹

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

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

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

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

1.1 研究區(qū)概況

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

圖1 研究區(qū)地理位置及冰川地貌

1.2 研究數(shù)據(jù)

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

表1 Sentinel-2影像參數(shù)

2 研究方法與數(shù)據(jù)處理

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

圖2 數(shù)據(jù)處理流程

2.1 影像配對策略與預(yù)處理

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

表2 基線閾值與批處理相關(guān)參數(shù)

2.2 偏移量估計(jì)與誤差后處理

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

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

2.3 時(shí)序偏移量計(jì)算

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

B·v=δd

(1)

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

v=(BTB)-1BTδd

(2)

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

(3)

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

(4)

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

3 結(jié)果與分析

3.1 冰川運(yùn)動結(jié)果可靠性分析

為驗(yàn)證SAD與SED對其相關(guān)性結(jié)果的影響程度,在研究區(qū)域內(nèi)選取3個(gè)穩(wěn)定區(qū)域(如圖1(b)中的Roi1、Roi2、Roi3),對3個(gè)區(qū)域東西向和南北向的形變均值進(jìn)行精度評定,最終獲得的相互關(guān)系圖如圖3所示。通過統(tǒng)計(jì)發(fā)現(xiàn),當(dāng)SAD與SED較大時(shí),其均值較高;當(dāng)SAD與SED較小時(shí),其均值較低。這充分說明在適當(dāng)?shù)臅r(shí)間基線下,選擇較小的SAD與SED,能夠有效提高結(jié)果精度。

圖3 太陽高度角差值和太陽方位角差值與穩(wěn)定區(qū)域形變均值的關(guān)系

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

表3 穩(wěn)定區(qū)域運(yùn)動速度統(tǒng)計(jì)結(jié)果(m/a)

圖4 穩(wěn)定區(qū)域殘余運(yùn)動頻率分布

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

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

圖5 冰川年均表面流速

圖6 崗普冰川剖線表面流速與地形因子關(guān)系

(1)剖面線AA′:剖面線AA′為冰川主流中心線,以南北向形變?yōu)橹?。沿該剖面線的高程和坡度變化較平緩,流速波動較小,但也可以看出3個(gè)波動現(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ì)積累和地形變化,也有一個(gè)加速現(xiàn)象。

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

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

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

3.3 冰川表面流速時(shí)序變化特征

通過OICC-SBAS技術(shù)處理得到的2015—2022年的冰川時(shí)序形變結(jié)果如圖7所示?？梢悦黠@看出,2017年2月—2018年3月、2020年11月—2021年12月,這兩時(shí)段內(nèi)的累計(jì)形變明顯大于其他時(shí)期。一方面是由于這兩時(shí)段的時(shí)間間隔較長,另一方面則是受氣候因素的影響。

圖7 時(shí)序形變處理結(jié)果

為了研究崗普冰川的變化規(guī)律,通過美國國家海洋和大氣管理局(NOAA)網(wǎng)站獲取距離研究區(qū)域最近的林芝站點(diǎn)降水量的逐日統(tǒng)計(jì)結(jié)果。在崗普冰川上選取3個(gè)區(qū)域(分別為圖7中的A、B、C),分析冰川形變與氣候因素的關(guān)系(如圖8所示)?？梢钥闯?研究時(shí)段內(nèi)年降水量具有明顯差異,2017年2月—2018年3月、2020年11月—2021年12月均具有較高的累計(jì)降水量。其中,2016年3月—2017年2月具有最高的累計(jì)降水量。該時(shí)段時(shí)間間隔較短,導(dǎo)致累計(jì)形變不是最大,但該時(shí)段形變?nèi)匀伙@著高于其他累計(jì)降水量較小的時(shí)段?？傮w而言,冰川的形變量與降水變化較為一致,在累計(jì)降水量高且時(shí)間間隔較長的時(shí)段內(nèi),冰川累計(jì)形變較高,而在累計(jì)降水量較低的時(shí)段內(nèi),冰川的形變量有明顯減少。

圖8 崗普冰川累計(jì)形變與氣溫降水關(guān)系

4 結(jié) 語

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

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