魏英楠， 馬 龍， 孫柏林， 張 晶

（內蒙古農業(yè)大學水利與土木建筑工程學院，內蒙古 呼和浩特 010018）

氣候變暖導致冰川持續(xù)退縮，高緯度地區(qū)、高海拔山區(qū)多年凍土層融化［1-2］，水資源供需矛盾加劇，對全球和區(qū)域尺度上降水量、徑流量等水文要素產(chǎn)生深刻影響。大興安嶺地處高緯寒旱地區(qū)，是我國半濕潤區(qū)與半干旱區(qū)、季風區(qū)與非季風區(qū)的分界線［3］。獨特的氣候條件和地理位置導致其成為全球氣候變化敏感地帶之一，揭示該地區(qū)長歷史時期水文變化特征的研究亟待開展。然而器測數(shù)據(jù)往往少于百年，極大地阻礙了對歷史水文變化的深入了解。

樹木年輪定年準確、易于獲取及分辨率高等優(yōu)勢［4］使其成為延長水文資料的重要手段之一。自20世紀30 年代美國內華達州首次利用樹輪延長歷史年徑流量資料后［5］，樹木年輪在全球范圍被廣泛應用于重建長時間尺度水文資料。北美［6-8］和歐洲地區(qū)［9-11］研究成果較多，國內主要集中在新疆［12-14］、青藏高原［15］及大興安嶺［3,16］等水文氣候要素限制作用顯著的高緯度寒旱區(qū)域。大興安嶺地區(qū)樹輪年代學的研究起步較晚，多以興安落葉松（Larix gmelinii）和樟子松（Pinus sylvestris）為代用材料，例如分析興安落葉松和樟子松徑向生長與氣候變化的響應及差異，發(fā)現(xiàn)氣候變化改變了樹木生長與氣候間的響應關系，樟子松受氣候變暖的影響更加顯著［17-18］；基于樟子松寬度年表重建伊敏河1868—2002 年年徑流總量［19］，重現(xiàn)了伊敏河過去135 a 經(jīng)歷的6 個干旱時期和4 個濕潤時期；利用樟子松首次建立海拉爾河整個流域及11個支流過去202～216 a年徑流量序列［20］，發(fā)現(xiàn)海拉爾徑流量的變化與厄爾尼諾-南方濤動（El Ni?o-Southern Oscillation，ENSO）密切相關，并預測隨著ENSO 活動的增強，海拉爾河流域將面臨極端水文事件。

白音敖包自然保護區(qū)地處大興安嶺山地南麓，位于東北、華北和蒙古植物區(qū)系交匯地帶［21］，原始林木資源豐富，且處于半干旱地區(qū)，發(fā)源于大興安嶺最高峰的貢格爾河是保護區(qū)內主要河流，水分條件對該區(qū)域樹木徑向生長限制作用顯著［12］，在此開展樹木年代學研究具有得天獨厚的優(yōu)勢。紅皮云杉（Picea koraiensis）是松科云杉屬常綠喬木，耐蔭耐寒耐干旱，在不同立地條件下均能生長［22］，被列入《世界自然保護聯(lián)盟》（IUCN）瀕危物種［23］，是保護區(qū)內樹齡較長、樹輪邊界明顯且對水文氣候要素響應敏感的珍貴樹種，具有較高的樹木年輪學研究價值。前人利用紅皮云杉重建了白音敖包過去160余年的降水量序列［24］，但關于徑流量歷史變化的研究較少。基于此，本文利用紅皮云杉樹輪寬度年表，重建大興安嶺南麓近2 個世紀的徑流量歷史序列，研究結果重現(xiàn)了大興安嶺南麓貢格爾河流域無數(shù)據(jù)記載時期的徑流量變化，為其今后水文相關研究提供參考。

1 資料與方法

1.1 研究區(qū)概況

采樣點白音敖包國家級自然保護區(qū)（圖1）地理位置介于大興安嶺南端和錫林郭勒高平原之間［21］，地形為大興安嶺山地向蒙古高原的過渡地帶。氣候類型屬于大陸性寒溫帶半干旱森林草原氣候，寒暖分明，水熱同期［24-25］。貢格爾河橫臥于保護區(qū)北側，境內流經(jīng)長度14 km，敖包河是其主要支流。

圖1 研究區(qū)采樣點和水文氣象站分布Fig.1 Distributions of sampling point and hydrometeorological stations in the study area

1.2 水文與氣候數(shù)據(jù)

本文使用的氣溫、降水量資料為林西氣象站1951—2016 年逐月平均、平均最高、平均最低氣溫和降水量（圖2a），數(shù)據(jù)來源于中國國家氣象信息中心。徑流量資料為距離研究區(qū)最近的錫林浩特水文站1972—2016年逐月徑流量（圖2b），數(shù)據(jù)來源于當?shù)厮目睖y局水文年鑒。

圖2 研究區(qū)多年月平均氣溫、降水量與徑流量變化Fig.2 Changes in monthly average temperature,precipitation and runoff for many years in the study area

研究區(qū)具有明顯的大陸性氣候特征，冬夏氣溫較差大，降水量年內分配極不均勻。多年平均氣溫4.73 ℃，7 月氣溫值最高，1 月氣溫為全年最低值。多年平均降水量約為373.17 mm，呈單峰分布，7 月降水量為全年最高，5—9 月降水量占全年的90%。多年平均徑流量約為0.17×108m3，4 月徑流量最大，12月至次年2月為枯水期。

1.3 樹木年輪采樣及年表建立

年輪樣本采集時間為2017 年8 月，采樣時嚴格遵循國際樹木年輪數(shù)據(jù)庫（ITRDB）標準［26］，共采集37 株樹的80 根樣芯。待樹輪樣本在陰涼通風處自然晾干后，使用由粗至細目數(shù)不同的砂紙逐級打磨，直至年輪邊界清晰可見。利用WinDENDROTM年輪分析系統(tǒng)［27］測定每一生長輪的日歷年。使用COFECHA 程序［28］檢驗交叉定年和寬度測量結果，對缺輪、假輪、定年錯誤等問題進行檢查，最終78根樹輪樣本用于年表的構建。

樹木的徑向生長不僅受自身遺傳因子的控制，還與其生長過程中所處環(huán)境密切相關。為了解樹木承載的水文氣候變化信息，需去除樣本自身生長趨勢［4］。使用ARSTAN 程序［29］負指數(shù)函數(shù)［30］去除樹木本身遺傳因子及競爭干擾產(chǎn)生的生長趨勢，最后利用雙權重平均法合成標準化年表（圖3a）。

圖3 紅皮云杉標準化年表特征Fig.3 Standardized chronological characteristics of Picea koraiensis

基于樹木年輪學研究的復本原理，利用子樣本信號強度（Subsample signal strength，SSS）確定年表的最低復本量（圖3b）［31］，為保證重建序列可靠性并保留年表最大有效長度，限定SSS＞0.85 時為年表的可靠時段。

通過分析年表的統(tǒng)計特征值［32］，可以判斷其是否適用于樹木年輪學研究。年表的主要統(tǒng)計特征值見表1。平均敏感度（MS）是判斷一個年表質量優(yōu)劣的重要指標之一，其數(shù)值越大，樹木徑向生長對周圍環(huán)境變化越敏感。一般認為MS＞0.2 的序列適合進行樹木年輪的研究［33］；標準差的數(shù)值越大，樹木群體受到氣候因素的控制越劇烈；一階自相關系數(shù)反映氣候條件對樹木生長的持續(xù)性影響，其數(shù)值越小年表質量就越高；信噪比是用來衡量樣本中包含環(huán)境信息量大小的統(tǒng)計量，其數(shù)值較高時，年表中含有較多氣候信息，適用于進行樹輪氣候學研究；樣本總體代表性數(shù)值越大，樣本的代表性就越好，一般樣本總體代表性＞0.85 時所建立的年表比較可靠。紅皮云杉年表統(tǒng)計特征值均滿足樹木年輪學標準，說明年表質量較好，包含較多氣候信息，可用于樹木年輪氣候學的研究之中。

表1 紅皮云杉年表主要統(tǒng)計特征Tab.1 Main statistical indicators of Picea koraiensis standardized chronology

1.4 研究方法

（1）使用Pearson 法對年表與氣溫、降水量、徑流量等要素進行相關分析［34］。

（2）采用逐一剔除法檢驗重建模型的穩(wěn)定性和可靠性。

（3）為研究徑流量變化趨勢，利用Z-score 法對重建徑流量歷史序列進行標準化處理，定義標準化后的徑流量連續(xù)2 a 以上位于0m3之上的時段為變豐期，連續(xù)2 a 以上位于0m3之下的時段為變枯期。計算公式如下：

式中：R為徑流量重建值（108m3）；mean 為重建徑流量序列的平均值（108m3）；SD 為重建徑流量序列標準差（108m3）。

（4）使用變異系數(shù)（Cv）反映徑流量歷史重建序列變化的劇烈程度，計算公式如下：

（5）定義重建徑流量高于mean+σ的年份定義為豐水年，低于mean-σ的年份定義為枯水年，連續(xù)經(jīng)歷2 a及以上豐水年的時期為豐水期，連續(xù)經(jīng)歷2 a以上枯水年的時期為枯水期。其中：

式中：σ為重建徑流量序列的標準差（108m3）；Xi為第i年對應的徑流量（108m3）；μ為重建徑流量序列的平均值（108m3）；N為樣本數(shù)量（個）。

（6）徑流量重建序列的周期變化分析，使用Morlet 小波分析法［35］進行分析，在Matlab 中，其計算公式為：

式中：ψ(t)為小波函數(shù)；t為時間（年）。

Morlet小波伸縮尺度（a）與周期（T）有如下關系：

式中：ω0為無量綱頻率，通常取6.2附近的經(jīng)驗值。

（7）采用Mann-Kendall 檢驗法對徑流量突變年份進行確定。

2 結果與分析

2.1 樹木徑向生長與水文氣候要素間響應關系的分析

考慮到樹木的生長在10月底之前已經(jīng)結束，冬季水文氣象因子對樹木年輪生長具有滯后效應［36］，因此利用皮爾遜相關法分析樹輪寬度年表與研究區(qū)前一年11—12 月和當年1—10 月水文氣候資料的響應關系（圖4）。

圖4 年表與水文氣候要素的響應關系Fig.4 Response relationship between chronology and hydroclimatic elements

除降水量外，樹輪寬度年表與當年的徑流量、氣溫關系更加密切。年表與各月平均氣溫、平均最高氣溫和平均最低氣溫均呈現(xiàn)顯著負相關關系，即過高的氣溫不利于紅皮云杉的生長。年表與6月氣溫負相關關系極顯著，夏季研究區(qū)氣溫處于全年最高值，樹木正值生長季旺期，此時氣溫升高，加強了樹木的呼吸、蒸騰作用，導致土壤中本就匱乏的水分不能滿足樹木的生長所需，從而起到負面效果［37］；年表與降水量表現(xiàn)為正向響應關系，但除前一年11月以外，其余月份的降水量均未能通過0.05 的顯著性檢驗，可能是降水年內分配極不均勻，對樹木徑向生長限制作用微弱。

紅皮云杉寬度年表與徑流量呈顯著正相關關系，3—10 月徑流量通過了0.05 顯著性檢驗，所以年表對徑流量的響應更加敏感。徑流并非直接作用于樹木生長，而是通過氣溫、降水等因子與其連接起來［38］。研究區(qū)屬于半干旱地區(qū)，水分對紅皮云杉生長過程起促進作用，但12 月研究區(qū)氣溫較低，樹木處于休眠狀態(tài)，此時水分過多將導致土壤凍結，樹根發(fā)生凍害，不利于樹木生長［39］；春季氣溫回升，氣候適宜樹木生長，樹木開始萌發(fā)幼芽，但根系尚處于休眠狀態(tài)，吸收功能較弱，徑流可促進根系吸收水分；夏季樹木進入生長旺盛時期，形成層細胞活躍，分裂活動旺盛，充足的水分條件有利于植物進行光合作用，形成寬輪；秋天氣溫下降，樹木蒸騰作用減弱，準備進入休眠狀態(tài)，此時充足的水分有助于植物為次年的生長儲存能量［40］。

2.2 徑流量重建

基于年表與水文氣候要素的響應分析結果，結合紅皮云杉的生長特性，本次選擇以3—10 月組合月徑流量為因變量，樹輪寬度年表為自變量，利用冪指數(shù)回歸方程重建了研究區(qū)1845—2016 年徑流量序列，重建方程如下：

式中：yi為第i年組合月徑流量（108m3）；xi為樹輪寬度標準化年表序列；r為組合月徑流量與年表的相關系數(shù)；N為月徑流量實測資料長度（年）；R2adj為調整后r2；F為重建方程方差檢驗；P為重建方程顯著性值。

利用國際樹木年輪研究中的逐一剔除法對此重建模型進行檢驗［30］。逐一剔除法主要采用的統(tǒng)計參數(shù)有相關系數(shù)（r）、誤差縮減值（RE）、有效系數(shù)（CE）、符號檢驗、一階差符號檢驗、乘積平均值檢驗（t）等。其中RE大于0.3代表擬合程度較好，越接近于1，方程可靠性越高［41］；符號檢驗利用每對樣本的數(shù)據(jù)差值符號進行檢驗，可檢測出樣本間是否存在顯著性差異；t不但包含了對符號的檢驗并增加了序列數(shù)值的檢驗，能夠直接反映出重建值是否包含有水文氣候信息，其值越大表明2 個序列越相近。重建序列和實測序列的r=0.815，一階差r=0.769，達到了95%的置信水平。以1972—1994 年為校準期、1994—2016 年為檢驗期對重建方程進行檢驗，校準期r=0.677，檢驗期r=0.881，均通過了0.01 顯著性檢驗；重建方程的RE=0.411，CE=0.640，二者遠遠大于0，說明重建序列可靠；重建方程的原序列符號檢驗（32+/13-）通過了0.05 顯著性檢驗，但一階差符號檢驗（23+/21-）未能通過，說明方程對低頻變化的重建更加準確［42］；乘積平均值（t=3.756）通過了0.05 顯著性檢驗，表明徑流量重建序列包含較多水文信息。

通過對比錫林浩特水文站觀測時段（1972—2016年）徑流量實測值與重建值（圖5），二者的擬合效果較好，在低頻的變化上能較好對應，變化的趨勢總體一致，在個別時段，實測值與重建值基本吻合，進一步證明重建方程的可靠性。

圖5 徑流量重建值與實測值的對比Fig.5 Comparison of the reconstructed values and measured values of runoff

綜上所述，重建方程的各項檢驗參數(shù)均通過了檢驗標準，可用于重建大興安嶺南麓貢格爾河流域歷史徑流量。

2.3 重建徑流量歷史序列的變化特征

貢格爾河1845 年以來徑流量歷史序列存在明顯的枯水和豐水變化（圖6a～b），徑流量序列整體呈下降趨勢，變化頻率較高，M-K 檢驗顯示在1849、1854年和1993年發(fā)生過徑流量突變。近172 a里貢格爾河流域經(jīng)歷了7 個變枯期和6 個變豐期，出現(xiàn)了1853—1855 年、2000—2010 年2 個枯水期和1869—1873 年、1911—1918 年、1952—1958 年3 個豐水期，存在30 個豐水年及18 個枯水年，占比分別為17.4%和10.5%。

圖6 徑流量歷史重建序列變化特征Fig.6 Change characteristics of historical reconstruction runoff series

可以看出，豐水年、枯水年出現(xiàn)年份較集中，1845—1852 年徑流量表現(xiàn)為下降趨勢，變化波動較小（Cv=19.80%），年代際均值較低；1853 年出現(xiàn)首個連續(xù)了3 a 枯水年，1854 年徑流量由枯水期快速向平水期轉變，1856—1866 年豐、枯水年零星出現(xiàn)；1867—1873 年共經(jīng)歷了6 個豐水年，1874 年徑流量由豐水期向平水期轉變，變化波動較劇烈（Cv=34.42%）；20 世紀10 年代豐水年大量出現(xiàn)，研究區(qū)經(jīng)歷持續(xù)時間最長的豐水期，年代際均值達到最高，序列波動較?。–v=23.12%）；1919 年徑流量表現(xiàn)出變枯趨勢，直至1929 年緩慢回升，1934 年研究區(qū)迎來豐水年，年代際均值較上一時期有所提升，20世紀30 年代后期徑流量向平水期轉變；20 世紀50年代研究區(qū)經(jīng)歷最后一個連續(xù)7 a的豐水期后，徑流量表現(xiàn)出顯著的變枯趨勢，變化波動（Cv=16.40%）較上一階段（Cv=26.89%）輕緩，年代際均值明顯下降；在經(jīng)歷了約30 a 的平水期后，徑流量在1993 年發(fā)生突變，變化劇烈（Cv=39.10%）且下降速率較快，最終20世紀末期枯水年集中出現(xiàn)，導致貢格爾河流域歷經(jīng)了持續(xù)時間長達11 a 的枯水期，年代際均值達到近2 個世紀的最低值；21 世紀末期徑流量緩慢回升，2016年徑流量達到豐水年標準。

小波分析發(fā)現(xiàn)貢格爾河徑流量重建序列在95%的置信水平上存在著顯著的3 a、7～12 a、15～22 a及30 a周期（圖7a），小波方差圖（圖7b）顯示15～22 a周期貫穿整個歷史時段，信號最強，表現(xiàn)最顯著。15～22 a 周期可能與太平洋年代際濤動（Pacific decadal oscillation，PDO）有關，亞洲地區(qū)氣候變化與PDO 關系密切［43］，研究區(qū)氣候變化受其影響較大；7～12 a 周期與ENSO 變化周期基本一致，表明ENSO活動對研究區(qū)存在影響。徑流量可以表征大區(qū)域的水文氣候效應，是檢測太陽活動對氣候影響的指標之一［44］，研究區(qū)11 a 周期與太陽黑子活動周期對應；準3 a 周期是我國半干旱地區(qū)重要周期［45］，研究區(qū)水文變化與其密切相關。

圖7 徑流量重建序列周期變化特征Fig.7 Periodic variation characteristics of reconstruction runoff series

2.4 與歷史事件及其他重建結果的比較

貢格爾河發(fā)源于3000 a 前，是研究區(qū)第一大內流河，其豐富的水源為沿岸居民與動植物提供良好的生存條件，是北方極為重要的生態(tài)屏障。為驗證徑流量重建序列的準確性，將重建結果與當?shù)貧v史典籍［24,46］等資料進行對比，發(fā)現(xiàn)重建豐枯時期與歷史氣象災害存在對應。

19 世紀50 年代至70 年代，研究區(qū)降水量較少導致貢格爾河徑流量呈變枯趨勢，連續(xù)經(jīng)歷3 a枯水年；民國元年至九年，赤峰多地接連暴雨，河水暴漲沖毀堤壩，導致山洪暴發(fā)，貢格爾河流域迎來持續(xù)時間最長的豐水期。20 世紀50 年代夏、秋雨澇，數(shù)千公頃農田被毀，徑流量重建序列最后一個豐水期在此時出現(xiàn)。20 世紀60 年代內蒙古降水量時空分布極不均勻，1965 年研究區(qū)遭受建國以來最嚴重的旱災，降水量比常年少50%～90%，連年干旱和上游牧民筑壩引水導致貢格爾河的水量大幅度減少，徑流量序列呈變枯趨勢，1968 年研究區(qū)處于枯水年。1970 年修建五道石門水電站，下游牧場灌溉供水減少，徑流量顯著回升，水電站發(fā)電記錄顯示在1992年發(fā)電量達到歷史最高峰，1992—1993 年貢格爾河徑流量連續(xù)2 a處于豐水年標準之上。然而20 世紀末期全球氣候變暖，大興安嶺冰川融化和連年干旱導致貢格爾河水量驟減，同時河流上游黃崗梁地區(qū)開發(fā)鐵礦，清洗礦砂需要大量河水，工業(yè)廢水對河水造成污染，導致生態(tài)嚴重退化，河水水位持續(xù)下降，貢格爾河歷經(jīng)長達10 a的枯水期；當?shù)卣疄榻鉀Q水面縮減、水位下降等問題，實施了保護治理、生態(tài)補水等一系列修復措施，研究區(qū)環(huán)境得到改善，21 世紀初徑流量開始緩慢回升，在徑流量重建序列終點已達到豐水年標準。

為增強此次重建結果的可信度，將大興安嶺西部利用樟子松重建的伊敏河流域1868—2002 年年徑流量［19］和利用落葉松重建的海拉爾河流域近2個世紀徑流量［20］與本次重建結果的重疊時段進行比較，貢格爾河與伊敏河（r=0.278，N=136，P＜0.01）及海拉爾河徑流量重建序列（r=0.216，N=162，P＜0.01）均呈顯著正相關關系，通過了0.01 顯著性水平檢驗，說明三者關系密切，具有相似的變化趨勢，驗證了本次重建結果的可靠性。

3 結論

（1）紅皮云杉生長過程與徑流量關系密切，受其變化影響顯著，利用年輪寬度標準化年表重建大興安嶺南麓貢格爾河流域近172 a徑流量序列，重建方程準確可靠。

（2）貢格爾河流域1845 年以來徑流量重建序列整體呈下降趨勢，存在明顯的豐枯變化，豐水期出現(xiàn)頻率更高，豐水年占比為17.4%，20世紀末以枯水期為主，但在人為調控下徑流量呈上升趨勢。

（3）小波分析發(fā)現(xiàn)徑流量序列存在顯著的3 a、7～12 a、15～22 a 及30 a 變化周期，15～22 a 周期貫穿整個歷史時段，研究區(qū)水文變化與全球大尺度氣候環(huán)流活動關系密切。