王紀(jì)碩, 張齊民, 孫一丹, 張豪磊, 張子彥, 閆世勇

(1.中國礦業(yè)大學(xué) 環(huán)境與測繪學(xué)院,江蘇 徐州 221116; 2.自然資源部國土環(huán)境與災(zāi)害監(jiān)測重點(diǎn)實(shí)驗(yàn)室,江蘇 徐州 221116; 3.山東省國土測繪院 基礎(chǔ)測繪中心,山東 濟(jì)南 250102)

在當(dāng)前全球變暖的背景下[1],以青藏高原為核心的“亞洲水塔”正在發(fā)生以失衡為特征的劇烈變化,青藏高原山地冰川加速退縮,對(duì)人類社會(huì)文明和生態(tài)安全的影響不可低估[2]。青藏高原喀喇昆侖山是全球山岳冰川最發(fā)達(dá)的高大山脈,也是中緯度地區(qū)現(xiàn)代冰川最發(fā)育的地帶,而冰川作為當(dāng)?shù)刂匾牡Y源庫,其厚度與儲(chǔ)量的變化對(duì)周邊區(qū)域的生態(tài)平衡起著至關(guān)重要的作用[3]。

獲取冰川厚度的方法主要有實(shí)地測量和間接估計(jì)2類。采用人工鉆孔法進(jìn)行實(shí)地測量,會(huì)受到地理環(huán)境限制,無法進(jìn)行大范圍作業(yè)。根據(jù)厚度-面積比例關(guān)系式可間接估計(jì)冰川的平均厚度,但由于不同冰川各自的屬性不同,使得該方法針對(duì)區(qū)域性冰川有效,對(duì)于單一冰川并不適用。20世紀(jì)80年代初期,我國首次采用探地雷達(dá)(ground penetrating radar,GPR)對(duì)天山烏魯木齊河源1號(hào)冰川進(jìn)行冰厚測量并取得一些成果,隨著GPR技術(shù)的不斷發(fā)展,該方法被廣泛應(yīng)用在我國青藏高原山地冰川的冰厚提取與冰儲(chǔ)量估算研究中,但GPR探測冰厚的準(zhǔn)確性受到電磁波穿透冰體能力和雷達(dá)系統(tǒng)圖像處理能力的制約,并且冰層的細(xì)微變化會(huì)分散GPR的信號(hào),從而降低圖像清晰度,當(dāng)前GPR技術(shù)多應(yīng)用于小型冰川厚度提取,在青藏高原大規(guī)模冰川厚度測量中的應(yīng)用很少[4]。衛(wèi)星遙感技術(shù)的發(fā)展使得大規(guī)模監(jiān)測冰川成為可能,在山地冰川應(yīng)用較多的是美國航天飛機(jī)雷達(dá)地形測繪使命(Shuttle Radar Topography Mission,SRTM) 數(shù)字高程模型(digital elevation model,DEM)數(shù)據(jù)和Landsat衛(wèi)星光學(xué)影像。文獻(xiàn)[5]根據(jù)淺冰近似(shallow ice approximation,SIA)理論,采用SRTM DEM數(shù)據(jù)進(jìn)行冰川厚度反演。SIA理論是基于冰川中流線平均基底剪切應(yīng)力經(jīng)驗(yàn)公式和冰川中流線平均坡度進(jìn)行冰川厚度反演,然而該方法只能求得冰川的平均冰厚和總體積,并且由于DEM數(shù)據(jù)獲取時(shí)間具有單一性,難以獲取冰川坡度的年際變化,同時(shí)冰川表面運(yùn)動(dòng)、溫度、降水等因素也都嚴(yán)重影響冰川厚度的變化,因而不能為冰川厚度的持續(xù)性變化提供參考。文獻(xiàn)[6]采用Landsat衛(wèi)星影像獲取青藏高原東南部冰川的表面狀態(tài),根據(jù)質(zhì)量平衡分布數(shù)據(jù)、冰川表面速度和冰流原理估算出冰川冰厚的空間分布;文獻(xiàn)[7]采用Landsat 7、Sentinel-2衛(wèi)星影像對(duì)喜馬拉雅山脈冰川的流速進(jìn)行反演,探索冰川動(dòng)態(tài)變化與質(zhì)量守恒的相互作用機(jī)制。但光學(xué)影像在青藏高原冰川監(jiān)測中的應(yīng)用極大地受限于惡劣的天氣狀況[8],而雷達(dá)影像具有穿透力強(qiáng)、覆蓋范圍廣、不受云雨天氣影響等優(yōu)點(diǎn),利用不同時(shí)期的雷達(dá)影像,基于冰川流速對(duì)冰川厚度進(jìn)行估算,能夠簡便高效地通過數(shù)值模擬的方法實(shí)現(xiàn)冰川厚度的反演。本文基于Sentinel-1A衛(wèi)星干涉寬幅(Interferometric Wide swath,IW)模式下的單視復(fù)數(shù)影像和SRTM DEM數(shù)據(jù),采用像素跟蹤算法和層流理論進(jìn)行冰川厚度反演與冰儲(chǔ)量估算,以認(rèn)識(shí)喀喇昆侖地區(qū)冰川厚度空間分布及冰厚與氣候之間的響應(yīng)關(guān)系;將估算結(jié)果與法國南部比利牛斯天文臺(tái)數(shù)據(jù)服務(wù)中心的全球冰川厚度數(shù)據(jù)產(chǎn)品進(jìn)行對(duì)比分析,從而驗(yàn)證本文方法的可行性。研究結(jié)果可為研究喀喇昆侖地區(qū)冰川厚度變化提供參考。

1 研究區(qū)概況和數(shù)據(jù)來源

喀喇昆侖山是全球山地冰川最發(fā)達(dá)的高大山脈,從阿富汗最東部向東南延伸約480 km,平均海拔超過6 000 m,山區(qū)終年低溫,年平均氣溫低于0 ℃,由于冬春兩季受到西風(fēng)大環(huán)流的影響,冰川成冰區(qū)域年降水量超過1 000 mm,獨(dú)特的地理環(huán)境使得喀喇昆侖山有超過37%的面積被冰川覆蓋。在當(dāng)前氣候變暖的形勢下,全球絕大多數(shù)冰川出現(xiàn)萎縮現(xiàn)象,而喀喇昆侖地區(qū)冰川對(duì)氣候變化響應(yīng)不夠顯著,在某些情況下出現(xiàn)末端前進(jìn)、面積增加的現(xiàn)象,被稱為“喀喇昆侖異?！盵9],本文選取喀喇昆侖地區(qū)卻哥隆瑪冰川、錫亞琴冰川和巴爾托洛冰川3條規(guī)模較大且具代表性的冰川作為研究對(duì)象。

為研究喀喇昆侖地區(qū)冰川厚度分布,采用6景Sentinel-1A衛(wèi)星降軌單視復(fù)數(shù)影像,時(shí)間跨度為2019年1月—2020年1月,影像距離向和方位向空間分辨率分別約為5、20 m,觀測范圍覆蓋研究區(qū)目標(biāo)冰川,影像數(shù)據(jù)全部來自歐洲航天局(European Space Agency,ESA) (https://search.asf.alaska.edu),具體影像信息見表1所列。另外,美國地質(zhì)調(diào)查局90 m空間分辨率的SRTM DEM數(shù)據(jù)主要用于研究區(qū)地形位移量補(bǔ)償、地理編碼及坡度計(jì)算;第6版蘭多夫冰川編目 (Randolph Glacier Inventory 6.0,RGI 6.0)用于提取研究區(qū)冰川矢量邊界,數(shù)據(jù)來源于全球陸地冰空間測量計(jì)劃(http://www.glims.org/RGI/);冰川厚度數(shù)據(jù)產(chǎn)品(https://doi.org/10.6096/1007)[10]作為驗(yàn)證數(shù)據(jù)用于結(jié)果的精度評(píng)定。

表1 研究區(qū)影像信息

2 研究方法

為獲取目標(biāo)冰川流速,并在其基礎(chǔ)上獲取冰川厚度信息,本研究采用像素跟蹤算法和層流理論進(jìn)行冰川流速提取與冰川厚度的反演,通過與全球冰川厚度產(chǎn)品進(jìn)行對(duì)比分析,驗(yàn)證該方法的反演精度,相應(yīng)技術(shù)流程如圖1所示。圖1中,NCC表示歸一化互相關(guān)(normalized cross-correlation)。

圖1 冰川厚度反演技術(shù)流程

2.1 像素跟蹤算法

對(duì)于流速較快的山地冰川,像素跟蹤算法具有良好的適用性[11],本研究利用覆蓋研究區(qū)的典型冰川多期Sentinel-1A衛(wèi)星遙感影像,采用改進(jìn)的像素跟蹤算法進(jìn)行位移量計(jì)算,該方法以NCC算法為基礎(chǔ),通過主、輔影像的非冰川區(qū)穩(wěn)定像素的亞像素級(jí)配準(zhǔn)結(jié)果精確解算主、輔影像之間的幾何變換模型,進(jìn)而估計(jì)冰川表面運(yùn)動(dòng)速度分布。NCC的計(jì)算準(zhǔn)則公式為:

(1)

通過上述復(fù)雜的匹配過程,能夠獲取研究區(qū)總體位移量,總體位移量可分為冰川運(yùn)動(dòng)位移量和非冰川位移量,其計(jì)算公式為:

D=Dglacier+Dnonglacier

(2)

其中:D為研究區(qū)總體位移量;Dglacier為冰川運(yùn)動(dòng)位移量;Dnonglacier為非冰川位移量,其中包含合成孔徑雷達(dá)(synthetic aperture radar,SAR)影像軌道不同引起的全局位移量、電離層位移量、噪聲位移量和地形高低起伏產(chǎn)生的位移量。為去除SAR影像軌道不同引起的全局位移量,采用隨機(jī)抽樣一致性算法來估算變換模型,獲得變換模型參數(shù)去除軌道相關(guān)位移量;由于喀喇昆侖地區(qū)處于中低緯度,而電離層引起的位移量對(duì)于高緯度地區(qū)較為明顯,因此電離層影響可以忽略不計(jì);為去除圖像配準(zhǔn)過程中產(chǎn)生的噪聲點(diǎn),在圖像后處理過程中采用自適應(yīng)中值濾波進(jìn)行濾波處理;冰川區(qū)地形陡峭引起的位移量受到空間基線長度的影響,采用DEM數(shù)據(jù)和基線信息對(duì)產(chǎn)生的地形位移量進(jìn)行補(bǔ)償。在考慮非冰川位移量的基礎(chǔ)上進(jìn)行位移合成和地理編碼,即可得到冰川區(qū)表面流速分布。

2.2 層流理論

采用層流理論進(jìn)行冰川厚度的反演[12],該方法的基本原理是綜合考慮冰川坡度、基底滑動(dòng)、自身重力等因素對(duì)冰川厚度的影響,根據(jù)輸入?yún)?shù)反演冰川厚度,如圖2所示。

圖2 層流理論示意圖

圖2中:H為冰厚;τ為基底剪應(yīng)力;F為冰體沿斜面方向向下滑動(dòng)的力,即重力沿斜面的分力;Flim為冰體靜止時(shí)的最大摩擦限度;m為冰體質(zhì)量。隨著冰體質(zhì)量增加,當(dāng)斜面方向向下滑動(dòng)的力超過摩擦限度,就會(huì)產(chǎn)生運(yùn)動(dòng)。

假設(shè)冰體底面積為S,冰體在斜面固定的力為τS,基底剪應(yīng)力等于冰體重力沿斜面方向的分力,即

τS=mgsinα,m=ρHS

(3)

其中:g為重力加速度;α為冰川表面坡度[6];ρ為冰密度。

冰-巖界面條件變化產(chǎn)生的底部摩擦力不可忽略,需要增加形狀因子f′,則有:

τ=f′ρgHsinα

(4)

其中,f′與山谷冰川橫斷面結(jié)構(gòu)有關(guān),山谷冰川f′取值范圍[13]為0.5～0.9。

根據(jù)層流理論,冰川表面運(yùn)動(dòng)速率us計(jì)算公式為:

(5)

其中:ub為冰川基底運(yùn)動(dòng)速率;A為蠕變參數(shù),與環(huán)境溫度、冰體雜質(zhì)等因素相關(guān)[14];n為格倫定律指數(shù),在1.5～4.2范圍內(nèi)變動(dòng),與冰川所受的應(yīng)力有關(guān)[15]。

結(jié)合式(4)、式(5)可以得出:

(6)

其中,k為基底滑動(dòng)系數(shù)。

3 數(shù)據(jù)處理與結(jié)果分析

3.1 冰川流速提取與結(jié)果分析

利用2019年1月—2020年1月的6景降軌SAR數(shù)據(jù),采用像素跟蹤算法進(jìn)行冰川流速提取,設(shè)置主影像模板窗口的大小為96,輔影像搜索窗口的大小為32,步長為2,可以在冰川區(qū)達(dá)到較高的配準(zhǔn)精度。為削弱噪聲的影響,采用自適應(yīng)中值濾波進(jìn)行降噪處理,設(shè)置最小濾波半徑為7,最大濾波半徑為15,得到目標(biāo)冰川平均流速如圖3所示。圖3中,冰川運(yùn)動(dòng)方向如紅色箭頭所示。

圖3 3條冰川流速分布圖

由于冰川底部基巖本身的差異及冰川受山谷拖滯作用影響,冰川的形態(tài)呈現(xiàn)多樣性、復(fù)雜性。從圖3可以看出:卻哥隆瑪冰川和巴爾托洛冰川最大速率位于冰川上游(約為38 cm/d),隨著冰川坡度不斷減緩,冰川流速逐漸放緩,但均呈上游積累區(qū)相對(duì)于下游消融區(qū)流速較快、冰川中心相對(duì)于冰川邊緣流速較快的特點(diǎn),這是冰體下滑過程中與周邊山體發(fā)生摩擦導(dǎo)致速度減緩,另外,表面冰磧物的阻熱特性一定程度上減緩了冰川的消融、阻礙冰川運(yùn)動(dòng);而錫亞琴冰川末端出現(xiàn)流速激增的情況(約為40 cm/d),相對(duì)于上游流速明顯增大,期間可能發(fā)生冰川躍動(dòng),這也是喀喇昆侖地區(qū)躍動(dòng)冰川常見現(xiàn)象[16],此外冰川流速變化符合“上游向下游逐漸減緩、中央向兩側(cè)逐漸減緩”的運(yùn)動(dòng)規(guī)律。

3.2 冰川厚度反演與結(jié)果分析

為分析冰川厚度分布,在獲取流速的基礎(chǔ)上,采用層流理論進(jìn)行冰川厚度反演,其中蠕變系數(shù)A=2.4×10-24/(Pa3·s),滑動(dòng)系數(shù)k=0.2,冰密度ρ=917 kg/m3,重力加速度[10]g=9.81 m/s2;綜合考慮喀喇昆侖山脈所處氣候條件及冰川結(jié)構(gòu)等影響因素,將格倫定律指數(shù)n取值為2.6,形狀因子f′取值為0.6。目標(biāo)冰川厚度分布如圖4所示。圖4中,a1—a2、b1—b2、c1—c2表示冰川剖面線始末位置。

圖4 3條冰川厚度反演圖

由圖4可知:卻哥隆瑪冰川的冰厚最大值出現(xiàn)在支流與干流交匯口,約為380 m,并且大部分區(qū)域厚度為220～250 m,占冰川面積的60.04%;錫亞琴冰川最大冰厚出現(xiàn)在冰川干流下游,約為380 m,大部分區(qū)域冰厚為180～220 m,占冰川面積的42.35%;而巴爾托洛冰川冰厚最大值集中在冰川干流上游,隨著海拔不斷升高,冰厚在210～300 m范圍內(nèi)不斷增加,且冰川大部分區(qū)域冰厚為170～200 m,占冰川面積的57.09%?？傮w而言,3條冰川均呈從上游到下游逐漸減薄、支流向干流不斷增厚的趨勢。

為進(jìn)一步分析地形因素對(duì)冰川厚度的影響,本文沿著圖4中的剖面線提取3條冰川的冰厚并進(jìn)行擬合,得到剖面圖,如圖5所示。

圖5 3條冰川干流高程和冰厚變化曲線

從圖5可以看出,卻哥隆瑪冰川冰厚存在2個(gè)峰值,分析其原因,有以下2點(diǎn):① 可能是冰川支流推動(dòng)作用的影響,使得上游和支流的物質(zhì)在支流交匯處釋放,造成支流交匯處冰體較厚;② 可能是由于冰川的搬運(yùn)作用,使得冰磧物在冰川下游拐角處堆積,一定程度上阻礙了冰體的運(yùn)動(dòng)和消融,促進(jìn)物質(zhì)積累。

錫亞琴冰川由于冰下水壓力的增加導(dǎo)致冰川末端基底滑動(dòng)增強(qiáng)與冰形變的加速[16],致使上游積累區(qū)物質(zhì)向下游推進(jìn),出現(xiàn)冰川下游相對(duì)上游較厚的現(xiàn)象。

巴爾托洛冰川為典型的樹狀型冰川,主干呈東西走向,冰川整體坡度較緩,厚度自東向西逐漸減薄,與高程變化一致,冰川下游海拔較低、氣溫較高等因素促進(jìn)冰雪消融并抑制冰雪積累[17]。

3.3 冰儲(chǔ)量估算與統(tǒng)計(jì)分析

為進(jìn)一步估算冰川儲(chǔ)量,將冰川厚度柵格影像重采樣至50 m分辨率,與冰川邊界矢量疊置分析,得到完整的冰川區(qū),利用RGI 6.0全球冰川編目數(shù)據(jù),通過統(tǒng)計(jì)冰川區(qū)的柵格像元數(shù)目,將柵格像元大小與每個(gè)像元厚度值的乘積之和作為冰川區(qū)冰儲(chǔ)量。卻哥隆瑪冰川、錫亞琴冰川和巴爾托洛冰川統(tǒng)計(jì)面積分別為283.08、1 053.88、781.06 km2,平均厚度值分別為237.27、229.66、172.39 m,卻哥隆瑪冰川、錫亞琴冰川、巴爾托洛冰川的冰儲(chǔ)量估算結(jié)果分別為67.42、242.18、149.72 km3。

為分析上述結(jié)果的可靠性,采用文獻(xiàn)[10]中的法國南部比利牛斯天文臺(tái)數(shù)據(jù)服務(wù)中心的冰川厚度產(chǎn)品作為驗(yàn)證數(shù)據(jù),該產(chǎn)品經(jīng)過精度檢驗(yàn),已被應(yīng)用于全球冰川流速與厚度特征分析等研究中。本文反演結(jié)果與厚度產(chǎn)品對(duì)比見表2所列。

表2 3條冰川厚度、冰儲(chǔ)量本文結(jié)果與厚度產(chǎn)品對(duì)比

從表2可以看出,卻哥隆瑪冰川、錫亞琴冰川、巴爾托洛冰川的平均厚度與厚度產(chǎn)品的誤差分別為14.05%、13.86%、11.31%,而冰儲(chǔ)量的誤差分別為14.05%、14.30%、13.85%,表明本文冰厚反演方法具有較好的可行性,可用于輔助分析青藏高原冰川厚度時(shí)空特征差異。

4 討 論

本文采用層流理論進(jìn)行冰川厚度反演,該方法對(duì)于光學(xué)影像和雷達(dá)影像的流速結(jié)果都具有良好的適用性,但模型中不同參數(shù)的組合易導(dǎo)致“異參同效”現(xiàn)象的發(fā)生。本研究根據(jù)冰川表面流速、坡度、形狀因子等模型參數(shù)的不確定性,考慮冰川不同區(qū)域的氣象條件差異,結(jié)合不同山谷冰川所處氣候類型、外界應(yīng)力及橫斷面結(jié)構(gòu)等屬性的差異,設(shè)置合理的參數(shù),盡可能減少“異參同效”現(xiàn)象,使得該方法獲得較為理想的結(jié)果。

模型反演存在不確定性,其中輸入數(shù)據(jù)、模型參數(shù)是造成反演結(jié)果不確定性的主要來源。

1) 根據(jù)式(6)可得,若冰川表面流速增加10倍,則會(huì)對(duì)冰厚高估超過75%[10],會(huì)導(dǎo)致冰厚模擬值在冰川的活躍階段與靜止階段之間不一致,這進(jìn)一步增加了厚度值的不確定性。

2) 對(duì)于山地冰川,可變參數(shù)τ和形狀因子f′的取值分別與冰川谷壁橫向阻力和冰川橫截面的形狀有關(guān),但實(shí)測值與理論值往往有很大差異[18]。在區(qū)域尺度上,由于冰川類型、氣候、地形等因素的顯著差異,且統(tǒng)一特征參數(shù)難以確定以及缺少有關(guān)物質(zhì)平衡的區(qū)域性數(shù)據(jù),該方法在區(qū)域研究中具有較大的局限性。

3) 雖然本文中冰川流速有具體的時(shí)間信息,但由于DEM數(shù)據(jù)的匱乏,無法獲取目標(biāo)時(shí)間段的DEM數(shù)據(jù),采用了2003年公開發(fā)布的SRTM DEM數(shù)據(jù)進(jìn)行坡度提取,時(shí)間的不一致最終也會(huì)直接導(dǎo)致冰厚反演結(jié)果的誤差[19]。

目前對(duì)于冰體流動(dòng)的相關(guān)研究中,關(guān)于冰體在給定應(yīng)力作用下的變形規(guī)律尚認(rèn)識(shí)不清,且冰體的變形速率與所處溫度和外界施加的壓力密切相關(guān)。雖然眾多研究者設(shè)計(jì)了盡可能接近冰川的冰厚實(shí)驗(yàn)來建立冰厚模型,但是實(shí)驗(yàn)條件均為理想狀態(tài),相對(duì)于自然條件存在一定差別。

影響冰川厚度變化的因素有很多,降水和溫度在氣象因子中占有重要地位[20]。喀喇昆侖地區(qū)降水較為充沛,但3條冰川所在的區(qū)域年均降水量有較大差異,卻哥隆瑪冰川年均降水量約為1 000 mm,而錫亞琴冰川、巴爾托洛冰川年均降水量分別約為600、700 mm,并且卻哥隆瑪冰川、錫亞琴冰川、巴爾托洛冰川年平均氣溫分別為-13.79、-22.93、-20.61 ℃,降水量和氣溫的差異會(huì)影響冰川消融與積累的轉(zhuǎn)化效率,進(jìn)而影響冰川的補(bǔ)給與冰儲(chǔ)量的變化,增加了冰厚反演結(jié)果的不確定性。

5 結(jié) 論

本文利用6景Sentinel-1A衛(wèi)星影像和90 m 空間分辨率SRTM DEM數(shù)據(jù),應(yīng)用層流理論模型對(duì)卻哥隆瑪冰川、錫亞琴冰川、巴爾托洛冰川3條冰川進(jìn)行厚度反演及冰儲(chǔ)量估算,分析冰厚空間分布主要特點(diǎn)及其原因,結(jié)論如下:

1) 基于冰川表面流速和層流理論能夠反演復(fù)雜地形區(qū)山地冰川的厚度信息,通過與法國南部比利牛斯天文臺(tái)數(shù)據(jù)服務(wù)中心的冰川厚度數(shù)據(jù)產(chǎn)品對(duì)比驗(yàn)證了結(jié)果的有效性。本文方法可為青藏高原地區(qū)山地冰川厚度反演與冰儲(chǔ)量估算研究提供參考。

2) 卻哥隆瑪冰川、錫亞琴冰川、巴爾托洛冰川的平均厚度分別為237.27、229.66、172.39 m,同一地區(qū)不同冰川厚度略有差異,冰川的搬運(yùn)作用和堆積作用與冰厚的空間分布存在著必然的聯(lián)系,但冰川整體上呈從上游到下游逐漸減薄、支流向干流不斷增厚的趨勢。

3) 層流理論物理意義明確,相對(duì)于傳統(tǒng)經(jīng)驗(yàn)公式更具有穩(wěn)健性,基于層流理論進(jìn)行冰厚估算時(shí),由于模型過于簡化,估算結(jié)果誤差偏大[15],實(shí)測資料的稀缺一定程度上限制了該方法結(jié)果的驗(yàn)證,如何有效利用豐富的冰川編目資料及其他輔助數(shù)據(jù)改進(jìn)厚度估算模型是以后進(jìn)一步研究的方向。