基于Kriging插值的校園內(nèi)澇模擬與模型優(yōu)選

楊亞紅，楊興峰，閆俊江，蔚 陽

1)蘭州理工大學土木工程學院,甘肅蘭州730050；2)清華大學深圳研究生院，廣東深圳518055

隨著全球氣候變化加劇，突發(fā)性強降雨引發(fā)的城市內(nèi)澇問題日趨嚴重，給人們的生命財產(chǎn)和社會經(jīng)濟發(fā)展構(gòu)成巨大威脅[1]．據(jù)統(tǒng)計，2008—2018年期間，中國平均每年因洪澇災害死亡人口約為760人，直接經(jīng)濟損失為2 119.50億元[2]．對暴雨引發(fā)的城市內(nèi)澇進行風險評估，可針對性地進行預防和治理，是降低城市內(nèi)澇災害損失的一條有效途徑．

國內(nèi)外對于洪澇災害風險評估的研究較多．BENITO等[3]綜合古代洪水信息(100～10 000 a)、歷史洪水數(shù)據(jù)(1 000 a)以及水文站數(shù)據(jù)(30～50 a)，提出了基于歷史洪災數(shù)據(jù)的洪水風險評估方法，并將其應(yīng)用于歐洲國家．WASHAKH等[4]綜合分析1990—2018年Arun流域49個冰湖的冰湖潰堤洪水(glacial lake outburst flood, GOLF)的氣溫、降水趨勢與所記錄GLOF事件發(fā)生的相關(guān)性，提出了一種簡單易行的GOLF風險評估模型．20世紀初，隨著地理信息系統(tǒng)(geographical information system, GIS)的不斷發(fā)展，基于情景模擬的綜合性風險評估與分析開始受到普遍關(guān)注．應(yīng)用該方法時，給定災害初始條件，根據(jù)災害演化機理進行模擬推演，得到災害的發(fā)展演化過程，從而估計災害損失．該分析方法能反映出災害的動態(tài)演化，并對以前從未發(fā)生過的非常規(guī)災害進行后果評價，其評估對象涉及多個方面，可預測建筑物、經(jīng)濟與社會等各方面的損失情況．AMEER等[5]利用情景模擬演化的方法對水土評估模型(soil and water assessment tool， SWAT)的模擬精度進行改進．此外，國內(nèi)學者也嘗試開展了一些基于情景模擬的洪澇災害風險評估．尹占娥等[6]基于情景模擬的方法針對上海市靜安區(qū)小尺度暴雨內(nèi)澇進行了風險評估．胡蓓蓓等[7]根據(jù)未來規(guī)劃，對天津市濱海新區(qū)2020年的暴雨內(nèi)澇風險進行了預測評估．

為使評估結(jié)果更接近實際，評估的時空精細化程度能夠達到事半功倍的效果．如何提高評估的時空精細化，成為目前風險內(nèi)澇災害評估的發(fā)展趨勢．目前，該模擬評估方法主要針對沿海城市小尺度流域研究較多，而針對內(nèi)陸地區(qū)研究較少且不夠深入．本研究建立一種基于暴雨洪水管理模型-地理信息系統(tǒng)(storm water management model with geographic information system， SWMM-GIS)的水文學模型的城市內(nèi)澇情景模擬風險評估方法，以中國西北地區(qū)某校園節(jié)點積水數(shù)據(jù)為例進行內(nèi)澇風險評估，模擬該地區(qū)在特定重現(xiàn)期降雨情景下的城市內(nèi)澇分布，并對可能造成的災害損失情況進行分析．旨在篩選出適合于該區(qū)域的內(nèi)澇插值模型方法，為利用有限觀測點更有效地識別出內(nèi)澇區(qū)域分布提供基礎(chǔ)．

1 材料和方法

1.1 基于Kriging法的情景模擬方法

1.1.1 基于Kriging插值的城市內(nèi)澇情景模擬

降雨徑流模型考慮降雨的物理過程，通過水力學公式可計算降雨引發(fā)內(nèi)澇的時空分布，而且模型對城市降雨積水的具體特征做了大幅簡化，在保證一定的計算精度的前提下大幅降低了計算時間．目前，世界范圍內(nèi)有許多可用的城市降雨徑流模型，例如，SWMM和伊利諾城市排水區(qū)域模擬(Illinois urban drainage area simulator, ILLU-DAS)等[8]．這類模型基本可以滿足城市暴雨排水系統(tǒng)模擬的需要．但由于是一維水動力模型，所以不適合做內(nèi)澇評估分析．

Kriging插值法的提出為內(nèi)澇時空分布評估提供了可能．該方法通過對積水點采樣數(shù)據(jù)及地理特征合理分析其可能出現(xiàn)的演化過程，進而選擇合適的擬合模型創(chuàng)建淹沒分析圖．Kriging插值方法不僅能量化已知點之間的空間自相關(guān)性，并可說明采樣數(shù)據(jù)在預測區(qū)域范圍內(nèi)的空間分布情況[9]．

與反距離權(quán)重法類似，Kriging插值法對周圍的測量值進行加權(quán)以得出未測量位置的預測．其計算公式[10]為

(1)

其中，si為第i個積水點位置；Z(si)為第i個位置的測量值；λi為第i個位置處的測量值的未知權(quán)重；s0為預測位置；N為測量值數(shù)．

1.1.2 半變異函數(shù)模型優(yōu)選及方法檢驗

1)半變異函數(shù)及擬合模型：Kriging法插值過程中，有限積水點擴散區(qū)域化模型是隨機且規(guī)律性的．因此，模型建立需要由合適的半變異函數(shù)實現(xiàn)從有限至區(qū)域化的相似擴散演化，使半變異函數(shù)曲線與積水點數(shù)據(jù)達到擬合，最終用于模型預測．

半變異函數(shù)計算公式[11]為

(2)

其中，h為距積水點si的距離；γ(h)為h處的積水深度變異函數(shù)；N(h) 為分割距離為h時樣本點的總個數(shù)；Z(si)、Z(si+h)分別為積水點si和si+h處的積水深度測量值．

2)半變異函數(shù)的模型優(yōu)選：Kriging法中，式(2)只能反映出積水點深度本身的變異特征和空間分布，不能反映出在不同空間滯后距離上的空間變異特征以及積水點的空間異質(zhì)性定量化特征，故需要通過常用半變異函數(shù)的擬合模型對已知數(shù)據(jù)點離散關(guān)系進行擬合[11]．擬合模型由半變異函數(shù)根據(jù)已知測量點的空間自相關(guān)性確定，如圖1[12]．目前常用的半變異模型有穩(wěn)定、球面、指數(shù)和高斯模型等．選擇合適的擬合模型是空間插值的前提條件，直接影響空間插值的精度結(jié)果．

圖1 半變異函數(shù)擬合曲線示例[12]Fig.1 Examples of semi-variation fitting curves[12]

圖1中的C0是塊金值，為半變異函數(shù)曲線在γ(h)軸上的截距，表示積水深度受到隨機因素影響的空間變異；C為偏基臺值，即半變異函數(shù)值達到穩(wěn)定時的增量，用于反映積水深度因空間結(jié)構(gòu)引起的變異；C0+C為基臺值，表示積水深度總的空間變異；C0/(C0+C)為基底效應(yīng)，表示空間的變異特征，該值越大，說明積水深度的空間變異更多是由隨機成分引起的；a為變程，即半變異函數(shù)第1次達到穩(wěn)定時所對應(yīng)的采樣距離，表示積水深度的空間自相關(guān)范圍．

最優(yōu)模型主要通過無偏性和一致性來判斷．前者常采用偏差均值(Kriged reduced mean error, KRME)指標來衡量；后者則采用一致性系數(shù)(Kriged reduced mean square error, KRMSE)指標判定．KRME越接近0，表示無偏估計的效果越佳．KRMSE越接近1，表示模型越好[9,12]．KRME與KRMSE計算公式分別為

(3)

(4)

3)檢驗方法：為保證插值精度，在模型優(yōu)選的基礎(chǔ)上進一步采用交叉檢驗(cross-validation)進行插值方法驗證．原理是在給定的樣本模型中匯總之后留出部分樣本用于建立模型進行預測，并求出這小部分樣本的預測誤差[9]．

1.2 研究區(qū)域概況

研究區(qū)域地處蘭州市靠近黃河南岸的七里河區(qū)，如圖2．該城區(qū)位于市內(nèi)中南部，地理坐標為103°36″E—103°54′54″E，35°50′32″N—36°06′12″N.地勢南高北低，平均海拔2 321 m．南部為石質(zhì)山地，地形特征為山高、谷深、坡陡，巖石裸露．七里河區(qū)境內(nèi)屬大陸性半干旱氣候．主要特征是：四季分明，冬夏長，春秋短．雨熱同季，垂直氣候變異顯著．氣溫、熱量和光照隨海拔由南向北升高；年平均降雨量為327 mm，由南到北逐漸降低．校區(qū)占地87 hm2，在開發(fā)之前多為平頂屋面，其數(shù)量少且占地面積有限；周圍多分布山楊、樺樹群體及少量灌木．這類植物滲透系數(shù)高，對水分蒸散與入滲作用顯著．校園開發(fā)后土地利用類型多樣化，其中，綠地面積占比為42%，不滲透面積占比為57%(包括新開發(fā)地塊4%、建筑物10%、開闊空地14%和道路29%)．

圖2 研究區(qū)域范圍(線內(nèi))Fig.2 The study area (in line)

利用校園雨水管網(wǎng)計算機輔助設(shè)計(computer aided design, CAD)圖將研究范圍劃分為子匯水區(qū)226個．經(jīng)統(tǒng)計，研究區(qū)域概化管道總長為8 804.55 m，匯水區(qū)總面積64.11 hm2，概化結(jié)果如圖3．

圖3 研究區(qū)域概化結(jié)果Fig.3 The regional generalization results of study area

1.3 暴雨強度公式的確定及設(shè)計降雨過程線

因缺乏蘭州市歷年降雨氣象水文資料，所以本研究降雨數(shù)據(jù)由暴雨強度公式生成，該公式由王杰等[13]根據(jù)蘭州市國家基本氣象站1960—2014年的降雨資料，通過皮爾遜-Ⅲ型分布擬合得出．降雨強度公式為

(5)

其中，i為降雨強度(單位：mm/min)；P為降雨重現(xiàn)期(單位：年)；t為降雨歷時(單位：min)． 在該模型中，采用芝加哥降雨模式描述降雨過程．這種降雨模式可以概括大多數(shù)降雨類型，能反映出地表徑流與出流量峰值，因此適用于區(qū)域和雨水系統(tǒng)設(shè)計的預評價[14]．研究中采用雨峰系數(shù)為0.291，其主要目的是分析早期降雨高峰對研究區(qū)域徑流和積水的影響．選取1、2、5和7 a重現(xiàn)期(降雨歷時短，為2 h)和10、20、30和50 a重現(xiàn)期(降雨歷時長，為12 h)繪制芝加哥降雨過程線如圖4．本研究中實測降雨數(shù)據(jù)采用2018年6月、7月及2019年6月、7月共4個月的降雨特征和對應(yīng)的流量數(shù)據(jù)，降雨編號依次為201806、201807、201905和201907(表1)．

圖4 不同重現(xiàn)期下芝加哥降雨過程線Fig.4 The precipitation pattern of Chicago in different recurrence periods

表1 實測降雨數(shù)據(jù)

2 結(jié)果與分析

基于以上數(shù)據(jù)，將降雨數(shù)據(jù)導入率定后的SWMM模型中，模擬得到每個重現(xiàn)期下的雨水井節(jié)點積水數(shù)據(jù)．分別取短歷時和長歷時下重現(xiàn)期最大(P=7 a與P=50 a)時的積水情形進行基于Kriging 插值的內(nèi)澇情景模擬，分析其可能出現(xiàn)的積水分布情況．

2.1 節(jié)點積水影響因子Pearson相關(guān)性分析

雨水管網(wǎng)中，管網(wǎng)節(jié)點(雨水井)的內(nèi)澇積水情況與雨水井的分布、內(nèi)底高程及最大深度和井內(nèi)污染物的空間濃度分布等因素密切相關(guān)．不合理的雨水井及管道的空間布置不當均會損壞管網(wǎng)，造成局部區(qū)域的內(nèi)澇積水與污染[15-16]．選取與影響積水深度相關(guān)的節(jié)點內(nèi)底高程與節(jié)點最大深度等節(jié)點結(jié)構(gòu)屬性變量，利用SPSS 19.0軟件進行Pearson相關(guān)性分析，得到研究區(qū)域P=7 a與P=50 a節(jié)點積水的主要影響因子相關(guān)性大小(表2)．其中，E為節(jié)點內(nèi)底標高；H為節(jié)點最大深度；Q1為P=7 a時節(jié)點積水深度；Q2為P=50 a時節(jié)點積水深度．結(jié)果表明，與內(nèi)底標高對積水深度的影響程度相比，節(jié)點最大深度對積水深度影響更大．重現(xiàn)期為7 a和50 a的積水深度標準差分別為1.034和1.043，它們分別與節(jié)點最大深度的標準差均非常接近；Pearson相關(guān)性系數(shù)分別為0.605與0.766，均大于0.5．分析其原因：重現(xiàn)期為7 a(短歷時)和50 a(長歷時)下降雨量分別達到28.33 mm和61.03 mm，按照氣象學上對降雨的劃分屬于大到暴雨[17-18]．當降雨強度最大時，到達地面形成地表徑流可能在校區(qū)南高北低的地勢作用驅(qū)使下直接匯入雨水管網(wǎng)并到達節(jié)點的最大深度，加上后續(xù)徑流的持續(xù)補充，導致節(jié)點可以在短時間內(nèi)保持著與最大深度相差無幾的積水深度，情況嚴重時可能出現(xiàn)溢流現(xiàn)象．

表2 節(jié)點積水深度與各影響因子的Pearson相關(guān)性系數(shù)

2.2 基于Kriging插值的節(jié)點積水深度情景模擬

通過節(jié)點積水影響因子Pearson相關(guān)性分析，選用節(jié)點最大深度作為節(jié)點積水深度的主要影響因子，基于校園數(shù)字高程模型(digital elevation model，DEM)進行Kriging插值情景模擬．

2.2.1 數(shù)據(jù)統(tǒng)計與空間差異分析

通過對原始節(jié)點積水深度數(shù)據(jù)進行統(tǒng)計分析，并繪制出正態(tài)分布圖，見圖5．由圖5可知，P=7 a與P=50 a時的節(jié)點積水深度數(shù)據(jù)均存在明顯的線性關(guān)系且基本符合正態(tài)分布．P=7 a重現(xiàn)期下211個積水節(jié)點深度變化為0～5.40 m，平均深度為1.39 m；P=50 a重現(xiàn)期下220個積水節(jié)點深度變化范圍為0～5.38 m，平均深度為1.50 m．此外，同一降雨情境下易發(fā)生積水的節(jié)點多集中在深度較淺的范圍內(nèi)：P=7 a時，積水深度為0.76～1.52 m內(nèi)有87個積水節(jié)點，占比39.4%；P=50 a積水深度在0.81～1.62 m內(nèi)有127個積水節(jié)點，占比57.7%．結(jié)果表明雨水管網(wǎng)節(jié)點深度設(shè)置不恰當，容易造成局部區(qū)域大面積積水．

圖5 積水點數(shù)據(jù)概括性統(tǒng)計分析Fig.5 Summary statistical analysis of water accumulation point data

2.2.2 半變異函數(shù)模型的對比分析與選擇

通常在Kriging插值模型方法的選擇中，半變異函數(shù)作為統(tǒng)計采樣點間的距離函數(shù)，其步長直接反映插值精度[14,19]．基于以上統(tǒng)計分析，利用ArcGIS10.4中空間分析工具確定點與相鄰元素之間的平均距離，保證參與計算的半方差的的點對數(shù)不少于30．針對P=7 a與P=50 a降雨情境下的節(jié)點積水深度數(shù)據(jù)，分別構(gòu)建穩(wěn)定、球面、指數(shù)與高斯4種半變異函數(shù)模型，并計算各自的特征參數(shù)進行模型優(yōu)選對比．其模型參數(shù)如表3所示．

表3 四種待選模型參數(shù)對比

通過模型擬合對比分析，在兩種降雨情境下得到的積水深度數(shù)據(jù)中，模型曲線與實際值平均化擬合度較好的均為高斯模型．在表3的模型參數(shù)分析中，高斯模型的KRME值較其他模型最接近0，其KRMSE值最接近于1．

綜上所述，4個待選擬合曲線待選模型中，高斯模型最適合P=7 a和P=50 a降雨情境下的節(jié)點積水深度數(shù)據(jù)做插值研究，故優(yōu)選高斯模型．

2.2.3 交叉驗證結(jié)果分析

為檢驗高斯模型用于研究區(qū)域節(jié)點積水深度Kriging插值模擬的預測精度，在ArcGIS10.4的地統(tǒng)計模塊中將研究區(qū)域P=7 a與P=50 a的246個節(jié)點積水深度數(shù)據(jù)隨機構(gòu)建為2個子集要素：65%(N=160)的數(shù)據(jù)作為訓練數(shù)據(jù)集，35%(N=86)的數(shù)據(jù)作為測試數(shù)據(jù)集．以進行節(jié)點的預測值與實測值的交叉驗證分析，得出研究區(qū)域節(jié)點積水深度的空間插值結(jié)果的預測值與實測值散點分布(圖6)．據(jù)圖6，兩個重現(xiàn)期降雨下的節(jié)點積水深度測量值與預測值擬合效果均較好，相關(guān)系數(shù)R2分別為0.784和0.857．擬合曲線(藍色線)與1∶1標準直線相差不大．通過將原始數(shù)據(jù)觀測值和交叉驗證所得到的預測值進行分析比較，發(fā)現(xiàn)在考慮節(jié)點最大深度為主要因素時，通過Kriging高斯模型插值，P=7 a預測節(jié)點積水深度均值為1.395 m，實測節(jié)點積水深度均值1.396 m；P=50 a預測節(jié)點積水深度均值為1.655 m，實測節(jié)點積水深度均值1.654 m．研究結(jié)果表明，在考慮節(jié)點最大深度因素的Kriging插值模擬高斯模型預測可取得較好的擬合效果．

圖6 交叉驗證結(jié)果Fig.6 The result of cross-validation

2.2.4 節(jié)點積水深度空間分布預測分析

選用高斯Kriging插值模型對研究區(qū)域P=7 a和P=50 a節(jié)點積水深度數(shù)據(jù)進行空間插值模擬(圖7)．從圖7可以看出，兩個重現(xiàn)期下節(jié)點積水深度變化范圍基本一致，分別為0.627～7.577 m與0.639～8.100 m．在空間上呈現(xiàn)東北與西南向積水分布較良好，東南與西北向積水分布較嚴重，這可能是高強度降雨下，東北與與西南方向的管網(wǎng)布置較稀疏且均是管徑較小的管段，而東南與西北方向的管網(wǎng)布置較密集且均是管徑較大的干管及排放口排放高重現(xiàn)期形成的即時徑流量累積導致的．此外，圖8中西北方向積水深度紅色區(qū)域逐漸擴大，且綠色區(qū)域逐漸被黃色區(qū)域占據(jù)，開始向紅色區(qū)域過渡，呈現(xiàn)出綠色區(qū)域先快后慢的減少，而紅色區(qū)域先快后慢的擴大趨勢．這很可能是因為高重現(xiàn)期早期較大強度降雨使管徑偏小的管段達到滿載負荷狀態(tài)，使得研究區(qū)域地勢較低容易積水的地方出現(xiàn)內(nèi)澇．

圖7 節(jié)點積水深度Kriging高斯模型插值結(jié)果Fig.7 Interpolation results of Kriging Gaussian model for node water depth

3 結(jié) 論

1)通過對蘭州市某校區(qū)雨水管網(wǎng)P=7 a和P=50 a節(jié)點積水深度數(shù)據(jù)進行Pearson相關(guān)性分析，結(jié)果表明：與節(jié)點內(nèi)底標高對節(jié)點積水深度的影響程度相比，節(jié)點最大深度對節(jié)點的積水深度影響更大．其Pearson相關(guān)系數(shù)分別為0.605(P=7 a)和0.766(P=50 a)，從而確定節(jié)點最大深度作為積水深度插值模擬的主要影響因子．

2)針對原始積水點數(shù)據(jù)構(gòu)建穩(wěn)定、球形、高斯和指數(shù)模型進行半變異函數(shù)模型的優(yōu)選及空間差異統(tǒng)計分析，發(fā)現(xiàn)在兩種降雨情境下得到的積水深度數(shù)據(jù)中，模型曲線與實際值平均化擬合度較好的均為高斯模型．在模型參數(shù)分析中，高斯模型的偏差均值較其他模型最接近0(P=7 a時KRME=-0.87×10-4；P=50 a時KRME=0.87×10-3)，其一致性系數(shù)最接近1(P=7 a時，KRMSE=0.939；P=50 a時，KRMSE=0.947)．

3)通過Kriging高斯模型插值，兩個重現(xiàn)期下節(jié)點積水深度變化基本處于0.627～8.100 m內(nèi)，且呈現(xiàn)東北與西南向積水分布較良好，東南與西北向積水分布較嚴重．同時積水深度紅色區(qū)域逐漸擴大，綠色區(qū)域逐漸被黃色占據(jù)，開始向紅色區(qū)域過渡．

4)在雨水管網(wǎng)中，節(jié)點(雨水井)積水內(nèi)澇除了與節(jié)點內(nèi)底高程和最大深度等節(jié)點本身的結(jié)構(gòu)屬性變量有關(guān)外，可能還與大氣濕度、氣溫、環(huán)境等因素密切相關(guān)．因此，為使得插值結(jié)果更加精確化，可考慮適當引入上述環(huán)境因子進行空間插值．

5)通過對常規(guī)開發(fā)情景下的研究區(qū)域進行節(jié)點積水模擬，得到不同降雨重現(xiàn)期下節(jié)點積水情況，結(jié)果表明局部區(qū)域存在內(nèi)澇風險．對于這種局部內(nèi)澇問題，國內(nèi)傳統(tǒng)的思想是“以排為主”，即采取加大管道管徑，提高排放能力的措施來解決．但隨著新型雨洪管理理念的引入，低影響開發(fā)理念越來越多地應(yīng)用在城市建設(shè)中．本研究為LID措施的開展提供一定的理論基礎(chǔ)，可以通過對比不同的LID措施對于研究區(qū)域雨洪控制的影響，篩選出最佳的LID措施組合來降低局部內(nèi)澇的風險．