周靜恬,邱玉寶,黃 琳,Juha LEMMETYINEN,石利娟, 李青寰, 施建成
(1.中國科學(xué)院空天信息創(chuàng)新研究院數(shù)字地球重點(diǎn)實(shí)驗室,北京 100094;2.可持續(xù)發(fā)展大數(shù)據(jù)國際研究中心,北京 100094;3.中國科學(xué)院大學(xué),北京 100049;4.中國科學(xué)院空天信息創(chuàng)新研究院-芬蘭氣象研究所北極觀測聯(lián)合研究中心(JRC-AO),芬蘭索丹屈萊FI-99660;5.芬蘭氣象研究所北極空間中心,芬蘭赫爾辛基FI-00560;6.中國科學(xué)院國家空間科學(xué)中心,北京 100190)
積雪是冰凍圈最為活躍的季節(jié)性要素,是全球氣候變化的靈敏指示器[1]。針對積雪的監(jiān)測具有重要的現(xiàn)實(shí)意義,當(dāng)前被動微波遙感在積雪監(jiān)測方面具有不可替代的作用。積雪微波輻射傳輸模型模擬可開展過程和機(jī)理研究,仿真環(huán)境可重現(xiàn)積雪演變、積雪與微波相互作用等重要過程,是被動微波積雪參數(shù)——雪深和雪水當(dāng)量反演算法發(fā)展的基礎(chǔ)。
針對積雪的物理特性是影響空間監(jiān)測的重要因素,其演變特性存在較大的空間異質(zhì)性[2-3],積雪物理特性的演化影響微波輻射[4-5],從而影響雪深/雪水當(dāng)量反演算法[6-7]。當(dāng)前計算雪深/雪水當(dāng)量的算法主要是半經(jīng)驗線性算法[8-9],其不確定性來源于雪物理特性變化及下墊面變化和大氣變化影響[4],其中雪的物理特性演變過程對算法影響較大[10]。由于積雪反演算法在時間上和空間上還存在很大的不確定性,導(dǎo)致精度受到影響[7,11-12]。通過改進(jìn)算法參數(shù)來更好捕捉動態(tài)積雪演變過程,或耦合發(fā)射率模型和陸表模型,或耦合發(fā)射率模型和由天氣驅(qū)動的積雪演變模型將很可能提高雪深/雪水當(dāng)量反演算法準(zhǔn)確性[13]。結(jié)合實(shí)驗測量的雪深等輔助信息[14]或了解積雪的先驗知識[15]有助于減少雪深/雪水當(dāng)量反演算法的不確定性。
為應(yīng)對上述挑戰(zhàn),論文采用在北歐地區(qū)的2009—2013年積雪地面和地基微波實(shí)驗觀測數(shù)據(jù)集[16-17],研究不同積雪期的分期判斷方法,探討了芬蘭積雪實(shí)驗場積雪特性的時間序列變化以及地基微波輻射計觀測亮溫差與雪深變化關(guān)系,最后采用MEMLS模型[18-20]開展積雪微波輻射的模擬、對比和分析研究。
北歐積雪實(shí)驗(Nordic Snow Radar Experiment,NoSREx)是ESA Earth Explorer 7候選任務(wù)CoReH2O[21](寒冷地區(qū)水文學(xué)高分辨率觀測)中Phase A研究組成部分,該實(shí)驗?zāi)康氖窃诒睒O寒帶森林地區(qū)的代表性地點(diǎn)提供整個冬季連續(xù)時間序列的地面、主被動微波積雪觀測。實(shí)驗區(qū)域位于芬蘭氣象研究所北極研究中心Sodankyl?的觀測站,地理位置為67.3618°N,26.6338°E(圖1)。實(shí)驗采用了包括10.65、18.7、37、90 GHz四個頻率的微波輻射計系統(tǒng),其中在2011—2012年實(shí)驗期間頻段21 GHz替換了90 GHz,四個入射角分別為30°、40°、50°、60°,獲得了水平和垂直極化下的觀測亮溫以及標(biāo)準(zhǔn)差。地面觀測數(shù)據(jù)包括積雪參數(shù)的人工和自動測量,人工測量是指每隔3~5天進(jìn)行雪坑(snow pit)觀測,其參數(shù)包括雪深、密度和雪水當(dāng)量、雪密度剖面、雪溫度剖面、雪層評估,觀測天數(shù)在2009—2013四年依次為91、31、23、32天。自動測量參數(shù)包括溫度、露點(diǎn)溫度、風(fēng)速、氣壓、雪深、地表濕度和溫度、氣溫、雪水當(dāng)量等,每天共測量8次,每次間隔3小時,觀測天數(shù)在2009—2013年依次為365、365、366、365天[16]。
圖1 NoSREx IOA的網(wǎng)絡(luò)攝像頭圖像以及主要微波儀器,SnowScat散射計和SodRad輻射計系統(tǒng)的照片[16]Fig.1 Webcam image of NoSREx IOA and photographs of main microwave instruments,the SnowScat scatterometer and the SodRad radiometer system[16]
該數(shù)據(jù)集已用于積雪相關(guān)模型模擬[22-23]、反演積雪參數(shù)[24-25]、改進(jìn)算法[26]、耦合雪物理和輻射傳輸模型[27-29]、分析觀測數(shù)據(jù)[29-31]等研究。
多層積雪微波輻射傳輸模型(MEMLS)是由Wiesmann等[18-19]開發(fā)的一個針對多層積雪的被動微波輻射傳輸模型,適用頻率范圍為5~100 GHz。它以輻射傳輸為基礎(chǔ),將積雪根據(jù)物理特性分為多層,利用六流近似理論來描述每個雪層內(nèi)部的多次散射與吸收,同時考慮了雪層之間的界面散射。由于模型在特定環(huán)境下開發(fā),在不同環(huán)境下的適應(yīng)性不同,為了將其有效應(yīng)用到芬蘭積雪試驗場,在該研究中需要調(diào)整模型的輸入?yún)?shù)地表均方根高度參數(shù)的值。
1.3.1 積雪積累期、穩(wěn)定期和消融期確定
根據(jù)氣溫和積雪厚度的變化,把積雪期劃為3個雪期,其中積雪“積累期”為出現(xiàn)降雪后積雪逐漸累積的過程,表現(xiàn)為雪深隨時間呈現(xiàn)顯著增加趨勢,氣溫基本已降低至0℃以下;積雪“穩(wěn)定期”為氣溫連續(xù)低于0℃時,雪深隨時間增加趨勢減小,即積雪壓實(shí)的過程,表現(xiàn)為雪深出現(xiàn)增加明顯變緩或甚至下降趨勢;積雪“消融期”為隨著溫度升高且連續(xù)5日大于0℃,積雪進(jìn)入融化或消融,可能伴隨再凍結(jié)過程,表現(xiàn)為雪深迅速減小,氣溫基本維持在0℃以上。
針對芬蘭Sodankyl?地區(qū),根據(jù)其地面觀測數(shù)據(jù)的溫度和雪深組合特征(圖2),研判3個時期開始和結(jié)束時間,方法為:積累期開始時間為雪深一周內(nèi)連續(xù)大于5 cm,氣溫處于從0℃以上過渡到0℃以下階段,此前氣溫不低于-20℃;積累期結(jié)束時間也是穩(wěn)定期開始時間,為連續(xù)3個周內(nèi)雪深變化范圍為±5 cm以內(nèi),氣溫至少連續(xù)5天低于0℃;穩(wěn)定期結(jié)束時間也是消融期開始時間,為連續(xù)5天氣溫大于0℃,雪深處于減小狀態(tài),一般為50~70 cm;消融期結(jié)束時間為雪深一周內(nèi)連續(xù)小于5 cm,氣溫處于0℃以下過渡到0℃以上階段。
圖2 2009—2013年氣溫和雪深時序變化圖Fig.2 Changes of temperature and snow depth from 2009 to 2013
1.3.2 模型輸入確定
MEMLS模型輸入包括頻率、入射角、天空背景亮度溫度、散射系數(shù)、積雪參數(shù)、地表溫度以及積雪-土壤界面反射率[19]。天空背景的亮度溫度所需的大氣透射率參照Pulliainen等[32-33]的提出的統(tǒng)計反演方法,天空背景亮度溫度所需的上下行大氣亮度溫度由Aschbache[34]提出的公式估算。散射系數(shù)參照Matzler等[18]的工作,采用玻恩近似計算。輸入的積雪參數(shù)包括積雪層數(shù)、積雪溫度、積雪濕度、積雪密度、積雪鹽度、積雪厚度,指數(shù)相關(guān)長度,指數(shù)相關(guān)長度采用Durand等[35]的經(jīng)驗公式通過積雪粒徑估算,其他參數(shù)由NoSREx實(shí)驗的地面測量獲得[16]。陸地積雪的鹽度通常認(rèn)為是0[36-37],本次模擬也將積雪鹽度設(shè)為0。地表溫度由NoSREx實(shí)驗的地面測量數(shù)據(jù)獲得。水平極化和垂直極化的積雪-土壤界面反射率,本研究是利用Wang&Choudhury(1981)半經(jīng)驗?zāi)P停ê喎QQHN)計算[38-39]。QHN模型中所需的土壤介電常數(shù)采用Dobson模型計算[40]。土壤-雪界面的粗糙程度通過土壤粗糙度來反映[41-42]。土壤粗糙度通常用均方根高度和表面相關(guān)長度兩個統(tǒng)計變量表示,其具體定義可見文獻(xiàn)[43-44]。土壤粗糙度較難直接測量,通常采用最小化成本函數(shù)(cost function)[45-48]的方法估算。本研究選用的代價函數(shù)如下:
式中:CF為代價函數(shù)的值;m為模擬次數(shù);TM為測量亮溫;TS為模擬亮溫。
氣溫和雪深隨時間變化如圖2所示。氣溫日波動性和月波動性較大,年波動性相似,一般在11—12月降低到0℃以下,在4—5月上升到0℃以上;每年雪深的變化趨勢相似,一般從10—11月逐漸累積,在3—4月達(dá)到最大值,在4月下旬至5月迅速融化。根據(jù)圖2以及本文1.3.1章節(jié)的分類標(biāo)準(zhǔn),判定芬蘭Sodankyl?的積雪期分為積累期、穩(wěn)定期和消融期,積累期一般為10月—次年2月,時長約為4~5個月;穩(wěn)定期為2—4月,時長約為2個月;消融期為4—5月,時長約為1個月,其具體結(jié)果時間如表1。
表1 不同雪期分類的時間(年-月-日)Table1 Time classification of different snow peroid(yyyy-mm-dd)
圖3展示了自然積雪的分層特性。如圖可知,融化狀態(tài)(Melt Forms,MF)主要出現(xiàn)在積累期早期(10—12月)和消融期(4—5月)。降水粒子(Precipitation Particles,PP)和分解碎片降水粒子(Decomposing and Fragmented precipitations particles,DF)主要出現(xiàn)在近雪表層(圖中均為綠色),且主要在4月前出現(xiàn)。在積累期后期和穩(wěn)定期,積雪分層顆粒形狀出現(xiàn)的種類主要是圓形顆粒(Rounded Grains,RG)、片狀顆粒(Faceted Crystals,F(xiàn)C)、深霜(Depth Hoar,DH)。在消融期主要類別為融化狀態(tài)。圖3(c)可見深藍(lán)色較少,表明2011—2012年缺乏深霜層,這與該年較其他年份更溫和的氣溫有關(guān),色帶參考季節(jié)性陸面積雪的國際分類[49]。
結(jié)合圖3~4可知:(1)在整個雪季,研究區(qū)域的垂直剖面上每個積雪層粒徑不同,一般底層積雪粒徑最大且為深霜層,即深霜層粒徑最大,底層粒徑會從小變大再變小,表層粒徑一直較小且變化不大。(2)在積累期早期,10—12月,雪深較淺時,雪的粒徑也較小,一般在1.5 mm以內(nèi)。(3)隨著雪深的積累,在積雪底層的粒徑逐漸增大。在積累期早期10—12月從較小的0.25~1 mm,在積累期后期(1—2月)增長為1~3 mm。(4)隨著雪深增加速度減緩,在積累期后期和穩(wěn)定期(1—4月),近地表的粒徑維持在2 mm左右,均出現(xiàn)了增加到至少2.5 mm的現(xiàn)象。(5)在融雪期4—5月,觀測雪層的顆粒形狀大多數(shù)為融雪狀態(tài),積雪融化,粒徑相較穩(wěn)定期的表層大底層小,為1 mm左右。(6)靠近積雪表面的粒徑值始終較小,基本維持在1.5 mm以內(nèi),但少數(shù)情況下會出現(xiàn)較大值。(7)觀測的平均粒徑的最大值出現(xiàn)在每年的穩(wěn)定期,2—3月,值為2.5~4 mm,均出現(xiàn)在近地表層。
圖3 2009—2013年自然積雪的分層特性Fig.3 Stratification characteristics of snow cover from 2009 to 2013
為研究線性亮溫梯度算法(18 GHz和37 GHz,V和H)在北歐實(shí)驗區(qū)域的適用性,分析了微波亮溫差對雪深變化的依賴性(圖5)。積累期(12月—次年2月)和穩(wěn)定期前期(2—3月),雪深和亮溫差(18~37 GHz)整體來看具有一致的變化趨勢,局部來看具有相反的波動性,即雪深小幅度減小亮溫差反向增加,如圖5中矩形框中的部分:黑線和鋸齒形狀的紅色和藍(lán)色曲線反向波動,這可能是由于融雪使得亮溫差減小,而后積雪重新凍結(jié)亮溫差增大。在穩(wěn)定期后期(3—4月),由于3月雪深均大于60 cm,亮溫差對雪深增加不再敏感趨于飽和,體現(xiàn)了算法的局限性[21];局部反向現(xiàn)象較之前不明顯,這可能是由于早期融雪形成的雪殼結(jié)構(gòu)逐漸松弛為更典型的冬末降雪[31],雪粒徑和雪顆粒形狀均改變。在消融期(4—5月),隨著積雪迅速融化雪深減小,亮溫差值波動性較大基本維持在±15 K以內(nèi),亮溫差與雪深線性相關(guān)性不明顯。
圖5 2009—2010年雪深和亮溫差(18 GHz和37 GHz)在50°入射角條件下隨時間的變化Fig.5 Variation of snow depth and brightness temperature difference(18 GHz and 37 GHz)with time at 50°incident angle from 2009 to 2010
雪深的變化和亮溫差變化在積累期和穩(wěn)定期總體上具有相似性,但局部來看沒有明顯的一致性,這可能是由于積雪亮溫還會受到如含水量、積雪粒徑以及雪顆粒形狀等的影響,而這些參數(shù)均隨著時間變化。由上一章可知整個雪季積雪演變較大,因此需在算法中考慮隨時間演化的積雪特征,已有研究考慮積雪動態(tài)變化的算法[8]。微波亮溫差對雪深和線性函數(shù)在不同時期依賴性不同,進(jìn)一步表明提高反演算法精度需考慮參數(shù)的動態(tài)變化。
模型模擬所需要的土壤粗糙度參數(shù)較難直接測量,且其對模型模擬結(jié)果有一定影響,在模擬前需要確定土壤粗糙度參數(shù)的值。地表均方根高度的值由土壤本身的性質(zhì)決定,不同情況下該數(shù)值可以在數(shù)毫米到幾十毫米之間變化[39]。采用1.3.2章節(jié)所述方法,得出的不同頻率和入射角的最優(yōu)值(表2),大多數(shù)情況下,在值為0.001 m時均方根誤差最小,因此在模擬時將垂直和水平極化下的地表粗糙度值均設(shè)為0.001 m。潘金梅等[18-20]的研究表明Sodankyl?地區(qū)的地表均方根高度大約為0.001 m,與本研究的估算結(jié)果吻合。
表2 地表均方根誤差參數(shù)最優(yōu)值Table 2 Optimal values of surface root mean square error parameters
圖4 2009—2013年雪季積雪粒徑演變Fig.4 Evolution of snow grain size in snow season from 2009 to 2013
調(diào)整參數(shù)值后,基于MEMLS模型模擬2009—2013年的微波輻射亮溫,用觀測亮溫和模擬亮溫的均方根誤差評估模擬的準(zhǔn)確性。如表3所示,在10.65 GHz頻率下,入射角度為50°和60°,極化方式為垂直時模擬結(jié)果最好,RMSE分別為4.79 K和4.72 K。其次模擬結(jié)果較好的是18.7 GHz頻率下,入射角為50°和60°,極化方式為垂直極化,RMSE分別為7.13 K和7.25 K。頻率為37 GHz時,模擬結(jié)果較差,當(dāng)水平極化入射角為60°時RMSE達(dá)到21.76 K。模擬結(jié)果最差的是在90 GHz下,極化方式為垂直極化,入射角為50°和60°,RMSE分別為22.37 K和22.14 K。因此在低頻波段的垂直極化的模擬結(jié)果較好,高頻模擬結(jié)果較差。
表3 模擬和觀測亮度溫度的均方根誤差Table 3 Root mean square errors of simulated and observed brightness temperatures
比較不同積雪期的模擬準(zhǔn)確度(表4),在10.65 GHz且垂直極化下,除入射角為40°,穩(wěn)定期的模擬亮溫的RMSE均小于其他兩時期;在18.7 GHz且垂直極化下,在入射角為40°和60°,穩(wěn)定期的模擬亮溫的RMSE小于其他兩時期;在37 GHz且垂直極化下,MEMLS模型在穩(wěn)定期的模擬效果較積累期和消融期更好;在90 GHz的模擬結(jié)果較差,模擬亮溫和觀測亮溫的RMSE均大于20 K。將2009—2013年模擬出的亮溫和觀測亮溫進(jìn)行了對比分析(圖6)??傮w上,水平極化較垂直極化的觀測和模擬亮溫波動性更大;10.65 GHz和18.7 GHz相較37 GHz和90 GHz的模擬亮溫更為接近觀測亮溫;穩(wěn)定期較積累期和消融期的模擬亮溫更為接近觀測亮溫。
圖6 2009—2013年入射角50°、垂直(V)和水平極化(H)下10.65、18.7、37 GHz和90 GHz的測量和模擬亮度溫度Fig.6 Measured and simulated brightness temperature at 10.65,18.7,37 GHz and 90 GHz at incident angle 50°,vertical(V)and horizontal polarization(H)from 2009 to 2013
表4 不同積雪期模擬和觀測亮度溫度的均方根誤差Table 4 Root mean square errors of simulated and observed brightness temperatures in different snow period
依據(jù)2009—2013年北歐積雪觀測數(shù)據(jù),分析了芬蘭Sodankyl?研究區(qū)域的雪季的不同時期的積雪物理演化過程,并利用MEMLS模型對積雪的微波輻射亮溫進(jìn)行了模擬,分析認(rèn)為:基于溫度和雪深的變化情況,可將整個雪季分為積累期(10月—次年2月)、穩(wěn)定期(2—4月)、消融期(4—5月)。
地面觀測數(shù)據(jù)分析表明,積雪物理演化過程有如下3點(diǎn)特征:1)氣溫和雪層溫度在穩(wěn)定期達(dá)到最低,雪深在穩(wěn)定期末積累到最厚;2)積累期和穩(wěn)定期顆粒形狀主要為圓形粒徑、片狀顆粒和深霜,深霜粒徑較大;3)積累期早期會出現(xiàn)融化狀態(tài)和降水粒子,消融期積雪顆粒形狀主要為融化狀態(tài)。
從3個積雪期與雪深、亮溫差的關(guān)系來看,積累期(10月—次年2月)和穩(wěn)定期前期(2—3月),雪深和亮溫差(18~37 GHz)整體來看變化趨勢一致,在穩(wěn)定期后期(3—4月),雪深均超過60 cm,亮溫差趨于飽和。在消融期(4—5月),積雪迅速融化雪深減小,亮溫差(18和37 GHz)基本維持在±15 K以內(nèi)且波動性較大,積雪微波輻射亮溫差和雪深兩者的關(guān)系隨著雪季的不同時期變化,難以用靜態(tài)半經(jīng)驗反演算法進(jìn)行描述。這可能是由于積雪物理特性本身(從而微波輻射)隨時間具有較大的變異性,后續(xù)可考慮針對不同積雪期的積雪演化動態(tài)過程來改進(jìn)積雪反演算法。
MEMLS模型輸入的地表均方根高度參數(shù)設(shè)為0.001 m[20]較為適合芬蘭實(shí)驗環(huán)境,校正后的MEMLS模型在較低頻率(10.65 GHz、18.7 GHz),較高入射角(50°和60°)且垂直極化下能更好的模擬地基輻射計觀測亮溫。3個雪期均在低頻(10.65 GHz和18.7 GHz)的垂直極化下,模擬結(jié)果較好。在穩(wěn)定期的10.65 GHz、50°入射角且垂直極化下模擬結(jié)果最佳,RMSE最低為2.49 K。對于37 GHz且垂直極化下,穩(wěn)定期的模擬效果較積累期和消融期更好,這表明根據(jù)氣溫和雪深變化來分類積雪期,可更好的表征積雪演化對微波輻射信號的影響。模擬結(jié)果在低頻段(10.65 GHz和18.7 GHz)較好,可考慮結(jié)合較低頻段來改進(jìn)反演算法。
MEMLS模型模擬和觀測亮溫的誤差是由模型本身的局限性以及觀測參數(shù)誤差等因素造成的。模擬結(jié)果受各個輸入?yún)?shù)的影響,不同參數(shù)對模擬結(jié)果影響程度不同[50]。驅(qū)動MEMLS模型模擬的參數(shù)通過積雪剖面調(diào)查獲取,該調(diào)查為了盡可能詳細(xì)地捕捉積雪變化,由專家每隔3~5天且連續(xù)四個冬季開展雪坑(snow pit)測量。根據(jù)該高質(zhì)量的積雪剖面調(diào)查結(jié)果進(jìn)行模擬,模擬和觀測亮溫仍存在誤差,這表明一些影響積雪輻射傳輸過程的重要變量可能在模型中并沒有得到有效的體現(xiàn),因此有必要開展新的積雪模型研究。此外,地面觀測參數(shù)在模型模擬和積雪反演中能夠起到重要的指導(dǎo)作用,但為模型和反演算法確定合理的積雪結(jié)構(gòu)等效參數(shù)從而讓模型更好的反映真實(shí)情況同樣重要。雖然輻射計觀測點(diǎn)和積雪刨面調(diào)查點(diǎn)的空間距離非常近,但積雪結(jié)構(gòu)的空間變異性依然會導(dǎo)致模擬的結(jié)果和輻射計觀測結(jié)果出現(xiàn)差異,這在一定程度上反映出將傳統(tǒng)地面調(diào)查技術(shù)獲取的積雪結(jié)構(gòu)信息直接運(yùn)用于模型模擬的局限性。而且,積雪剖面調(diào)查獲取積雪結(jié)構(gòu)信息的過程中會存在調(diào)查者主觀因素的影響(例如如何確定積雪中每一層的邊界以及如何選取粒徑的測量位置),這也會導(dǎo)致輻射計觀測結(jié)果和模擬結(jié)果出現(xiàn)偏差。一般而言,高頻波段比低頻波段對積雪粒徑更敏感,而水平極化比垂直極化對積雪分層結(jié)構(gòu)更敏感。本次模擬結(jié)果低頻優(yōu)于高頻,且低頻(10.65 GHz和18.7 GHz)的垂直極化優(yōu)于水平極化,高頻(90 GHz)的水平極化優(yōu)于垂直極化,這說明積雪結(jié)構(gòu)信息在模型模擬中的重要性。一方面,積雪剖面調(diào)查取樣過程中的細(xì)微偏差足以對模擬結(jié)果產(chǎn)生顯著影響;另一方面,模型模擬中的理想化條件和自然條件下積雪的狀況存在區(qū)別,物理理論和實(shí)際觀測的有效連接還迫切的需要更深入的研究。