王秀杰 ，齊喜玲 ，滕振敏，李丹丹

(1. 天津大學 水利工程仿真與安全國家重點實驗室，天津 300072； 2. 天津大學 建筑工程學院，天津 300072； 3. 崇左市左江治旱工程管理中心，廣西 崇左 532200)

洪災具有破壞性強、影響范圍廣、突發(fā)性強、頻率高等特點[1]，極易造成社會經濟和人民生命財產損失。據(jù)統(tǒng)計，洪澇災害損失在全球所有自然災害損失中占比高達40%[2]。因此，及時準確的洪水預報作為非工程措施[3]，在防洪減災中發(fā)揮著重要作用[4]。洪水預報中，氣象數(shù)據(jù)的可靠性是關鍵[5]，然而降雨資料短缺一直是個難題[6]。因此，插補移用降雨資料是非常必要的[7]。常用的方法有線性相關法、反距離權重法、泰森多邊形法及克里金插值法等[8-9]。國內外諸多學者在降雨資料空間插補方面做了大量研究[10-11]，但較少考慮降雨序列的非平穩(wěn)性。隨著全球氣候變暖，極端降雨事件的頻次和強度均發(fā)生了變化，變化環(huán)境導致極端降雨存在非平穩(wěn)性[12-13]。臺風雨的極端降雨強度大、分布范圍廣，一定影響范圍內的雨量站間的降雨資料存在線性協(xié)整關系。因此，考慮降雨數(shù)據(jù)的線性協(xié)整關系并建立模型插補短缺數(shù)據(jù)，是提高降雨數(shù)據(jù)插補質量的一種新思路。

因多種因素干擾，降雨序列往往是缺失的[14]。許多學者對無資料地區(qū)數(shù)據(jù)移用進行了研究，通常采用鄰近測站資料替代研究區(qū)域數(shù)據(jù)的方法[8]，其中，金倩芳[15]分別以衛(wèi)星降雨數(shù)據(jù)及測站降雨數(shù)據(jù)作為驅動建立水文模型，計算結果表明衛(wèi)星降雨數(shù)據(jù)驅動的水文預報模型在洪峰預報精度上低于實測降雨數(shù)據(jù)驅動的模型?？梢姳O(jiān)測儀器精度不高會降低衛(wèi)星雷達降雨產品的空間分辨率[16]，造成數(shù)據(jù)精度不高，使得衛(wèi)星降雨數(shù)據(jù)驅動的水文預報模型在流域的數(shù)據(jù)移用方面適用性較差[17]。彭安幫等[18]對流域劃分水文相似組并建立子流域洪水預報模型，對無資料地區(qū)進行參數(shù)移植，模擬中發(fā)現(xiàn)在洪水預報精度方面，總體上參數(shù)移植后的精度較移植前較差。

無資料地區(qū)在降雨數(shù)據(jù)移用進行流域水文預報時，除需要考慮地形、距離、下墊面等因素做出相似判斷外，降雨數(shù)據(jù)的同一度分析也至關重要。因此，存在由于忽略降雨中心時空分布而造成的數(shù)據(jù)移用精度不高現(xiàn)象[19-20]。非臺風影響下的降雨資料缺少測站引入集對分析（Set Pair Analysis, SPA）[21-22]，對相鄰站點降雨數(shù)據(jù)進行確定性和不確定性研究分析，從同、異、反3個角度對降雨序列集合進行定量特性分析測站之間的關聯(lián)性，可有效提高移用資料的精度。降雨序列是非平穩(wěn)時間序列，需要尋求一種能從非平穩(wěn)信號中感受到數(shù)據(jù)偏差或噪聲的數(shù)據(jù)影響、有效提取具有實際物理意義的特征信號的方法[23]，才能準確地進行降雨序列間的集對分析。研究表明，補充集合經驗模態(tài)（Complementary Ensemble Empirical Mode Decomposition, CEEMD）[24-26]相比于經驗模態(tài)分解（Empirical Mode Decomposition, EMD）及集合平均經驗模態(tài)分解（Ensemble Empirical Mode Decomposition, EEMD）不僅解決了模態(tài)混疊問題，而且對原始信號進行了精確重構，分解效果更好。利用CEEMD將非臺風雨降雨數(shù)據(jù)分解為多個本征模態(tài)函數(shù)（Intrinsic Mode Function,IMF）分量，分別建立集對分析，通過綜合聯(lián)系度判別測站之間的相似性，選擇最佳降雨數(shù)據(jù)移用，提高資料移用的準確性，可為無資料地區(qū)高精度的洪水預報提供參考依據(jù)。

本文針對廣西駝英水庫洪水預報過程中降雨數(shù)據(jù)短缺問題，提出一種基于線性協(xié)整的臺風雨降雨插補和集對分析的非臺風雨降雨資料移用方法，以期為洪水預報提供較準確可靠的降雨數(shù)據(jù)。

1 研究方法

1.1 線性協(xié)整理論

協(xié)整分析研究兩個及以上非平穩(wěn)時間序列間的均衡關系，主要步驟包括平穩(wěn)性檢驗、協(xié)整關系的估計與檢驗及誤差模型的建立，具體步驟見文獻[27]。

（1）平穩(wěn)性檢驗方法。對時間序列進行平穩(wěn)性檢驗，即單位根（Augmented Dickey-Fuller, ADF）檢驗。ADF檢驗中最小二乘法回歸模型有以下3種：

式中：Δyt為 變量y的一階差分；t為時間項；α為 截距；β為回歸系數(shù)；ρ為系數(shù)；i為時間滯后項；k為滯后階數(shù)（AIC和SC準則確定）；Δyt-i為 滯后i時刻的變量y的一階差分；δi為 時間項相關系數(shù)；εt為白噪聲。

在零單位根檢驗時，假設有單位根，若ADF的顯著性檢驗統(tǒng)計量小于相應的置信度（10%，5%，1%），則說明該時間序列為平穩(wěn)序列，否則對降雨序列先進行一階差分后再次進行檢驗，直至降雨序列差分為平穩(wěn)序列，并判斷兩個降雨序列的單整階數(shù)。若兩個時間序列具有相同的單整階數(shù)，則進行以下的協(xié)整關系估計與檢驗。

（2）協(xié)整關系的估計與檢驗。通過最小二乘法對2個降雨序列進行協(xié)整回歸，并對回歸方程的殘差是否具有單位根進行檢驗。若2個時間序列是協(xié)整的，則殘差將是平穩(wěn)的。

（3）誤差修正模型建立。若降雨序列間存在協(xié)整關系，則可根據(jù)兩變量的一階自回歸滯后模型建立誤差修正模型（Error Correction Model, ECM）：

式中：Δyt、Δxt分 別為y、x變量在t時刻的一階差分；β0、 β1為 回歸系數(shù)；et-1為兩序列協(xié)整關系的前一時刻計算殘差；γ為修正系數(shù)，即因變量拉回長期均衡狀態(tài)的速度。

1.2 集對分析（SPA）理論

集合對H(A，B)是指由兩個相互關聯(lián)的集合A和B組成的一對集合。這兩個集合的性質包括同一度、差異度和對立度[28]，其表達式如下：

式（5）經簡化后可得：

式中：μ為集合對的綜合聯(lián)系度；S為同一的個數(shù)；F為差異的個數(shù)；P為對立的個數(shù)；N為特征個數(shù)總數(shù)；a為同一度；b為差異度；c為對立度；m為差異不確定性系數(shù)；n為對立系數(shù)；在計算中一般取m=0.5，n=-1。

根據(jù)綜合聯(lián)系度判斷目標站點的場次降雨序列與相鄰站點降雨序列的相似性，進行綜合聯(lián)系度高的降雨數(shù)據(jù)的移用，并非距離近的站點間的降雨數(shù)據(jù)移用，這在很大程度上取決于降雨的時空分布。

1.3 補充集合經驗模態(tài)（CEEMD）分解方法

CEEMD算法通過對原始信號加入白噪聲，有效消除重構信號中的殘留噪聲，具體步驟[29-30]如下：

（1）在原始信號中分別加入n種固定比例的高斯白噪聲，對加噪的信號進行分解，得到n個IMF，將得到的各IMF進行加總平均得IMF1：

式中：IMF1為 原信號的第一階分量；N為分解次數(shù)； ωi(t) 為隨機高斯白噪聲；E1表示EMD分解結果中第一個模態(tài)函數(shù)分量。

從原始信號中減去IMF1得到一階殘差r1(t)：

（2）對殘差重構待分解信號重新分解，即可得到原始信號的第二模態(tài)函數(shù)IMF2：

式中：r1(t)+εE1(ωi(t)) 為重構分解信號。

對k=2, 3，···，k，分別計算第k階殘差：rk=rk-1-IMFk， 根據(jù)式（9）得出rk(t)+εEk(ωi(t))的第1個分量，通過對IMF平均得：

（3）反復迭代，直到殘差的極值個數(shù)不超過2：

式中：R為 最終殘差值；K為最終IMF分量數(shù)。

1.4 SCS-CN洪水預報方法

SCS-CN模型是基于水量平衡方程和兩個基本假設建立的[31]。由于其模型結構簡單，參數(shù)較少，被廣泛應用于無資料中小河流的洪水預報中[32-33]，詳細的洪水預報模型原理和過程見參考文獻[4]，此處不再贅述。

本文采用絕對誤差、相對誤差、相關系數(shù)和納什系數(shù)4個指標評定模型結果。相關系數(shù)用來反映模擬值和實測值之間的線性相關程度，納什系數(shù)可以描述洪水預報過程與實測過程之間的擬合程度。

2 實例應用

馱英水庫位于廣西寧明縣那堪鄉(xiāng)垌中村上游約6 km的珠江流域西江水系明江支流公安河上游河段（圖1），流域經緯度范圍為21°35′45′′N～21°55′5.9′′N，107°17′7.1′′E～107°34′3.7′′E，流域面積606 km2，屬南亞熱帶季風氣候[34]。流域內83.87%面積的土壤為壤土，飽和導水率在1.8～18 mm/h范圍內。夏末秋季，受南海及太平洋臺風影響，工程區(qū)多出現(xiàn)臺風雨，且臺風雨強度大、雨量集中。

圖1 馱英水庫上游流域示意Fig. 1 Schematic diagram of the upper reaches of Tuoying Reservoir

由于特殊的地理位置和氣候特征，馱英水庫流域降雨分為臺風雨和非臺風雨2種類型。馱英水庫未建前曾設有馱英水文站，水文站目前已停測，存有1987—1999年共計13 年的水文氣象資料。馱英水庫上游流域已建九特、那馱、板固、叫弄、枯強、潭昔、小平共7個雨量站，各站實測降雨是水庫來水預報的重要數(shù)據(jù)來源，現(xiàn)僅有2012—2018年實測資料。數(shù)據(jù)的不完整性給洪水預報模型的建立和率定帶來了一定挑戰(zhàn)。由于臺風雨與非臺風雨的時空異質性差異較大，采用線性協(xié)整理論與SPA法分別對臺風雨和非臺風雨兩種情況下的降雨進行插補移用計算?；诰€性協(xié)整理論的臺風雨選取板固和叫弄站2013年、2014年和2016年的3場降雨資料計算；非臺風雨選取那馱、九特、板固站2013年、2015年、2017年和2018年的4場降雨序列進行SPA分析降雨數(shù)據(jù)移用計算。

2.1 基于線性協(xié)整關系的臺風降雨計算

由于降雨序列為非平穩(wěn)序列，臺風影響下一定范圍內雨量站的原始降雨序列間存在“偽回歸”的線性關系。因此，需進行協(xié)整分析。首先對降雨序列進行ADF檢驗，以2013年某場降雨為例，表1為臺風雨過程下的流域上游板固站與叫弄站降雨序列的ADF檢驗結果。從表1可以看出，板固站與叫弄站的AIC及SC值均在無截距無趨勢模型達到最小，且兩站的降雨經一階差分后均無單位根，這說明經過一階差分后的降雨序列成為平穩(wěn)序列，且都為一階單整序列。

表1 板固站及叫弄站降雨序列ADF檢驗結果Tab. 1 ADF test results of rainfall series at Bangu Station and Jiaonong Station

用普通最小二乘法（OLS法）估計上述雨量站兩降雨序列的長期均衡方程：

對回歸方程的殘差進行單位根檢驗，檢驗結果見表2。殘差平穩(wěn)性檢驗按照第1種情形（有截距有趨勢）接受原假設、第2種情形（有截距無趨勢）接受原假設、第3種情形（無截距無趨勢）拒絕原假設的流程說明殘差是平穩(wěn)性序列，從而進一步證明了在該場臺風雨影響下板固站與叫弄站的降雨序列滿足長期線性協(xié)整關系。

表2 殘差序列ADF檢驗結果Tab. 2 ADF test results of residual sequence

基于二者的線性協(xié)整關系建立ECM誤差修正模型，其中，ect-1為兩個降雨序列協(xié)整關系的前一時刻計算殘差，并利用OLS法進行誤差修正模型的回歸分析：

根據(jù)板固站的降雨序列，利用ECM模型計算叫弄站的降雨序列，計算結果見圖2。場次雨量峰值的絕對誤差為0.7 mm，計算結果和實測降雨的最大絕對誤差也僅為4.0 mm，且降雨峰值出現(xiàn)時間相同；計算結果與實測數(shù)據(jù)的納什系數(shù)為0.901，相關系數(shù)為0.979，主雨峰相對誤差1.5%，場次降雨總量相對誤差0.3%，誤差較小滿足要求，這表明該方法可以用于降雨序列間的插補計算。

圖2 叫弄站協(xié)整計算結果與實測降雨對比Fig. 2 Comparison of co-integration calculation results and measured rainfall at Jiaonong Station

利用協(xié)整理論對2014和2016年的兩場降雨計算分析，模擬結果見圖3。可以看出，總體計算結果較好。2014年場次降雨協(xié)整計算結果和實測數(shù)據(jù)相比，最大絕對誤差9.4 mm，納什系數(shù)0.892，相關系數(shù)0.955，主雨峰相對誤差6.25%，場次降雨總量相對誤差8.8%；2016年某場降雨與真實數(shù)據(jù)對比，最大絕對誤差6.9 mm，納什系數(shù)0.886，相關系數(shù)0.958，主雨峰相對誤差14.69%（實測降雨47.0 mm，模擬降雨40.1 mm），場次降雨總量相對誤差8.1%；協(xié)整計算雨量與實測數(shù)據(jù)擬合效果較好。由此可見，雨量信息缺失值可采用測站間協(xié)整理論插值得出，提高降雨序列完整性。

圖3 叫弄站協(xié)整計算結果與實測降雨對比Fig. 3 Comparison of co-integration calculation results and measured rainfall at Jiaonong Station

2.2 基于SPA的非臺風降雨移用

那馱、九特、板固3個雨量站間彼此相鄰，2018年某一場降雨實測過程如圖4所示?？梢钥闯瞿邱W站的降水峰值處于板固站和九特站之間；該場降雨中，板固站降水峰值時刻更接近那馱站雨量峰值時刻，板固站和那馱站相關系數(shù)為0.907，九特和那馱的相關系數(shù)為0.725，板固站降雨序列和那馱站的相關性更好。從降水時空分布初步考慮，板固站的降水數(shù)據(jù)移用到那馱站更為合理。

圖4 各雨量站場次降雨數(shù)據(jù)對比Fig. 4 Comparison of rainfall data of each rainfall station

選取非臺風雨影響下的那馱、九特、板固站記錄的場次降雨序列分別進行CEEMD分解，結果見圖5。CEEMD將降雨原序列分解為4個IMF分量，每個IMF的頻率、振幅、周期都不同；從IMF1到IMF4頻率和振幅變小，周期變長；IMF1和IMF2的振幅和頻率較大，周期較短，IMF3和IMF4的振幅和頻率較小，周期較長。

圖5 CEEMD分解結果Fig. 5 CEEMD decomposition results

每個IMF分量的集對分析結果見表3。從表3可以看出，本場降雨過程中，那馱-九特的同一性較高，而差異性和對立性低于那馱-板固站，具體表現(xiàn)在那馱站-板固站的平均綜合聯(lián)系度0.893，那馱-九特的平均綜合聯(lián)系度0.732，在數(shù)據(jù)移植時更適合將板固站的資料移用到那馱站。盡管九特站與那馱站的距離更近，但場次降雨相似性在很大程度上取決于降雨的時空分布。

表3 IMF分量集對分析結果Tab. 3 Analysis results of IMF component set pairs

假設將九特站與板固站記錄的該場降雨過程分別移植至那馱站，并參與面雨量計算進行洪水預報。建立的SCS-CN模型經率定后納什系數(shù)為0.87，滿足預報精度要求，可用于洪水模擬研究。移用數(shù)據(jù)的預報結果與真實那馱站降雨資料預報結果見圖6。結果表明：利用九特站移用的降雨數(shù)據(jù)的水文預測值與真實數(shù)據(jù)預測值之間的納什系數(shù)為0.847，自板固站移用數(shù)據(jù)的模擬結果與真實預報值間的納什系數(shù)為0.905。因此，板固站的移用數(shù)據(jù)模擬效果好于九特站。在洪峰預測精度上，實測洪峰為563 m3/s，采用板固站移用情況下預測洪峰為591 m3/s，相對誤差為4.97%；采用九特站移用情況下洪峰為503 m3/s，相對誤差為10.70%。這說明自板固站移用降雨數(shù)據(jù)更能保證洪峰流量的準確性，可避免因低估洪峰而發(fā)生重大生命財產損失。

圖6 洪水預報模擬結果對比Fig. 6 Comparison of flood prediction simulation results

那馱、板固、九特降雨資料中，存在多場降雨那馱-板固的綜合聯(lián)系度好于那馱-九特的現(xiàn)象。另外，分析了2013年、2015年和2017年3場降雨測站間集對關系。結果表明：2013年那馱-板固的綜合聯(lián)系度為0.863，略好于那馱-九特的綜合聯(lián)系度（0.836）；2015年那馱-板固的綜合聯(lián)系度（0.866）明顯高于那馱-九特的（0.741）；2017年那馱-板固的綜合聯(lián)系度（0.816）略高于那馱-九特的（0.786）。將板固站的降雨資料移用到那馱站作為洪水預報輸入，模擬結果都好于九特站移用結果。因此，非臺風雨影響下的目標站點降雨資料移用時，提前對相鄰站點同時段現(xiàn)有序列分別進行集對分析，根據(jù)綜合聯(lián)系度指標衡量移用對象在此場降雨中與目標站降雨序列的相似性，可提高移用數(shù)據(jù)的可靠性和合理性。

3 結 語

針對駝英水庫洪水預報過程中降雨數(shù)據(jù)短缺現(xiàn)象，提出了基于線性協(xié)整和集對分析分別對臺風雨和非臺風雨降水序列插補移用，并用此方法對馱英水庫進行了洪水預報研究。主要結論如下：

（1）基于線性協(xié)整關系的臺風雨計算，測站之間滿足線性協(xié)整關系，協(xié)整理論模擬值與實測降雨的納什系數(shù)均在0.85以上，相關系數(shù)均大于0.90，絕對誤差、主雨峰相對誤差和場次降雨總量誤差均較小，這表明該方法可以插補降雨序列，提高降雨序列的完整性。

（2）非臺風降雨移用采用SPA，各個雨量站經CEEMD分解，將IMF的分量采用集對分析，根據(jù)綜合聯(lián)系度計算結果移用最佳的降雨資料，利用移用降雨數(shù)據(jù)進行洪水預報，其納什系數(shù)達到0.85，洪峰預測值相對誤差較小。由此可見，SPA方法為無資料地區(qū)降雨數(shù)據(jù)移用提供了一種新思路。