南水北調(diào)中線工程冬季輸水冰情風險研究

練繼建，楊德明，趙 新

(1.天津大學水利工程仿真與安全國家重點實驗室，天津 300072； 2.天津大學建筑工程學院，天津 300350)

南水北調(diào)中線工程受冬季寒冷天氣的影響，渠道經(jīng)常會發(fā)生不同程度的冰情[1-4]，從而改變渠道的水力條件，造成水工建筑物等設(shè)施損壞，不僅嚴重威脅工程冬季輸水的運行安全，還會造成巨大的經(jīng)濟損失。冰情分為多個階段，其中初期是冰情演變的基礎(chǔ)，尤為重要。本文研究的冰情風險指渠道剛開始的產(chǎn)冰階段(即冰期初期階段，不是冰塞冰壩等冰害形成階段)的冰情風險，該階段由于受外界氣象條件的影響，具有較強的隨機性和不確定性，為渠道冰情的風險評估帶來了極大挑戰(zhàn)?！度珖鴼庀蟀l(fā)展“十四五”規(guī)劃》中提出“建立健全以氣象災(zāi)害預警信息為先導的全社會快速響應(yīng)和高效聯(lián)動機制”的方針，因此，通過氣象數(shù)據(jù)對輸水渠道進行冰情風險研究具有重要意義。

對冰情的演變過程已有大量研究。例如：Shen[5]提出了雙層冰輸運數(shù)學模型對表面冰和浮冰進行數(shù)值模擬；Wang等[6]基于二次成核和絮凝等多個過程，開發(fā)了水的過冷過程演變模型和水內(nèi)冰演變模型；宗全利等[7]根據(jù)熱平衡理論，綜合考慮水力、熱力、氣候等條件，推導得到渠道不凍長度的計算公式；唐偉[8]通過地下水抽灌的地熱能利用方式，結(jié)合熱平衡理論對輸水渠道無冰蓋輸水的可行性進行了分析；穆祥鵬等[1]根據(jù)傳熱學理論，對鋪設(shè)保溫蓋板的輸水渠道水體進行了數(shù)學模擬。由上述研究可以看出，熱平衡過程直接影響著冰的形成和演變。

對輸水工程冬季輸水潛在風險的研究主要在風險因素識別[9-12]和風險評價方法[11-16]兩個方面開展。姜蓓蕾等[17-18]利用層次分析法構(gòu)建評價指標體系，并采用模糊綜合評價法對南水北調(diào)東線工程和中線工程某段進行了風險評價。胡丹等[19]構(gòu)建了直覺模糊集理論評價模型，采用TOPSIS法對輸水渠道進行綜合評價并計算了相應(yīng)的風險等級。劉勇等[20]針對長距離輸水工程的特點，考慮風險的預測、可控和轉(zhuǎn)移等，構(gòu)建了多維功效函數(shù)對工程單元進行風險識別?，F(xiàn)有的輸水風險評價關(guān)注點主要集中在輸水建筑物上，缺乏輸水渠道冰情過程的風險評價，也忽視了采用概率理論對風險因素進行宏觀分析。概率分析可以較為直接地描述風險因素之間的內(nèi)在聯(lián)系，已得到廣泛應(yīng)用[21-24]。劉強等[25]針對寒區(qū)公路的風吹雪災(zāi)害問題，利用概率論思想，構(gòu)建聯(lián)合密度函數(shù)進行脆弱性評估，并在此基礎(chǔ)上建立風險評估體系進行風險分析。吳柳萍等[26]基于信息擴散理論，建立了福建省森林病蟲害風險概率模型，對森林病蟲害發(fā)生的風險范圍進行了評估。王磊之等[24]構(gòu)建了暴雨-潮位聯(lián)合分布模型分析深圳市風險概率，從而定量揭示了極端降水與潮位遭遇風險的內(nèi)在規(guī)律。

本文以渠道熱平衡理論為基礎(chǔ)，對評價指標進行概率分析，從多角度對冬季輸水冰情進行風險評估，以期為輸水工程安全運行和防災(zāi)減災(zāi)工作提供參考。

1 研究方法與研究數(shù)據(jù)

1.1 研究方法

針對輸水工程冬季輸水冰情風險，以熱平衡理論為出發(fā)點，計算出整個冰情的形成與演變過程中各部分熱交換量，推導出目標單元渠段的整體總熱量變化計算公式，并從中找出宏觀風險性評價指標；基于冰情風險評價指標與風險發(fā)生結(jié)果相對應(yīng)的原則，通過總熱量變化率(每米渠道單位水面面積的總熱量變化)和風險評價指標的聯(lián)合概率函數(shù)求出冰情的風險等級；根據(jù)風險等級對南水北調(diào)中線工程(河北段)進行風險評估，分析該渠段的熱量變化規(guī)律和風險等級分布。

1.2 研究數(shù)據(jù)

以南水北調(diào)中線工程河北段為研究區(qū)域，將其分為安陽節(jié)制閘—銘河節(jié)制閘、銘河節(jié)制閘—午河節(jié)制閘、午河節(jié)制閘—沙河節(jié)制閘和沙河節(jié)制閘—渠末4段，工程位置如圖1所示。氣象數(shù)據(jù)采用 2000—2018年邯鄲、保定、邢臺、石家莊4個城市的氣溫、風速和太陽總輻射數(shù)據(jù)，來源于國家氣象科學數(shù)據(jù)中心和國家青藏高原科學數(shù)據(jù)中心。地圖矢量數(shù)據(jù)和柵格數(shù)據(jù)來源于地理空間數(shù)據(jù)云和國家地理信息系統(tǒng)網(wǎng)。

圖1 南水北調(diào)中線工程河北段工程位置示意圖

2 熱平衡理論計算

對于沒有封凍的輸水渠道，不考慮水體與渠壁和渠床的交換以及外界特殊天氣變化(降雨、下雪等)，只考慮水面的熱交換和水流的動力增熱，其中水面的熱交換包括短波輻射、凈長波輻射、蒸發(fā)散熱和對流散熱。根據(jù)熱平衡理論，單元渠段的總熱量變化可表示為

W總=W水+W熱交換+W動力=W水+

W短波-W長波-W蒸發(fā)-W對流+W動力

(1)

式中：W總為水體總的熱量變化，kJ；W水為水體因水溫而產(chǎn)生的自身熱量，kJ；W熱交換為水面的總熱交換，kJ；W動力為水流的動力增熱，kJ；W短波、W長波、W蒸發(fā)、W對流分別為水體的短波輻射熱量、凈長波輻射熱量、蒸發(fā)散熱和對流散熱，kJ。

當輸水渠道在冬季運行時，假設(shè)整個渠道恒定運行且初始沿程水溫近似為0℃，即認為W水=0 kJ，同時可以忽視水力因素(較小的流量、水位變化等)對水體熱量變化的影響，則式(1)可表達為

W總=W短波-W長波-W蒸發(fā)-W對流+W動力=

DLσ總=DL(σ短波-σ長波-σ蒸發(fā)-σ對流+σ動力)

(2)

其中σ短波=(1-φ)(1-0.65C2)Rs

σ長波=0.97α[εa(273+Ta)4-(273+Tw)4]

σ蒸發(fā)=f(vw)(es-ea)

σ對流=0.47f(vw)(Tw-Ta)

此外，一個長距離輸水渠道一般由多個單元輸水渠段組成，則整個輸水渠道的總熱交換量可由各單元輸水渠段的總熱交換量求和得出。當各單元輸水渠段總熱量變化為正值時，則整個單元渠道基本不會產(chǎn)生水內(nèi)冰和流冰，從而保證整個輸水渠道安全運行。

3 風險評估

3.1 評價指標選取

輸水渠道出現(xiàn)冰情是多方面原因造成的，既與渠道自身的工程特性和水流的運行狀態(tài)等因素密切相關(guān)，同時也受到外界天氣環(huán)境的影響。而當渠道在某一恒定流工況無冰運行時，其水力條件基本沒有發(fā)生改變，而外界環(huán)境變化是造成渠道出現(xiàn)大量水內(nèi)冰和流冰而產(chǎn)生冰情的主要原因。結(jié)合式(2)及實際工程運行資料，可知氣溫Ta、風速vw、太陽總輻射Rs3個氣象參數(shù)是影響熱量變化的關(guān)鍵評價指標，也是渠道冰情風險主要致災(zāi)因子。

本文采用單因素局部分析法[27-28]，依據(jù)式(3)分析判斷氣溫、風速、太陽總輻射3個參數(shù)在熱量變化計算公式中的敏感性，結(jié)果如表1所示。

表1 3個參數(shù)的敏感性分析結(jié)果

(3)

式中：Sr為參數(shù)的敏感性指標；P1、P2分別為參數(shù)變化前后的取值；σ1、σ2分別為采用P1、P2計算得到的熱量變化率。其中Sr絕對值越小，敏感程度越小，代表著該參數(shù)對計算結(jié)果的影響程度越小。

從表1可以看出，3個參數(shù)的敏感性按照從大到小順序排列為氣溫、風速、太陽總輻射，由此得出太陽總輻射對渠道熱量變化影響最小。同時還發(fā)現(xiàn)太陽總輻射的Sr不會隨太陽總輻射值變化而改變，這是因為太陽總輻射與熱量變化呈線性關(guān)系；而氣溫的Sr變化較小，風速的Sr會隨著風速的增大而減小，說明風速變化越大，對熱量變化的影響越小。

3.2 概率風險分析

通過概率分析初步計算致災(zāi)因子同時發(fā)生的事件概率，再利用式(2)求出總熱量變化率從而進一步判斷渠道冰情風險等級?？紤]到氣象數(shù)據(jù)在時間和空間上的有效性和連續(xù)性，選取2000—2018年每年的11月到次年2月的日數(shù)據(jù)作為總體樣本，然后根據(jù)樣本數(shù)據(jù)繪制出風險評價指標的經(jīng)驗分布曲線和相應(yīng)的概率密度直方圖，通過二者相互匹配擬合出風險評價指標的概率密度函數(shù)。南水北調(diào)中線工程河北段主要途經(jīng)邯鄲、邢臺、石家莊、保定4個城市，分別對其氣溫、風速和太陽總輻射數(shù)據(jù)的概率密度曲線進行擬合。以邢臺為例，其氣溫、風速和太陽總輻射經(jīng)驗分布曲線和概率密度直方圖如圖2和圖3所示，利用Spass數(shù)學分析軟件對經(jīng)驗分布曲線和概率密度直方圖進行擬合，方差檢驗結(jié)果表明，氣溫和風速的概率密度曲線呈一階指數(shù)形式，因此采用Gauss函數(shù)作為擬合曲線的表達形式，4個城市的擬合表達式及其擬合度如表2所示。

圖2 邢臺市氣溫、風速和太陽總輻射經(jīng)驗分布曲線

圖3 邢臺市氣溫、風速和太陽總輻射概率密度直方圖

由于季節(jié)影響，12月太陽總輻射值小于11月和1月的輻射值，通過預處理和方差分析發(fā)現(xiàn)其概率密度曲線呈多峰形態(tài)，由多個一階指數(shù)組成。以邢臺為例，其概率密度曲線由兩個Gauss一階指數(shù)函數(shù)疊加組成，如式(4)所示；其余3個城市疊加函數(shù)過多，不在此贅述。

表2 4個城市氣溫和風速的概率密度函數(shù)及其擬合度

F(Rs)=F(Rs)1+F(Rs)2=

(4)

式中：F(Rs)為太陽總輻射的總概率密度函數(shù)；F(Rs)1、F(Rs)2為太陽總輻射的子概率密度函數(shù)。

在求得氣溫、風速、太陽總輻射3個評價指標的概率密度函數(shù)后，基于指標之間的相關(guān)性，利用其復相關(guān)系數(shù)構(gòu)造聯(lián)合概率密度函數(shù)[29-30]，并利用積事件的發(fā)生概率來代表3個評價指標在不同強度下同時發(fā)生的可能性：

(5)

式中：F(Ta,vw,Rs)為氣溫、風速、太陽總輻射3個評價指標不同強度下同時發(fā)生的聯(lián)合概率密度函數(shù)；r為風速、太陽總輻射對氣溫的復相關(guān)系數(shù)。通過數(shù)據(jù)的線性回歸計算，得出邯鄲、邢臺、石家莊、保定4個城市風速、太陽總輻射對氣溫的復相關(guān)系數(shù)分別為0.388、0.429、0.475和0.243，同時相對應(yīng)的ANOVA方差分析回歸系數(shù)分別為0、0、0和0.028，均小于顯著性水平0.05，表明分析的線性模型是成立的。

對于渠道冰情風險，總熱量變化率與冰情的風險等級呈負相關(guān)關(guān)系，一般情況下總熱量變化率越小，代表著其冰情的風險等級越高，風險性越大。圖4為4個城市氣溫、風速、太陽總輻射3個評價指標共9 638個工況組合計算得到的總熱量變化率概率曲線。依據(jù)GB/T 27921—2011《風險管理 風險評估技術(shù)》，按照等級劃分標準將總熱量變化率劃分為5個等級：>7 750 kJ/m2為極低風險，>2 000～7 750 kJ/m2為低風險，>-1 000～2 000 kJ/m2為較高風險，>-10 500～-1 000 kJ/m2為高風險 ，≤-10 500 kJ/m2為極高風險，與冰情風險等級相對應(yīng)。

圖4 輸水渠道總熱量變化率概率曲線

為了更好地將氣象事件發(fā)生的可能性與風險性相關(guān)聯(lián)，根據(jù)9 638個工況組合繪制出不同渠段的總熱量變化率-概率散點圖，并對散點進行了曲線擬合，如圖5所示。通過對該擬合曲線進行積分，可以求出不同渠段不同風險等級發(fā)生的概率，如表3所示。從表3可以看出風險等級為極高和高的發(fā)生概率很小，以風險等級為低和較高的發(fā)生概率為主。

圖5 4個渠段總熱量變化率-概率散點圖及其擬合曲線

表3 4個渠段不同冰情風險等級發(fā)生概率

考慮氣溫、風速和太陽總輻射3個評價指標發(fā)生概率與冰情風險的內(nèi)在聯(lián)系，將3個評價指標也按照工程標準進行風險等級劃分(表4)，并分析4個渠段處在極高和高風險冰情情況下，氣溫、風速和太陽總輻射3個評價指標所屬的風險等級，結(jié)果如表5所示。從表5可以看出當渠道處在極高和高風險冰情時，氣溫和太陽總輻射都處在高風險以上，而風速則處在較高風險以上。結(jié)合表1的敏感性分析結(jié)果，可能是因為風速變化越大，對熱量變化的影響越小，而氣溫于對于渠道冰情演變一直具有較大影響。

表4 3個評價指標的風險等級劃分

表5 4個渠段在極高和高風險冰情下3個評價

4 實例分析

以2021年1月6—8日氣象條件為基礎(chǔ)，對南水北調(diào)中線工程河北段的冬季輸水進行初步的風險分析。

4.1 整體風險分析

表6為4個渠段熱量變化率及冰情風險等級計算結(jié)果。從表6可以看出，3天中4個渠段的總熱量變化率都是負值，說明渠段水體一直在散失熱量，有發(fā)生冰情的潛在風險。沙河節(jié)制閘—渠末段的熱量變化率最大，遠大于其他3個渠段，說明該渠段的冰情風險最大。在整個熱交換過程中，動力增熱和凈長波輻射占比很小，幾乎可以忽略不計，而短波輻射、蒸發(fā)散熱和對流散熱明顯占據(jù)主導地位，是影響熱量變化的主要因素。

根據(jù)實際監(jiān)測結(jié)果，1月6日前渠道無流冰現(xiàn)象，自7日開始渠道出現(xiàn)表面流冰層，產(chǎn)生冰情，這可能是由于渠道受“拉尼娜”寒流影響發(fā)生強降溫，渠道水體熱量散失嚴重所致。從表6可以看出，6、7日4個渠段的冰情風險等級基本都為高風險，而到了8日隨著天氣回暖，冰情風險等級也隨之變小。但是，6日和8日沙河節(jié)制閘—渠末段的總熱量變化率明顯高于其他渠段，以致風險等級為極高和高，其風險性高于其他渠段，因此應(yīng)加強對該段的實時監(jiān)測，采取相關(guān)冰情防護措施。

表6 4個渠段熱量變化率和風險分析結(jié)果 單位：kJ/m2

根據(jù)式(5)計算得出這3天的氣象條件聯(lián)合發(fā)生概率很小，這既體現(xiàn)了極端天氣發(fā)生的可能性較小，又間接驗證了表3中極高和高風險發(fā)生概率所占比例最小這一結(jié)論。結(jié)合表5，可得出各評價指標風險等級與冰情風險等級的對應(yīng)關(guān)系如表7所示。從表7可以看出7日各渠段氣溫和風速的風險等級基本屬于較高風險以上，而太陽總輻射都屬于極高風險。在太陽總輻射風險等級不變的情況下，雖然風速風險等級波動較大，但是各渠段的冰情風險等級都屬于高風險，這也證明了氣溫影響冰情風險程度遠超其他兩個評價指標，與表1和表7結(jié)果一致。

表7 評價指標和冰情風險等級的對應(yīng)關(guān)系

由此可見，冬季平均氣溫越低、風速越大、太陽總輻射越小，則渠道水體散失熱量也就越多，而冰情的風險等級也就越來越高，對渠道運行造成的危害也就越大。此外，通過對渠道風險分析可知，即使氣溫、風速和太陽總輻射3個評價指標中某個指標處于低風險等級，但渠道仍有發(fā)生高風險冰情的可能性，因此在實際管理中應(yīng)時刻密切關(guān)注氣象條件的動態(tài)變化。

4.2 熱交換的時程變化分析

圖6為4個渠段熱交換量的時程變化曲線。從圖6(a)可以看出冬季輻射熱交換量基本在白天，光照時間在9～10 h之間，在1:00左右達到熱交換量最大值。由于這3天的氣溫和輻射變化幅度不大，因此短波輻射和凈長波輻射的熱交換量近似周期性變化，如圖6(a)(b)所示。從圖6(c)可以看出蒸發(fā)散熱變化趨勢與風速日過程變化曲線相似，這是因為蒸發(fā)熱交換量為風速函數(shù)。同時，從圖6(c)(d)還可以發(fā)現(xiàn)夜間蒸發(fā)和對流熱交換量較大，且沒有輻射熱交換，因此夜間渠道熱量散失嚴重，極易出現(xiàn)流冰。為了進一步了解熱交換量過程變化，繪制出各渠道總熱交換量時程變化曲線和累積總熱交換量時程變化曲線，如圖7所示。

圖6 4個渠段熱交換量的時程變化曲線

圖7 4個渠段總熱交換量和累積總熱交換量時程變化曲線

從圖7(a)可以看出，由于這3天天氣條件逐漸變好，因此各渠段總熱交換量時程變化趨勢是波動向上的。同時可以看出中午和夜間是曲線的峰值，熱交換量達到極值。由圖7(b)可見，各渠段的累積總熱交換量是逐漸變大的，這是因為在這3天總熱交換量基本是負值，安陽節(jié)制閘—銘河節(jié)制閘、銘河節(jié)制閘—午河節(jié)制閘和午河節(jié)制閘—沙河節(jié)制閘這3段的累積總熱交換量差別不大。對比圖7(a)(b)，隨著8日天氣條件相對變好，熱交換量由散失變?yōu)槲?，因此圖7曲線的末端出現(xiàn)短暫的上升趨勢?？梢?，在這3天中各渠段的冰情風險持續(xù)增大，必要時需采取相應(yīng)的措施進行防護。

綜上所述，氣象條件是造成輸水渠道冰情風險的主要因素，直接影響冰情風險等級。當冷空氣來臨時，各渠段冰情風險的可能性比平時要高；而沙河節(jié)制閘—渠末段發(fā)生冰情的敏感性和風險性比其他渠段要高，因此在正常防護過程中，應(yīng)加大對沙河節(jié)制閘—渠末段的監(jiān)測，并加強該渠段的防冰措施，以有效降低冰情風險。

5 結(jié) 語

通過對南水北調(diào)中線工程河北段冬季輸水的冰情風險研究，發(fā)現(xiàn)短波輻射、蒸發(fā)散熱和對流散熱是影響熱量變化的主要因素。在對冰情風險進行概率分析和風險評價時，不僅可以判斷當下氣象事件發(fā)生的可能性，而且還可以迅速識別出冰情的風險等級，并分析出氣象事件發(fā)生概率與冰情風險等級的內(nèi)在聯(lián)系，從而為今后工程的冰情防治提供科學的決策依據(jù)。同時，針對實際輸水工程，冰情風險評估是一個相對動態(tài)的過程，會隨著運行工況和氣象情況變化而變化。因此，面對動態(tài)化的風險評估，如何能夠更迅速準確量化和測評冰情風險是今后需要深入研究的課題。