門寶輝，張 騰

（華北電力大學(xué) 水利與水電工程學(xué)院，北京 102206）

1 研究背景

全球氣候變暖和高頻人類活動加速了全球及區(qū)域尺度的水文循環(huán)過程，并導(dǎo)致降水量時空變化及其區(qū)域差異的不確定性增大［1-3］，而降水量時空分配不均勻除了直接體現(xiàn)在極端天氣現(xiàn)象外，還對農(nóng)業(yè)生產(chǎn)、城市洪澇、經(jīng)濟(jì)社會發(fā)展等多個方面產(chǎn)生重要影響［4-5］。

北京地區(qū)屬溫帶季風(fēng)氣候，常年干旱少雨，水資源人均占有量低于300 m3，是我國嚴(yán)重水資源短缺的地區(qū)之一。全球溫室效應(yīng)和人口快速增長等問題加劇了該地區(qū)綠水資源匱乏問題，水資源已成為掣肘北京地區(qū)經(jīng)濟(jì)高質(zhì)量發(fā)展和社會生態(tài)文明環(huán)境建設(shè)的關(guān)鍵因素，因此明確該區(qū)域降水量的未來發(fā)展態(tài)勢具有深遠(yuǎn)意義?，F(xiàn)階段我國學(xué)者對包含北京在內(nèi)的多個地區(qū)的降水序列變化特征等方面開展了許多研究。

李鵬程等［6］采用M-K檢驗法和小波變換對北京、天津地區(qū)1958—2008年降水序列變化特征及其與汛期降水周期關(guān)聯(lián)進(jìn)行分析，結(jié)果表明研究時域內(nèi)北京、天津地區(qū)的年降水量與汛期降水量均呈現(xiàn)出減少趨勢；宋曉猛等［7］分析了北京地區(qū)1980—2012年降水歷時和降水等級的時空變化特征；尤煥苓等［8］對1981—2010年間北京地區(qū)20個常規(guī)氣象觀測站日降水?dāng)?shù)據(jù)進(jìn)行統(tǒng)計計算，發(fā)現(xiàn)從1980年代以來北京地區(qū)極端類型降水量呈現(xiàn)出減少趨勢；朱龍騰等［9］借助經(jīng)驗?zāi)B(tài)分解北京地區(qū)1951—2009年共59年降水序列，分解結(jié)果表明該地區(qū)降水量在短期內(nèi)呈波動性減少趨勢；魯向暉等［10］采用累積距平法和Morlet小波變換等方法對贛江流域1966—2015年不同時間尺度下的降水時空分布特征進(jìn)行分析，結(jié)果表明該流域在年尺度上汛期降水量緩慢增加，在其他時間尺度下呈不顯著減少勢態(tài)；徐東坡等［11］對全國1956—2018年降水序列進(jìn)行Pettitt檢驗和小波分析，發(fā)現(xiàn)全國降水量呈現(xiàn)非顯著性增長趨勢，并且降水量增加趨勢從西北向東南減弱。

以上研究工作剖析了北京及全國各地降水序列在不同時空尺度下的變化規(guī)律和演變趨勢，具有良好的借鑒經(jīng)驗，但其降水?dāng)?shù)據(jù)多截止在2015年以前并且序列長度局限在50～70 a左右，其對氣候變化規(guī)律研究和未來降水演變趨勢不具有較強(qiáng)說服力。鑒此，本文研究過程如圖1所示，以北京地區(qū)1724—2019年近三百年超長降水序列深入研究該地區(qū)降水變化特征和未來演變趨勢，首先使用BS-Pettitt耦合模型［12］對長降水序列進(jìn)行多突變點檢測，然后采用極點對稱模態(tài)分解算法［13］分析多時間尺度周期性以及趨勢性變化規(guī)律，并選取多段代表性短降水時序進(jìn)行對比分析；最后基于隨機(jī)森林算法［14］構(gòu)建不同長度時序下的北京地區(qū)年降水序列隨機(jī)模擬模型。在降水序列組成成分分析的基礎(chǔ)上引入隨機(jī)模擬，增強(qiáng)研究完整性，有助于更好地研究水文時序變化情況、對比分析揭示長時間尺度序列演變規(guī)律，以期為該區(qū)域水資源系統(tǒng)規(guī)劃設(shè)計和提高城市適應(yīng)氣候變化能力提供支撐。

圖1 技術(shù)路線

2 方法與數(shù)據(jù)

2.1 BS-Pettitt突變檢驗法Pettitt檢驗法是由Pettitt［15］新發(fā)展的一種非參數(shù)檢測法，其在有效診斷水文時序突變情況的同時，還可檢驗突變是否具備數(shù)學(xué)統(tǒng)計意義上的顯著性，擬定統(tǒng)計量Ut，N：

式中：xt，xi為序列中的隨機(jī)變量；N為序列總長度；Ut，N為新序列，由逐次統(tǒng)計序列中前樣本數(shù)值大于后樣本數(shù)值構(gòu)成。若有t時刻，滿足：

式中：kt，N為最大統(tǒng)計量；p為顯著性檢驗指標(biāo)。即可判定該點為此序列的突變點，并且其在數(shù)學(xué)統(tǒng)計意義上表現(xiàn)顯著［16］。Pettitt檢驗在單突變點檢測方面具備高效率和高精準(zhǔn)度，但其用于多突變點檢測境況時表現(xiàn)欠佳。因此利用二元分割法（Binary Segmentation，BS）的迭代功能，通過將水文序列分段進(jìn)行Pettitt檢驗達(dá)到多突變點識別的目的，即構(gòu)建BS-Pettitt耦合模型［12］。耦合模型具體步驟如下：

（1）設(shè)水文序列Xt（t＝1，2，…，N），采用Pettitt法對序列X［t1∶t2］（t1＝1，t2＝296）進(jìn)行突變檢測，并將該突變點記為k，并進(jìn)行下一步，否則認(rèn)為該序列無突變點。

（2）以上一步突變點k為分段點，將水文序列分成左右兩部分分別重新進(jìn)行Pettitt突變檢測。將左半段中突變點位置記為kl，并重復(fù)檢測，直到無突變點檢出，記錄末次分段序列t2，令kup＝t2。

（3）與（2）相似，將右半段中突變點位置記為kr，并重復(fù)檢測，直到無突變點檢出，記錄末次分段序列t1，令kdown＝t1-1。

（4）若kup＝kdown，表明原始水文序列只有一個突變點k，否則，重新構(gòu)建序列X［t1∶t2］（t1＝kup＋1，t2＝kdown），并重復(fù)（1）—（3），直到檢測出水文序列全部突變點，即突變點集合Mi（i＝1，2，…，m）。

（5）將突變點集合排序后，添加1和N重新構(gòu)建序列Sj（j＝1，2，…，m＋2），并依次對X［t1∶t2］（t1＝Sj，t2＝kdown）進(jìn)行Pettitt突變檢測，獲得新的突變點集合，并重復(fù)檢測，直到突變點數(shù)目恒定時，此刻的突變點集合即為多突變點最終檢測結(jié)果。

2.2 去趨勢預(yù)置白算法去趨勢預(yù)置白算法（Trend-free Pre-whiting，TFPW）是2002年Yue等［17］提出的一種預(yù)處理技術(shù)，該方法能高效剔除序列分析時的自相關(guān)性干擾，目前已在多領(lǐng)域得到廣泛使用［18-20］。本次研究中，TFPW 主要用于多點突變檢驗之前，借助剔除趨勢和預(yù)置白化兩部分充分削減原始序列自相關(guān)性對多突變點檢測的干擾。具體步驟如下：

（1）計算水文序列為Xt（t＝1，2，…，N）的線性趨勢β，則有：

（2）去除水文序列趨勢成分后形成新序列Yt：

（3）計算新序列Yt的一階自相關(guān)系數(shù)r，并對自相關(guān)系數(shù)進(jìn)行顯著性檢驗（顯著性水平取0.1），如果通過顯著性檢驗，則直接進(jìn)行BS-Pettitt突變檢驗，若未通過，進(jìn)行第（4）步的預(yù)處理。

（4）去除序列自相關(guān)項形成新序列Y′t；重新補(bǔ)加趨勢成分形成不受自相關(guān)性干擾的新序列Y″t：

2.3 極點對稱模態(tài)分解算法自1998年由Huang等［21］首次明確提出經(jīng)驗?zāi)B(tài)分解（Empirical Mode Decomposition，EMD），其改進(jìn)算法例如集合經(jīng)驗?zāi)B(tài)分解（Ensemble Empirical Mode Decomposition，EEMD）［22］、局部均值分解（Local Mean Decomposition，LMD）［23］、極點對稱模態(tài)分解（Extreme-point Symmetric Mode Decomposition，ESMD）等在信號診斷［24］、時序分析［25］等方面均得到廣泛使用。

ESMD算法是王金良等［13］基于EMD算法的新發(fā)展，其實質(zhì)是以內(nèi)部極點對稱插值法取代EMD算法外包絡(luò)線的三次樣條插值法對原始研究序列進(jìn)行插值處理，隨后優(yōu)化最末殘余分量為 “最小二乘”意義下的最優(yōu) “自適應(yīng)性全局平均線”，并憑此提供分解過程中的最優(yōu)篩查次數(shù)，該算法在高效抑制EMD處理過程中 “模態(tài)混疊”等現(xiàn)象時表現(xiàn)良好。ESMD算法分解過程如下：

（1）標(biāo)識研究序列X所有極值點，記為Ei（i＝1，2，3，…，n）。

（2）使用直線連接緊鄰極值點，標(biāo)記其線段中點為Zi（i＝1，2，3，…，n-1）。

（3）采用線性插值法延展序列首尾中點Z0、Zn。

（4）利用上述步驟中產(chǎn)生的n＋1個中點構(gòu)建m條內(nèi)插曲線Li（i＝1，2，3，…，m）（m≥1），并計算上述曲線的均值曲線：

（5）對X-L*新序列循環(huán)上述（2）—（4）步，截至符合以下任一條件：（a）均值曲線（ε為允許誤差）；（b）迭代循環(huán)次數(shù)達(dá)到設(shè)定的允許上限值F。此時成功分解第一個本征模函數(shù)（Intrinsic Mode Function，IMF）分量M1。

（6）對X-M1新序列循環(huán)上述（2）—（5）步，依序分解得到Mi（i＝1，2，3，…，N），直至殘余分量R滿足其為單調(diào)序列或小于等于設(shè)定的剩余極值點數(shù)，此時通過分解得到了所有IMF分量以及殘差分量。

（7）讓篩選次數(shù)最大值限定在有效區(qū)間［Fmin，F(xiàn)max］內(nèi)取值，之后重復(fù)以上所有過程。然后計算不同篩選次數(shù)最大值下對應(yīng)的X-R序列方差σ2和方差比率σ/σ0（σ0為原始信號序列的標(biāo)準(zhǔn)差），找到方差比率最小時的篩選最大值F0。更改迭代次數(shù)允許最大值為F0，重復(fù)上述全部步驟，便可分解得到殘差分量，亦是 “自適應(yīng)性全局平均線”。

至此經(jīng)由ESMD算法分解獲得全部IMF分量以及殘差分量，故原始序列可表示為：

2.4 隨機(jī)森林算法Breiman基于Bagging算法改進(jìn)回歸樹模型，創(chuàng)新地提出了隨機(jī)森林算法［14］，該算法主要由子訓(xùn)練樣本集和子回歸模型（決策樹）組成，其通過Bootstrap重抽樣從原樣本集中有放回地抽取多個樣本組成樣本容量與原樣本集合一致的子訓(xùn)練樣本集，并對每個子訓(xùn)練樣本集構(gòu)建子回歸模型，即組成隨機(jī)森林模型，綜合所有子回歸模型的子模擬結(jié)果得到最終模擬結(jié)果。該算法在大量的理論和實例中均表現(xiàn)出較高模擬精準(zhǔn)度，目前已得到廣泛應(yīng)用［26-28］，模型計算流程見圖2。

圖2 隨機(jī)森林模型計算流程圖

在隨機(jī)森林模型構(gòu)建后，選定均方根誤差RRMSE、平均相對誤差絕對值MMRE和決定系數(shù)R2三個指標(biāo)來評估模型的模擬精度，當(dāng)前兩者越小、后者越大時模型模擬精度越好。各指標(biāo)計算公式分別為：

式中：xi為模型擬合值，mm；yi為年降水實測值，mm；ˉx為模型擬合系列值的均值，mm；ˉy為年降水實測系列值的均值，mm。

2.5 數(shù)據(jù)資料本次研究采用的數(shù)據(jù)資料是北京地區(qū)長時序面降水序列，序列長度為296年（1724—2019）（見圖3），前250年（1724—1973）數(shù)據(jù)來自1970年代我國氣象學(xué)者依據(jù)古籍 《晴雨錄》記載數(shù)據(jù)并關(guān)聯(lián)有儀器觀測以來的北京氣象站實測降水?dāng)?shù)據(jù)后推算延伸并整編的 《北京250年降水（1724—1973）》①中央氣象局研究所.北京250年降水（1724—1973）.1975.，故認(rèn)為此匯編數(shù)據(jù)具有一致性，另外從北京氣象局獲取北京氣象站后續(xù)年份降水?dāng)?shù)據(jù)，依托北京氣象站得到296年北京地區(qū)1724—2019年長降水序列數(shù)據(jù)，本研究主要集中于降水序列成分分析和隨機(jī)模擬，旨在探明降水波動變化規(guī)律，而非精確量值變化，故該時序數(shù)據(jù)能夠滿足本次研究的需要，具有一定的空間一致性和可靠性。

圖3 北京地區(qū)1724—2019年降水序列

3 降水序列組成成分分析

3.1 突變檢驗首先采用TFPW 檢驗北京地區(qū)296年降水序列是否存在自相關(guān)并進(jìn)行預(yù)處理，隨后將序列代入BS-Pettitt模型進(jìn)行突變檢驗，檢驗結(jié)果初步表明北京地區(qū)296年降水序列存在7個突變點，突變點發(fā)生的年份分別是1742年、1770年、1813年、1871年、1893年、1947年和1999年。各時段內(nèi)各突變點的Pettitt統(tǒng)計量Ut，N變化過程如圖4（a）—（g）所示，根據(jù)突變檢驗的結(jié)果得到北京地區(qū)296年降水量序列的均值變化過程（見圖4（h）），可以看出降水序列分段均值變化與降水序列變化情況匹配良好，能較為完整地體現(xiàn)出降水序列變化趨勢，較為充分表明BS-Pettitt耦合模型具有較高準(zhǔn)確性，并且結(jié)果與以往研究結(jié)論大致相同［29］。

圖4 多點突變檢驗結(jié)果

3.2 周期項與趨勢項分析采用ESMD算法對北京地區(qū)296年降水序列進(jìn)行分析，分解得到7個IMF分量及1個殘差分量，IMF分量分解結(jié)果如圖5所示，從圖中可以發(fā)現(xiàn)，所有IMF分量的振蕩變化均線都在零值附近，其局部變化也基本對稱，并且各分量的波動都存在一定的周期性。通過初步觀察分量間幅值、頻率的相對變化，可以發(fā)現(xiàn)隨著分解尺度的增大，其振幅隨之減小，而波長逐步增長。殘差分量沒有表現(xiàn)出明顯的周期變化規(guī)律，但由于研究序列長度有限，故不排除其屬于更長周期變化的組成成分的可能性。原序列殘差分量可以表征北京地區(qū)年降水序列在整個研究時域上呈現(xiàn) “上升—平穩(wěn)—下降”的變化趨勢。

圖5 ESMD算法IMF分量結(jié)果

3.3 結(jié)果分析與討論水文時序突變即意味著時序某點前后統(tǒng)計規(guī)律發(fā)生顯著變化，能否準(zhǔn)確并全面診斷水文時序突變點與水文序列樣本容量、突變點數(shù)目、突變位置、樣本及總體分布等因素息息相關(guān)，針對北京地區(qū)降水序列突變點檢驗研究數(shù)據(jù)多局限于50～70 a的樣本容量小的現(xiàn)狀，依據(jù)BSPettitt模型在多點突變檢驗過程中選取含有顯著突變點的連續(xù)短時序降水時段1724—1769年、1966—2019年、1776—2000年、1813—2000年、1871—2000年、1896—2000年，并結(jié)合研究區(qū)現(xiàn)有文獻(xiàn)，選取以往研究時序1951—2015年、1961—2010年、1958—2008年、1980—2012年，采用K-S（Kolmogorov-Smirnov）研判初步檢驗結(jié)果中各不同尺度短時序是否與296 a長時序降水具有顯著性差異。設(shè)定原假設(shè)H0：兩個長短時序數(shù)據(jù)所服從分布不一致（或存在顯著性差異），檢驗結(jié)果見表1。

表1 不同時段降水時序K-S檢驗結(jié)果（顯著性水平α＝0.05）

結(jié)合K-S檢驗結(jié)果與現(xiàn)有短時序突變點研究［6］結(jié)果發(fā)現(xiàn)，雖然短時序檢測出的突變點如1961年、1995年，仍含在1724—2019年長時序突變檢出點中且具有統(tǒng)計顯著性，但與原時序相比仍存在不小誤差，時序長短是直接影響水文時序演變規(guī)律精準(zhǔn)度因素之一，在診斷降水序列突變點的位置及分布規(guī)律時，長時間序列更能探明更多的突變點和更長更精準(zhǔn)的突變分布規(guī)律，更具有優(yōu)勢和診斷精確性。另一方面，K-S檢驗結(jié)果同時也對BS-Pettitt初步檢驗結(jié)果進(jìn)行反饋研判，1724—1769年序列與原序列存在顯著性差異，但仍檢測出其內(nèi)含顯著性突變點，究其原因是因為Pettitt檢驗自身算法不適宜在序列首尾進(jìn)行檢測，存在誤差，故舍去突變點1742年，分析得到原序列存在1770年、1813年、1871年、1893年、1947年和1999年6個突變點。根據(jù)K-S檢驗結(jié)果中顯著性較高且兼具代表性的短降水時序1813—2000年、1896—2000年、1951—2015年、1961—2010年，同樣對其進(jìn)行ESMD分解，分解結(jié)果見圖6。

由圖5、圖6中不同時序長度降水序列ESMD算法分解結(jié)果可以看出，各IMF分量均存在明顯的周期性波動，其代表著水文序列中固有的不同時間尺度下的周期變化規(guī)律，并且隨著階數(shù)的增加，IMF分量振蕩幅度減小，波長逐步增大，反映出更大時間尺度的周期變化特征。對比分析不同時序降水序列結(jié)果，依次分解出7、6、5、4、4個IMF分量和1個殘差余項，結(jié)果表明長時間序列能分解出更多的IMF分量且分解結(jié)果更平穩(wěn)、周期波動更單一，利用周期圖法處理各IMF分量，并將各模態(tài)分量對時序整體數(shù)據(jù)變化影響以方差貢獻(xiàn)率呈現(xiàn)，各降水時序的IMF分量平均周期及方差貢獻(xiàn)率見表2。

圖6 短時序ESMD分解結(jié)果

從表2中可以看出，首先，北京地區(qū)降水周期規(guī)律存在2.5～4 a、7～15 a、25～35 a、74 a左右、95～100 a的周期變化尺度，并且1724—2019年序列分解結(jié)果中涵蓋其余4個時序的周期結(jié)果，從分解結(jié)果來看，更長時間序列通過ESMD算法能分解出降水序列隱含的更多更完整的周期尺度。其中，僅有1724—2019年、1813—2000年序列檢測出90 a尺度周期。其次，根據(jù)方差貢獻(xiàn)率計算結(jié)果可知，北京年降水序列在年際變化尺度以2.5～4.5 a的周期振蕩變化為主，在年代際尺度上，則以16～20 a的周期振蕩變化為主。1724—2019年長序列能顯示北京地區(qū)三百年時間跨度上更加完整的周期變化尺度，其中以2.5～4 a的周期振蕩最明顯，其方差貢獻(xiàn)率最大，年降水量以 “下降—上升—下降”的交替變化貫穿始終；年代際尺度中7～15 a和25～35 a的振蕩較為明顯，且呈現(xiàn)疏密相間的波動特征，其在1760—1790年和1835—1865年振蕩波動平穩(wěn)且特征相似，說明降水在此周期尺度下沒有明顯波動；74 a和95～100 a的周期變化尺度中，降水周期波動明顯減緩，降水整體呈現(xiàn)出 “下降—上升—下降—上升—下降—上升—下降”的循環(huán)變化規(guī)律，很顯然這是短時序無法檢測出的周期尺度特征，長時序ESMD分解在降水規(guī)律分析中體現(xiàn)出獨特的優(yōu)勢。時序長短不僅影響IMF分解數(shù)量和周期尺度分析，在降水序列趨勢分析上同樣起到關(guān)鍵作用。各時序殘差余項結(jié)果對比分析見圖7。

圖7 各時序ESMD分解殘差余項對比分析

表2 各時序ESMD周期分解結(jié)果對比

由圖7可見，所有序列的殘差余項均在近代后期開始呈現(xiàn)下降趨勢，但時序長短嚴(yán)重制約了降水序列趨勢項規(guī)律分析，通過旱澇等級和樹木年輪記錄等代用降水資料重建的降水時序表明我國東部降水序列存在1600—1910年上升，1900年之后下降的變化趨勢［30］，檢測序列中僅有1724—2019年、1813—2000年序列殘差分量可以表征北京地區(qū)年降水序列在整個研究時域上呈現(xiàn) “上升—平穩(wěn)—下降”的變化趨勢，其余短時序僅能代表當(dāng)前研究時域內(nèi)部分趨勢變化，均難以充分揭露序列演變規(guī)律，其中1961—2010年序列趨勢變化結(jié)果誤差明顯，究其原因可能與ESMD算法本身精準(zhǔn)度存在關(guān)聯(lián)。

在北京地區(qū)降水序列成分分析中，1724—2019年長時間序列能檢測出更多的突變點、更多更精準(zhǔn)的周期尺度和更完整的趨勢變化規(guī)律，表明長時間序列在降水序列成分分析體現(xiàn)出較強(qiáng)精準(zhǔn)度和完整性，更有利于探明該地區(qū)長期降水變化規(guī)律和演變趨勢，為后期降水模擬預(yù)測和防洪抗?jié)程峁┲匾碚撝巍?/p>

4 降水序列隨機(jī)模擬

4.1 基于隨機(jī)森林的隨機(jī)模擬模型構(gòu)建依據(jù)前文ESMD算法對北京地區(qū)年降水序列的周期項和趨勢項分析結(jié)果，初步掌握北京地區(qū)現(xiàn)狀和未來一段時期的降水量變化趨勢。本節(jié)通過隨機(jī)森林算法構(gòu)建北京地區(qū)年降水模擬模型。

降水模擬模型構(gòu)建可分為擬合因子選取和模型參數(shù)優(yōu)化兩部分內(nèi)容，本次模型構(gòu)建依次選用時序當(dāng)年四季降水量和前10年降水量共14個解釋變量作為備選擬合因子，構(gòu)建隨機(jī)森林模型對備選擬合因子進(jìn)行重要性估計。本次研究序列1724—2019年，為進(jìn)一步驗證長時序水文序列在時域分析中的優(yōu)勢，選取K-S檢驗中顯著性水平較高且具有代表性的時序1813—2019年和1951—2019年，考慮模型訓(xùn)練期和驗證期合適比例，分別將1724—1993年、1813—1993年、1951—1993年作為訓(xùn)練期樣本，其余作為驗證期樣本，三段時序擬合因子重要性估計結(jié)果見表3。

表3 備選擬合因子重要性排序

選取表3中重要性最高的5個解釋變量作為隨機(jī)森林模型擬合因子。隨機(jī)森林模型主要待定參數(shù)有兩個：子預(yù)報決策樹模型數(shù)Nt、決策樹節(jié)點劃分待選變量數(shù)M，前者愈大，愈能有效抑制隨機(jī)模擬過程中的過擬合問題；后者愈大，愈能縮小子回歸模型間的差異性。Nt通常取值應(yīng)較大，M取值約等于擬合因子總數(shù)的1/3［31］。另結(jié)合數(shù)據(jù)實際情況將本次模型構(gòu)建Nt取值2500，M取值3。

4.2 隨機(jī)模擬結(jié)果以北京地區(qū)不同時序年降水量和相應(yīng)的擬合因子系列作為模型的輸入，進(jìn)一步分別構(gòu)建北京地區(qū)降水模擬的隨機(jī)森林模型。各北京地區(qū)降水模擬模型訓(xùn)練期和驗證期模擬分析結(jié)果見圖8。

圖8 不同時序降水隨機(jī)模擬模型擬合結(jié)果分析

基于訓(xùn)練期和驗證期模擬及相應(yīng)年降水序列，選用決定系數(shù)R2、均方根誤差RRMSE和平均相對誤差絕對值MMRE三個指標(biāo)評定模型精準(zhǔn)度，指標(biāo)評定結(jié)果見表4。

表4 各時序年降水序列模擬模型訓(xùn)練期及驗證期模擬性能對比分析

由表4可知，由三個時序分別構(gòu)建的北京地區(qū)降水模擬模型在訓(xùn)練期均表現(xiàn)出較高水平的模擬精度，其中1724—2019年序列誤差最小，1951—2019年序列誤差最大，三個隨機(jī)模擬模型確定性系數(shù)R2均達(dá)0.98，平均相對誤差絕對值MMRE均小于5%，除1951—2019年模型外，其余時序模型均方根誤差RRMSE均小于25。驗證期精度雖低于訓(xùn)練期，但三個隨機(jī)森林模擬模型平均相對誤差絕對值小于15%，確定性系數(shù)均達(dá)0.80以上，其中1724—2019年序列誤差最小，1813—2019年序列在均方根誤差和平均相對誤差絕對值指標(biāo)均優(yōu)于1951—2019年序列。結(jié)合圖8可知，本次構(gòu)建的所有年降水隨機(jī)森林模型模擬精度均較高，在訓(xùn)練期和驗證期均表現(xiàn)出較高的模擬精準(zhǔn)度。

4.3 結(jié)果討論本次研究從定性和定量兩個視角實現(xiàn)北京地區(qū)年降水模擬，定性為定量提供依據(jù)，定量具體化表述定性，兩者結(jié)果相輔相成。在定性分析中，ESMD算法分解結(jié)果中各IMF分量表明北京地區(qū)近三百年降水序列存在2.5～4年、7～15年、25～35年、74年左右、95～100年的周期尺度變化規(guī)律，其中以2.5～4年的周期年際變化尺度為主；殘差余項表明降水序列存在 “上升—平穩(wěn)—下降”的整體變化趨勢，且近年來下降勢態(tài)加劇。

基于隨機(jī)森林模型構(gòu)建的北京地區(qū)年降水隨機(jī)模擬模型高精度地實現(xiàn)了北京地區(qū)年降水序列的隨機(jī)模擬，由1724—2019年原序列構(gòu)建的隨機(jī)森林模型在訓(xùn)練期和驗證期均誤差最小，1813—2019年序列在訓(xùn)練期和驗證期表現(xiàn)同樣優(yōu)于1951—2019年序列，表明長時間序列在提高降水序列隨機(jī)模擬精準(zhǔn)度方面同樣起到關(guān)鍵作用，但在模型驗證期依然存在一定誤差，需在后續(xù)工作進(jìn)一步豐富擬合因子涉及面和時序尺度，為北京地區(qū)長時期年降水序列演變趨勢研究提供科學(xué)依據(jù)。

5 結(jié)論

本次研究通過TFPW-BS-Pettitt法和ESMD算法對北京地區(qū)1724—2019年等不同降水序列組成成分進(jìn)行分析，并利用隨機(jī)森林模型對年降水序列進(jìn)行隨機(jī)模擬，得出結(jié)論如下：

（1）北京地區(qū)296年長降水序列存在 “上升—平穩(wěn)—下降”的趨勢變化特征和2.5～4年、7～15年、25～35年、74年左右和95～100年的準(zhǔn)周期波動規(guī)律，并且在整個研究時域內(nèi)存在1770年、1813年、1871年、1893年、1947年和1999年6個突變點。

（2）構(gòu)建隨機(jī)森林模型對北京地區(qū)年降水量進(jìn)行隨機(jī)模擬，其中，模型訓(xùn)練期和驗證期平均相對誤差絕對值分別小于5%和15%，結(jié)果表明隨機(jī)森林模型擁有較好的模擬效果。隨機(jī)森林模型作為長期降水序列隨機(jī)模擬的新途徑表現(xiàn)良好，但由于擬合因子是模型構(gòu)建的關(guān)鍵所在，因此如何選取和篩查科學(xué)、合理的擬合因子有待于在后續(xù)研究中不斷修改和完善。

（3）通過對比分析可知，長時間序列在降水序列成分分析和隨機(jī)模擬均表現(xiàn)更好，既能檢出降水序列隱含的多突變點、更精準(zhǔn)更完整的周期尺度和趨勢變化規(guī)律，更能有效提高隨機(jī)模擬精度，增強(qiáng)研究邏輯性和說服力。

（4）在傳統(tǒng)水文序列成分分析的基礎(chǔ)上引入隨機(jī)模擬，將研究時域從歷史時期延伸到未來短期內(nèi)，增強(qiáng)研究完整性，有助于更好地研究水文時序變化情況、揭示序列演變規(guī)律。

致謝：感謝水利部發(fā)展研究中心李淼博士在數(shù)據(jù)收集方面給予的支持和幫助！