房坤寶 ，崔克儉 ，孫青芳 ，王曙光 ，史雨剛 ，孫黛珍

（1.山西農(nóng)業(yè)大學(xué)農(nóng)學(xué)院，山西太谷030801；2.山西農(nóng)業(yè)大學(xué)文理學(xué)院，山西太谷030801；3.山西農(nóng)業(yè)大學(xué)工學(xué)院，山西太谷030801）

政府間氣候變化專門委員會（Intergovernmental Panel on Climate Change，IPCC）氣候變化 2014 年綜合報告指出，氣候變化會減少亞熱帶地區(qū)可再生的地表水和地下水資源，會加劇行業(yè)間對水的爭奪；在典型濃度路徑（RCP）8.5情景下，許多中緯度地區(qū)和亞熱帶干燥地區(qū)，平均降水可能會減少，極端降水事件很可能強(qiáng)度更大、頻率更高。華北地區(qū)地處北半球中緯度亞熱帶季風(fēng)氣候區(qū)，降水年際變化波動大，空間分布極不均勻，為干旱多發(fā)區(qū)[1-3]。研究表明，過去60 a華北地區(qū)糧食因干旱損失量不斷增加[4]。冬小麥?zhǔn)侨A北地區(qū)主要糧食作物[5]，小麥生長發(fā)育時期恰好是缺水季節(jié)——冬春季[6]。近年來，影響冬小麥生產(chǎn)的冬春季降水量的減少越來越成為制約小麥生產(chǎn)的主要因素，引起政府與社會各界的高度關(guān)注[7]。

20世紀(jì)70年代后期，氣候變暖背景下，全球降水量受東亞季風(fēng)影響逐漸減弱[8]，我國降水年際變化地域差異明顯，西北、東北、長江流域及華南等地降水呈明顯增加趨勢[9-10]，但我國的西北東部、東北大部和華北地區(qū)冬小麥主要生育期降雨明顯減少，其中，西北東部和華北地區(qū)的夏季降水減少趨勢明顯[11]，以華北為中心的冬麥區(qū)冬春季的降水總體呈減少趨勢，加劇了水資源和農(nóng)業(yè)發(fā)展的矛盾并嚴(yán)重威脅到冬小麥的安全生產(chǎn)和糧食增產(chǎn)[12]。

本研究利用多年趨勢線、3 a滑動平均值法、Hurst指數(shù)法、M-K突變分析等降水量變化特征分析方法，分析未來華北地區(qū)北部冬麥區(qū)降水量的變化趨勢。并在此基礎(chǔ)上，基于加權(quán)馬爾可夫鏈預(yù)測方法分析降水量與時間的關(guān)系[13]，構(gòu)建該地區(qū)降水量預(yù)測模型，預(yù)測未來降水量，并驗證降水變化趨勢的準(zhǔn)確性，從而為華北地區(qū)冬小麥生產(chǎn)應(yīng)對未來氣候變化、制定科學(xué)管理措施提供依據(jù)。

1 材料和方法

1.1 數(shù)據(jù)來源

選取華北地區(qū)北京、天津、太原、石家莊4個站點的氣象數(shù)據(jù)來進(jìn)行研究。數(shù)據(jù)來源于中國氣象數(shù)據(jù)網(wǎng)（http://data.cma.cn/）4個站點近60a的氣象資料。

1.2 研究方法

1.2.1 趨勢分析 通過曲線的趨勢線性方程斜率的正負(fù)來判斷變化趨勢，并利用其相關(guān)系數(shù)R2對其進(jìn)行檢驗。首先，對降水量數(shù)據(jù)進(jìn)行處理，得到降水序列，列出散點圖。然后對數(shù)據(jù)進(jìn)行滑動平均，得到新序列。最后，從散點圖中得出趨勢線，并計算趨勢線方程。

1.2.2 Hurst指數(shù)分析 Hurst指數(shù)可以判斷時間序列的持續(xù)性。當(dāng)H＝0.5時，表明時間序列可以用隨機(jī)游走來描述；當(dāng)0.5＜H＜1時，表明時間序列存在長期記憶性，即將來的趨勢與過去一致，H越接近1，持續(xù)性越強(qiáng)。當(dāng)0≤H＜0.5時，表明時間序列具有反持續(xù)性，即將來的趨勢與過去相反，H越接近0，反持續(xù)性越強(qiáng)。本研究采用最常見的R/S（RescaledRangeAnalysis）分析方法計算Hurst指數(shù)。

1.2.3 M-K突變分析法 M-K突變分析法是通過分析UFK和UBK曲線圖，檢測降水序列的變化趨勢，明確降水突變的年份。若UFK或UBK的值大于0，則表明序列呈上升趨勢，小于0則呈下降趨勢。當(dāng)曲線超過臨界直線時，表明該時間區(qū)域上升或下降趨勢顯著。若UFK和UBK這2條曲線出現(xiàn)交叉點，且交叉點在臨界線之間，則交點對應(yīng)的時刻便是突變開始時間。

1.2.4 加權(quán)馬爾可夫鏈降水預(yù)測模型 加權(quán)馬爾可夫鏈預(yù)測模型是國內(nèi)具有代表性的降水預(yù)測方法，是時間和狀態(tài)都離散的馬爾可夫過程。

預(yù)測模型建立如下。

第一，各站點降水序列利用滑動平均法得到新的降水序列。

第二，對新序列進(jìn)行狀態(tài)分配，根據(jù)分級標(biāo)準(zhǔn)確定各時段滑動平均降水量的狀態(tài)。利用樣本均方差對降水量進(jìn)行分級，通常按枯水年、偏枯水年、平水年、偏豐水年、豐水年劃分為5個等級，分別為。其中，為降水平均值，s為標(biāo)準(zhǔn)差。a1在［1.0，1.5］之間取值，a2在［0.3，0.6］之間取值。本研究中a1取1.0，a2取0.5。

第三，通過狀態(tài)分配得到各階到各狀態(tài)的轉(zhuǎn)移矩陣。對第二步所得到的結(jié)果進(jìn)行統(tǒng)計，得到不同滯時（步長1～5）的一步轉(zhuǎn)移概率矩陣。其中，轉(zhuǎn)移概率Pij是描述其演化進(jìn)程最重要的量。

式中，Pij為降水量處于狀態(tài)i時，下一情況變?yōu)闋顟B(tài)j的概率。

第四，確定新序列各階自相關(guān)系數(shù)，并且對自相關(guān)系數(shù)進(jìn)行加權(quán)分析，得到各階自相關(guān)的權(quán)重。其中，各階自相關(guān)系數(shù)計算公式如下。

其中，rk表示第k階（滯時步長為k）自相關(guān)系數(shù)，xt表示第 t年的年降水量，xt＋k表示第 t＋k 年的年降水量，表示降水量序列的多年平均值，n為年降水量序列的長度。

根據(jù)公式（2）確定各滯時權(quán)重系數(shù)pk，其計算公式如下。

其中，m為按預(yù)測需要計算得到的最大階數(shù)。

第五，分別以前年各時段的汛期降水量作為初始狀態(tài)，結(jié)合相應(yīng)的各階轉(zhuǎn)移概率矩陣，可預(yù)測出降水量的狀態(tài)概率。

第六，對同一狀態(tài)的各預(yù)測概率Pi加權(quán)求和，得降水量處于該狀態(tài)的預(yù)測概率。

max{Pi，i∈E}所在的狀態(tài)即是這個預(yù)測時段降水量的預(yù)測狀態(tài)，查詢步驟（2）即可得到降水量的預(yù)測范圍。將預(yù)測值加入到原序列中，重復(fù)以上步驟可進(jìn)行下一年的降水量狀態(tài)的預(yù)測。

第七，預(yù)測結(jié)果計算與驗證。對各狀態(tài)分別賦予相應(yīng)的權(quán)重，權(quán)重集 d＝ ｛d1，d2，d3，d4，d5｝，其中，η為最大概率作用，通常取2或者4，本研究取2；再計算模糊集理論中的級別特征值，H＝Σmii×di，式中，m＝5，i為狀態(tài)，di為權(quán)重。

確定最大概率的狀態(tài)后，根據(jù)公式（5）確定系統(tǒng)在預(yù)報時段的預(yù)報值。

其中，Ti，Bi分別是狀態(tài)區(qū)間值的上限與下限。

最后利用前幾年的預(yù)測值與實際值比較，由于中長期水文預(yù)報的復(fù)雜性，一般認(rèn)為相對誤差小于20%即可。預(yù)測誤差若滿足預(yù)測精度需要，則加權(quán)馬爾可夫鏈模型進(jìn)行降水量中長期預(yù)報有效可行。

2 結(jié)果與分析

2.1 降水量變化趨勢分析

由圖1可知，華北地區(qū)小麥生育期降水趨勢曲線和3 a滑動平均值曲線基本上都呈下降趨勢，下降幅度0.53～7.45 mm/10 a。各站點降水量平均值下降趨勢為4.12 mm/10 a。下降趨勢幅度從高到低排序為：天津＞北京＞太原＞石家莊。同時從3 a滑動平均值曲線可以看出，每個地區(qū)降水量有較強(qiáng)的波動性和規(guī)律性。其中，每個地區(qū)趨勢線的相關(guān)系數(shù)R2的絕對值都遠(yuǎn)小于1，也就說明降水量近年波動較大，下降趨勢不明顯。同時也說明用趨勢線代表降水量變化意義并不大，只能用來判斷變化趨勢。

2.2 突變分析

由圖2可知，北京地區(qū)自1962年以來，UF值一直小于0，降水量呈下降趨勢；在2006—2013年小于下臨界值，其下降趨勢十分顯著。據(jù)UF和UB曲線交叉點的位置，可知北京降水量下降是一個很明顯的突變，具體時間從1967年開始。天津地區(qū)自1978—1979年UF值在0附近波動，其他年份均小于0，降水量呈下降趨勢，但不顯著。太原地區(qū)從1984年至今，UF值小于0但都未超過臨界點，處于降水量下降趨勢，但不顯著；在1953—1983年，UF值大于0，降水量處于上升狀態(tài)，確定太原年降水量從20世紀(jì)80年代開始下降是一個突變現(xiàn)象，具體時間由交叉點可知從1981年開始。另外，幾個交叉點周圍都沒有降水趨勢的變動，所以，不是突變點。石家莊地區(qū)自1967年至今UF值一直小于0，且均未超過臨界值，故降水量一直呈下降狀態(tài)?？傮w來看，降水突變年份可能出現(xiàn)在20世紀(jì)七八十年代， 但近年來各地區(qū)生育期降水量基本均呈下降趨勢。

2.3 Hurst指數(shù)分析

對北方冬麥區(qū)各站點的降水量進(jìn)行R/S分析，得到各站點冬小麥生育期降水量Hurst指數(shù)（表1）。結(jié)果表明，降水量在各個站點的Hurst指數(shù)值均大于0.5且都接近0.5，表明在華北地區(qū)冬小麥生育期的降水量將持續(xù)下降，但其持續(xù)性不強(qiáng)。

表1 華北地區(qū)冬小麥生育期降水量的Hurst指數(shù)值

2.4 加權(quán)馬爾可夫鏈降水預(yù)測模型

為了對未來降水量有個更明確的了解，對各站點近60 a降水量使用馬爾可夫鏈進(jìn)行模擬。利用已知年份降水量驗證降水模型的準(zhǔn)確性，最后得出預(yù)測年份的降水量。以太原地區(qū)為例，過程為：基于各站點降水序列利用滑動平均法得到新的降水序列（表2）；對新序列進(jìn)行狀態(tài)分配，根據(jù)分級標(biāo)準(zhǔn)確定各時段滑動平均降水量的狀態(tài)（表3）；根據(jù)公式（1）通過狀態(tài)分配得到各階到各狀態(tài)的轉(zhuǎn)移矩陣；確定新序列各階自相關(guān)系數(shù)，并且對自相關(guān)系數(shù)進(jìn)行加權(quán)分析（由公式（2）計算得到各階自相關(guān)系數(shù)，k1＝0.332，k2＝0.108，k3＝0.333，k4＝0.171，k5＝0.054）；根據(jù)狀態(tài)轉(zhuǎn)移矩陣與權(quán)重系數(shù)計算加權(quán)求和值（表5）。

利用降水序列中的狀態(tài)，計算近60 a的滑動平均降水量的各種步長的狀態(tài)轉(zhuǎn)移概率矩陣如下。

由公式（3），（4）計算得到表 4。由表 4 可知，當(dāng)i＝5時，P5＝0.412，為最大值，說明 2014年的滑動平均降水量狀態(tài)為5，即數(shù)值區(qū)間為［228.5，＋∞）。

由公式（5）求出2013年滑動平均降水量為282.38 mm，而實測值為340.40 mm，相對誤差為17.04%。

同理，重復(fù)以上方法可得各站點其他時段的滑動平均降水量預(yù)測驗證（表5）。

從表5可以看出，2013—2017年這5 a的馬爾可夫鏈預(yù)測值與真實值誤差都在50 mm上下，預(yù)測狀態(tài)基本上相同，與實際情況基本吻合。從數(shù)值上講，太原、北京、天津的預(yù)測相對誤差都較小，都不超過20%。

由于中長期水文預(yù)報的復(fù)雜性，一般認(rèn)為相對誤差小于20%即可。預(yù)測誤差滿足預(yù)測精度需要，證明加權(quán)馬爾可夫鏈模型進(jìn)行降水量中長期預(yù)報是有效可行的。因此，加權(quán)馬爾可夫鏈模型可以進(jìn)行下一步預(yù)測。

重復(fù)上述方法，可得各站點未來5 a的預(yù)測降水量（表 6）。

表2 太原市小麥生育期降水序列及狀態(tài)

表3 太原滑動平均降水量等級劃分

由預(yù)測量可知，未來降水量基本上都在一個狀 態(tài)之間，存在下降趨勢，但是下降趨勢不明顯。

表4 2014年太原滑動平均降水量預(yù)測

表5 冬小麥生育期滑動平均降水量預(yù)測驗證

表6 降水預(yù)測

3 結(jié)論與討論

本研究利用華北冬麥區(qū)各站點近60 a的逐年降水?dāng)?shù)據(jù)資料，綜合運用了線性趨勢法、3點滑動平均法、Hurst指數(shù)分析、M-K非參數(shù)突變檢驗和加權(quán)馬爾可夫鏈預(yù)測模型，對華北地區(qū)各站點及站點平均降水量的變化規(guī)律進(jìn)行了研究與預(yù)測。

線性趨勢表明，近年來，各站點小麥生育期降水量均略呈下降趨勢，下降0.53～7.54 mm/10 a，與華北地區(qū)整體變化趨勢表現(xiàn)一致[14]，這與張皓等[15]的研究結(jié)果一致。然而，孫燕等[16-17]研究認(rèn)為，華北地區(qū)秋冬春季降水量無明顯下降趨勢。這可能是因為研究站點單一以及研究區(qū)域不同而引起的。

M-K非參數(shù)統(tǒng)計檢驗結(jié)果表明，近60 a來，太原、北京的小麥生育期降水量的突變可能發(fā)生在1981，1967年，與張磊[18]研究整個北方地區(qū)的降水突變年份1972年有一定的誤差，這是由于研究區(qū)域范圍不同引起的，同時由于小麥生育期降水量數(shù)據(jù)時間序列有限，以及未對降水量數(shù)據(jù)進(jìn)行拆分，使天津、石家莊等部分地區(qū)的突變年份未明確，還有待更長時間序列上的驗證。

華北地區(qū)小麥生育期各站點降水量的Hurst指數(shù)均大于0.5，表明未來下降趨勢將繼續(xù)持續(xù)，但由于各站點指數(shù)都接近0.5，說明持續(xù)性不強(qiáng)。

由加權(quán)馬爾可夫鏈預(yù)測結(jié)果可知，該方法誤差在50 mm左右，相對誤差不超過20%，與其他學(xué)者結(jié)果大致相同[19-20]，適用于本地區(qū)冬小麥生育期降水量的預(yù)測。未來冬小麥生育期降水量基本上都在一個狀態(tài)區(qū)間，大部分都在平水年間，比較預(yù)測年份之前的年降水量有小幅度的下降，以此驗證了下降趨勢但持續(xù)性不強(qiáng)的結(jié)論。