青藏高原納木錯湖冰物候變化遙感監(jiān)測與模擬

吳艷紅，郭立男，范蘭馨，文夢宣，遲皓婧，張兵

1.中國科學(xué)院空天信息創(chuàng)新研究院 數(shù)字地球重點(diǎn)實(shí)驗(yàn)室,北京100094;

2.中國科學(xué)院大學(xué),北京100049

1 引 言

湖冰物候（湖泊的凍結(jié)和消融時間）直接受湖泊能量平衡變化的影響，對區(qū)域氣候變化有很好的指示作用（Kouraev等，2007；Marszelewski和Skowron，2006）。湖冰物候在反映氣候變化方面某種程度上可能比氣溫更可靠（Livingstone，1997）。湖冰物候的變化影響著湖泊生態(tài)系統(tǒng)中生物的繁殖以及漁業(yè)、運(yùn)輸業(yè)等人類生產(chǎn)活動。因此，近年來中高緯度湖冰物候?qū)θ蜃兣捻憫?yīng)越來越為研究人員所關(guān)注（Salonen等，2009；Bengtsson，2011；Fu和Yao，2015；Gou等，2017；Latifovic和Pouliot，2007；Magnuson 等，2000）。比如，Magnuson 等（2000）利用實(shí)測數(shù)據(jù)分析了150年間北半球39 個湖泊的湖冰物候變化趨勢，發(fā)現(xiàn)北半球湖冰物候表現(xiàn)出結(jié)冰日延后和融冰日提前的趨勢，并為后續(xù)的研究所證實(shí)（Bernhardt 等，2012；Brown 和Duguay，2010；Ghanbari等，2009）。

基于AVHRR、MODIS 等光學(xué)傳感器（如：Latifovic和Pouliot，2007；Kropá?ek等，2013；Duguay等，2015）和微波遙感器（如：QuickSCAT，SSM/I，SMMR，AMSR-E）等對湖冰物候進(jìn)行動態(tài)監(jiān)測已成為當(dāng)前湖冰物候研究的重要手段。其中，光學(xué)遙感湖冰物候監(jiān)測常用的方法包括閾值法、指數(shù)法、反射率光譜曲線法等（龐毓雯等，2020；Guo等，2018）。在微波遙感方面，邱玉寶等（2017）利用微波亮溫數(shù)據(jù)分析了青藏高原湖冰的凍融變化。

遙感是獲取湖冰物候信息的有效途徑。然而，由于遙感平臺的服務(wù)年限相對較短，難以揭示湖冰物候的長期變化規(guī)律。另外，遙感信息反映了湖冰物候變化的事實(shí)而無法解釋其變化的物理機(jī)制。認(rèn)識湖泊凍融過程的物理機(jī)制，進(jìn)而量化湖冰物候?qū)夂虻捻憫?yīng)，通常則需借助數(shù)值模擬的方法來實(shí)現(xiàn)（任曉倩等，2014）。數(shù)理模型考慮了多種氣候強(qiáng)迫項和湖泊形態(tài)特征對湖冰凍融過程的作用，有助于從機(jī)理的層面揭示和預(yù)測湖冰的生消過程。現(xiàn)有湖泊過程模式如MIMLAKE96（Fang 和Stefan，1996）、LAKEoneD （Joehnk 和Umlauf，2001）、LAKE模式（Stepanenko等，2016）、SIMSTRAT （Goudsmit 等，2002） CLIM （Duguay等，2003）、LIMNOS（Vavrus等，1996）、HIGHTSI（Launiainen 和Cheng，1998）等。這些模型在世界各地區(qū)已經(jīng)得到了應(yīng)用和驗(yàn)證。中國湖泊研究起步較晚，現(xiàn)有模型多為外國學(xué)者針對當(dāng)?shù)丨h(huán)境開發(fā)，在青藏高原的適用性仍然有待驗(yàn)證（宋興宇等，2020）。方楠等（2017）改進(jìn)并評估了WRF模型在納木錯的適用性；蘇東生等（2018）利用“中國區(qū)域高時空分辨率地表氣象驅(qū)動數(shù)據(jù)集（CMFD）”驅(qū)動Flake 模型模擬青海湖熱力狀況對氣候變化的響應(yīng)；Huang 等（2019）評價了3 個一維湖泊模式在納木錯湖泊的適用性并通過敏感性試驗(yàn)對模型進(jìn)行率定。

青藏高原地區(qū)受人類活動影響較小，其氣溫升溫速率大于同期全球陸表氣溫的變化速率（Huang等，2017），是全球變化的敏感區(qū)域（杜娟等，2020）。由于歷史觀測數(shù)據(jù)的缺乏，青藏高原湖冰物候的變化研究還不夠充分。本文利用過程模型（LAKE 2.3）重建青藏高原納木錯湖泊長時序（1962年—2018年）的湖面溫度和湖冰物候，并將模擬結(jié)果與遙感監(jiān)測結(jié)果進(jìn)行對比驗(yàn)證。在此基礎(chǔ)上，分析近60 a 來納木錯湖冰物候的變化趨勢。本文所用的LAKE 2.3 模型（Stepanenko 等，2016）是綜合性的一維湖泊模型，涉及湖泊的熱力學(xué)過程以及生物地球化學(xué)過程。該模型已有許多應(yīng)用，但在青藏高原湖泊研究中的應(yīng)用還較少。

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

2.1 研究區(qū)概況

納木錯位于30°30′N—30°55′N，90°16′E—91°03′E（圖1），海拔4718 m，湖泊面積1981 km2。納木錯是西藏自治區(qū)第二大湖，也是世界上海拔最高的湖泊之一。湖泊集水面積10610 km2，2005年測得的最大深度超過90 m（Wang 等，2009）。納木錯所處流域內(nèi)冰川占比約1.5%（Wu 等，2014）。近幾十年來，隨著青藏高原氣溫的升高，納木錯湖泊的面積、水位、水溫、冰期以及生物化學(xué)過程也發(fā)生著深刻的變化（Guo 等，2020；Huang 等，2017；Zhang等，2011）。納木錯屬于高原亞寒帶季風(fēng)半干旱氣候帶，輻射強(qiáng)烈，年日照時數(shù)達(dá)2900—3200 h，年均溫1.9 ℃，年降水量約500 mm，湖面年蒸發(fā)量約1200 mm。

圖1 納木錯及鄰近氣象站點(diǎn)地理位置Fig.1 Locations of the Nam Co and the nearest meteorological station

2.2 遙感數(shù)據(jù)

本文首先基于多源遙感信息應(yīng)用不同方法提取納木錯的湖冰物候信息。主要的遙感數(shù)據(jù)源包括：ARC-Lake（V3），MODIS 以及基于微波輻射計（AMSR-E、AMSR2 和FY-3 MWRI）的湖泊亮溫信息。

ARC-Lake（V3）是基于沿軌道掃描輻射計（ATSR2/AATSR）獲取的湖泊表面溫度日尺度數(shù)據(jù)集（http：//www.laketemp.net/home_ARCLake/［2021-04-30］），包括日間和夜間的湖面溫度。本研究以日間溫度和夜間溫度的均值作為納木錯湖面溫度日值，以此獲取的湖面溫度序列作為模擬湖面溫度的驗(yàn)證數(shù)據(jù)之一，時間范圍1995-06-01—2012-04-05。

中分辨率成像光譜儀MODIS（moderate resolution imaging spectroradiometer）搭載于兩顆極軌衛(wèi)星Terra和Aqua上。Terra上午自北向南穿過赤道，Aqua下午穿過赤道。雙星系統(tǒng)1—2 d 就可以覆蓋全球一次。MODIS 傳感器包含36 個波段，已經(jīng)擁有應(yīng)用于各個領(lǐng)域的多種產(chǎn)品。本文應(yīng)用的是其L3級日尺度陸面溫度產(chǎn)品（MOD11A1），空間分辨率為1 km。數(shù)據(jù)集包含日間和夜間的陸表溫度及發(fā)射率。為計算湖泊面上的平均溫度，本文首先將納木錯湖泊矢量邊界圖向內(nèi)縮進(jìn)1 km 以避免湖岸混合像元的影響。然后，以日間溫度和夜間溫度的均值作為湖表的日均溫度，從而獲取湖表日均溫度序列。為保證數(shù)據(jù)集之間的可比性，本文選擇的MODIS產(chǎn)品的時間跨度與ARC-Lake一致。

此外，本文還采用了邱玉寶等（2017）生產(chǎn)的“2002年—2016年高亞洲地區(qū)中大型湖泊微波亮溫和凍融數(shù)據(jù)集”（http：//www.csdata.org/p/49/［2021-04-30］）。該數(shù)據(jù)集依據(jù)微波輻射計（AMSR-E、AMSR2和FY-3 MWRI）觀測像元內(nèi)湖泊和陸表的面積比例，并應(yīng)用混合像元分解方法獲取像元（亞像元級）的湖泊亮溫信息，是通過被動微波遙感對像元級湖冰凍融進(jìn)行監(jiān)測的一種方法。

2.3 LAKE模型驅(qū)動數(shù)據(jù)

LAKE 模型需要的輸入包括氣溫、比濕、氣壓、風(fēng)速、短波輻射、長波輻射和降水等氣象要素。短波輻射和向下的長波輻射根據(jù)緯度、實(shí)際日照時數(shù)、氣溫和水汽壓以計算而得。其他氣象變量以距納木錯最近的當(dāng)雄站（30°29′N，91°06′E，海拔4200 m）的氣象數(shù)據(jù)為依據(jù)。基礎(chǔ)氣象數(shù)據(jù)來自中國氣象數(shù)據(jù)網(wǎng)，包括1963年—2019年近60 a的氣溫、降水、風(fēng)速、氣壓、相對濕度和日照時數(shù)。

考慮到氣溫隨高程的變化，本文首先以0.65℃/100 m 的氣溫大氣直減率（胡芩 等，2014）對氣象站氣溫Tstation進(jìn)行高程校正得到湖泊氣溫Tlake，即

式中，Hstation為氣象站海拔，Hlake為湖泊海拔。其次，利用氣壓階公式對氣象站氣壓Pstation進(jìn)行高程校正，獲取湖面氣壓Plake，即

式中，比濕q的計算公式為

式中，ε= 0.622，e為水汽壓，是相對濕度RH 與飽和水汽壓ew的乘積，而飽和水汽壓可用Bolton公式（Bolton，1980）計算：

3 方 法

3.1 湖冰物候指標(biāo)

本文基于遙感與模擬提取的湖冰物候指標(biāo)包括：開始結(jié)冰日FUS（湖面由完全是水到開始有像元變?yōu)楸鶎?yīng)的日期）、完全融化日BUE（湖面由存在冰到所有像元為水對應(yīng)的日期），以及湖泊結(jié)冰期（ID），即BUE與FUS二者之間的天數(shù)。這3個指標(biāo)與湖氣界面的水熱交換過程密切相關(guān)。

3.2 湖冰物候遙感監(jiān)測

本文參照Guo等（2018）的方法，基于不同的遙感信息源提取納木錯2000年—2015年的湖冰物候，并分析其變化規(guī)律。湖冰物候的遙感監(jiān)測方法具體包括：（1）利用MOD09GA 產(chǎn)品生成全湖平均的反射率時間序列，通過閾值法在反射率時間序列上提取湖冰物候（M1）；（2）利用MOD10A1冰雪產(chǎn)品，計算湖面水體占湖泊面積的比，以水占 比10% 為 物 候 提 取 閾 值（M2）；（3） 利 用MOD11A1 產(chǎn)品生成全湖平均湖面溫度時間序列，通過閾值法提取湖冰物候（M3）。上述3 種方法的詳細(xì)過程可參考Guo 等（2018）。此外，本文從“2002年—2016年高亞洲地區(qū)中大型湖泊微波亮溫和凍融數(shù)據(jù)集”中獲取納木錯2002年—2016年的湖冰物候作為對比數(shù)據(jù)（M4）。

3.3 基于LAKE模型的湖冰物候模擬

本文采用的LAKE模型（Stepanenko等，2016）其熱力學(xué)過程表述為

式中，T為水溫，垂直坐標(biāo)z以湖面為z=0，垂直向下為正，h為湖泊深度，并定義為ξ=z/h為坐標(biāo)。cw是水的比熱容，ρ是水的密度，λ渦流擴(kuò)散率，是湖泊表面水量平衡，S是太陽輻射通量，穿透到深度z，M是浮力混合（對流）。穿透到深度z的太陽輻射計算：

式中，αe是消光系數(shù)。

湖泊水溫的上邊界條件是湖泊表面能量平衡：

式中，S是太陽輻射，Ea向下的大氣輻射，Es向上的長波輻射，H和LE分別是顯熱和潛熱通量，α是表面反照率。

本研究中用到的是LAKE 2.3 模型。LAKE 模型包含的參數(shù)較多，但參數(shù)大多具有明確的物理意義。模擬過程中，除了湖泊基本形態(tài)參數(shù)（湖泊經(jīng)緯度、湖泊面積、湖泊深度等），本文參考2005年—2007年實(shí)地考察的納木錯湖盆形態(tài)估算出各個深度的水平截面面積，作為湖泊形態(tài)參數(shù)輸入模型。在模擬中，本文取納木錯的平均深度為40 m，水體的消光系數(shù)為0.1，冰的消光系數(shù)為100。主要的模型參數(shù)設(shè)定如表1。這些參數(shù)取值與Huang等（2019）的研究中對納木錯消光系數(shù)設(shè)置相一致。

表1 納木錯LAKE模型主要參數(shù)設(shè)定Table 1 Model parameter configuration for the Nam Co

本文設(shè)定LAKE 模型時預(yù)熱時間為20 a，并以模型最終輸出的表層溫度（skin temperature）的年內(nèi)變化特征獲取湖泊結(jié)冰和融化的日期。

4 結(jié)果分析與討論

4.1 基于遙感的湖冰物候

圖2 給出了納木錯2003年—2014年基于不同遙感監(jiān)測方法的湖冰物候?？v軸是年積日（DOY）（從當(dāng)年01-01 起開始計算的天數(shù)），大于365 d 的表明物候日期發(fā)生在第2年。由圖2 中可以看出，納木錯湖泊結(jié)冰日多始于12 月底至次年1 月，次年4—5 月湖冰則完全融化。多年平均情況下，開始結(jié)冰日（FUS）為01-07，完全化凍日（BUS）為03-10。年均湖泊結(jié)冰期（ID）的長度大致為124 d，標(biāo)準(zhǔn)差為10 d。研究時段內(nèi)，湖冰物候的變化總體趨勢并不顯著。

圖2 遙感監(jiān)測的湖冰物候（每個柱狀圖的上下端分別對應(yīng)湖冰物候中的FUS和BUE；M1—M4分別為對應(yīng)的4種遙感監(jiān)測方法）Fig.2 Lake ice phenology from remote sensing（The upper and lower ends of each histogram correspond to FUS and BUE in lake ice phenology respectively；M1—M4 represents the four different methods used）

不同方法之間的比較可以看出，基于微波亮溫數(shù)據(jù)（M4）獲取的冰期通常比其他方法監(jiān)測所得的冰期短，而其他3種遙感方法在湖冰物候監(jiān)測上的差異沒有明顯的規(guī)律。這一結(jié)果表明，受遙感數(shù)據(jù)時空分辨率的影響（如：數(shù)據(jù)缺失），湖冰物候的遙感監(jiān)測還存在較大的不確定性。

4.2 湖冰物候模擬

由于遙感的數(shù)據(jù)序列相對較短，難以反映納木錯湖冰物候近幾十年來的變化規(guī)律，因此本文借助LAKE模型以重建納木錯的湖冰物候序列，并揭示湖冰物候的長期變化規(guī)律及其對氣候變化的響應(yīng)。

圖3 首先給出了LAKE 模型湖面溫度的模擬結(jié)果。由圖3中可以看出，湖泊表面溫度的模擬值與遙感監(jiān)測的湖面溫度一致性較高，模擬結(jié)果與ARC-Lake 水溫數(shù)據(jù)的納什效率系數(shù)（Nash-Sutcliffe efficiency coefficient）達(dá)0.82，均方根誤差為1.75℃。基于LAKE 模型的湖冰物候模擬結(jié)果如圖4 所示。由圖4 可以看出，模擬所得的3 個物候指標(biāo)中，融化日期BUE 與遙感監(jiān)測的結(jié)果一致性最高，均方根誤差為8.5 d。模擬的開始結(jié)冰日FUS普遍早于遙感監(jiān)測的結(jié)果，模擬的結(jié)冰期總體上大于遙感監(jiān)測的結(jié)冰期。

圖3 基于模型和遙感的湖面溫度對比Fig.3 Comparison of lake surface temperature based on modelling and remote sensing

圖4 模擬湖冰物候與遙感湖冰物候（每個柱狀圖的上下端分別對應(yīng)湖冰物候）Fig.4 Comparison of lake ice phenology based on modelling and remote sensing（The upper and lower ends of each histogram correspond to FUS and BUE in lake ice phenology respectively）

LAKE 模擬結(jié)果與遙感監(jiān)測結(jié)果之間的差異，一方面受模型的不確定性的影響（包括模型結(jié)果、參數(shù)以及輸入的不確定性）；另一方面也受遙感監(jiān)測的不確定性的影響。特別是遙感監(jiān)測的湖面溫度通常為衛(wèi)星過境的瞬時值（或日間夜間瞬時值的均值），而不能真正代表一天內(nèi)的平均溫度，且MODIS 監(jiān)測的表面溫度受云等影響波動較大。此外，確定湖冰物候時相關(guān)閾值的設(shè)定也會影響湖冰物候的判別。因此，雖然遙感和模擬都可以獲取湖冰物候信息，但湖冰物候監(jiān)測與模擬的準(zhǔn)確性還需要進(jìn)一步的地面觀測加以驗(yàn)證。

4.3 近60年來湖冰物候變化趨勢

LAKE 等有物理基礎(chǔ)的數(shù)值模型刻畫了氣候因子與湖冰物候之間的定量關(guān)系，有助于揭示湖冰物候?qū)夂蜃兓捻憫?yīng)?；贚AKE模型，本文重建了近60 a來納木錯的湖冰物候如圖5所示。

圖5 1962年—2018年納木錯湖冰物候變化趨勢（FUS、BUE、ID分別為開始結(jié)冰日、完全融化日和冰期）Fig.5 Trends in lake ice phenology during the period 1962—2018 of Nam Co（FUS，BUE，and ID are date of freeze-up start，break-up end and ice duration respectively）

基于Mann-Kendall（MK）非參數(shù)檢驗(yàn)法（Kendall，1955；Mann，1945）的湖冰物候變化趨勢分析表明，近60 a 來納木錯湖泊結(jié)冰推遲、融化提前、冰期縮短，且所有物候指標(biāo)的變化趨勢均達(dá)到95%的顯著性水平。結(jié)冰日期（FUS）平均每10 a延后3.3 d，融化日期（BUE）每10 a提前3.0 d。冰期（ID）的變化速率最快，以6.4 d/（10 a）的速率縮短。

4.4 湖冰物候?qū)夂蜃兓拿舾行?/h3>

已有研究表明，至21世紀(jì)末，青藏高原氣溫相對于1980年—1999年將升高2℃（Chen和Frauenfeld，2014）。為了定量分析湖冰物候?qū)夂蜃兓捻憫?yīng)特征，本文基于LAKE模型開展了湖冰物候?qū)夂蜃兓拿舾行苑治觥Ｔ诿舾行苑治鲋?，假設(shè)納木錯氣溫以0.1℃的間隔持續(xù)升高至2℃，而后基于LAKE 模型的模擬結(jié)果，計算分析湖冰物候隨氣溫升高相較于歷史時期（1963年—2018年）的可能變化，結(jié)果見圖6。

圖6 湖冰物候?qū)鉁刈兓拿舾行裕‵US、BUE、ID分別為開始結(jié)冰日、完全融化日和冰期）Fig.6 Sensitiveness of lake ice phenology in response to increasing air temperature（FUS，BUE，and ID are date of freeze-up start，break-up end and ice duration respectively）

模擬結(jié)果表明，當(dāng)氣溫升高達(dá)0.7℃時，開始結(jié)冰日（FUS）顯著推遲。如果氣溫升高2℃，結(jié)冰日期平均推遲7.3 d。湖冰融化過程對氣溫的響應(yīng)比結(jié)冰過程更顯著。當(dāng)氣溫升高2℃，融化日期平均提前12.4 d。因此，冰期隨著溫度升高將顯著縮短，平均縮短19.7 d。湖冰物候的變化會影響湖泊水量平衡、能量平衡過程、湖泊及周圍的生態(tài)系統(tǒng)和人類活動。因此，在氣候變暖的背景下，需要對湖冰物候變化情況進(jìn)行持續(xù)的監(jiān)測和預(yù)測。

5 結(jié) 論

本研究以納木錯湖泊為研究區(qū)，基于4種不同的遙感監(jiān)測方式獲取了2002年—2015年納木錯的湖冰物候指標(biāo)，結(jié)果表明各方法之間的年際變化情況比較一致，但基于微波亮溫數(shù)據(jù)獲取的冰期通常比其他方法監(jiān)測的時間短。

在遙感監(jiān)測的基礎(chǔ)上，利用一維湖泊模型（LAKE）重建了1962年—2018年納木錯湖冰物候序列，揭示了近60 a 來納木錯湖冰物候的變化趨勢及其對氣候變化的響應(yīng)特征。研究結(jié)果表明，近60 a 來，納木錯結(jié)冰推遲、融化提前、冰期縮短（6.7 d/（10 a））。湖冰物候（特別是融冰日期）對氣候變化十分敏感。在氣溫升高2℃的氣候變化情景下，未來納木錯冰期將縮短19.7 d。

本研究將遙感與數(shù)值模擬相結(jié)合，發(fā)揮兩種手段的優(yōu)勢，為湖冰物候研究提供了新的思路和方法。但是，在青藏高原湖冰物候的研究中，現(xiàn)有的遙感方法和數(shù)值模擬方法都存在不同程度的不確定性。在未來的研究中，除了需要進(jìn)一步提高遙感信息的精度和分辨率，完善模型的結(jié)構(gòu)和參數(shù)化方案外，在條件適宜的地方開展持續(xù)的湖冰物候地面觀測，將有助于提高青藏高原湖冰物候監(jiān)測與模擬水平。