鄭州“7·20”暴雨移植子牙河流域方法及成果

□楊豐源

2021年汛期，鄭州發(fā)生了歷史罕見的“7·20”特大暴雨，造成了重大人員傷亡和財(cái)產(chǎn)損失。許多地區(qū)、部門以期通過鄭州暴雨在本地區(qū)重演，來檢驗(yàn)防范措施是否有效，應(yīng)急預(yù)案是否完備可行，達(dá)到補(bǔ)齊短板、降低損失的目的。現(xiàn)以子牙河流域?yàn)槔?，從天氣形?shì)上分析目標(biāo)流域發(fā)生鄭州“7·20”暴雨的可能性，然后分析暴雨中心位置，進(jìn)而介紹移植技術(shù)、移植結(jié)果，為子牙河流域防范鄭州“7·20”級(jí)別特大暴雨提供參考。

1.移植可能性分析

歷 史上1939年、1956年、1963年、1996年暴雨，均在子牙河系出現(xiàn)大洪水，天氣形勢(shì)與此次相近。從形成原因來看，1939年夏天赤道附近系統(tǒng)活躍，副高位置偏西北，臺(tái)風(fēng)北上次數(shù)多，海河流域汛期暴雨次數(shù)多；1956年受臺(tái)風(fēng)外圍及副高影響，輸送充足水汽；1963年暴雨成因?yàn)槎啻蔚蛪翰酆偷蜏u共同影響結(jié)果；1996年暴雨成因?yàn)?608號(hào)臺(tái)風(fēng)北上減弱為低氣壓與西風(fēng)槽和冷空氣結(jié)合的影響?！?·20”暴雨為集合南亞高壓東伸、副熱帶高壓西伸北抬、環(huán)流形勢(shì)阻塞、低渦緩慢西移、強(qiáng)臺(tái)風(fēng)遠(yuǎn)程輸送水汽、低空急流發(fā)展等多種形勢(shì)和因素共同影響。因此，從天氣形勢(shì)上分析，子牙河系存在發(fā)生“7·20”暴雨的可能性；從地形特征來看，子牙河系西部太行山迎風(fēng)坡抬升特性明顯，“7·20”暴雨鄭州市西北側(cè)山地起到了相同的作用，因此將鄭州“7·20”暴雨移植至子牙河系是可行的。

2.暴雨中心選擇

結(jié)合典型歷史暴雨洪水的暴雨中心位置，如1939年滏陽河一帶，1956年和1963年暴雨中心分別發(fā)生在滹沱河上游井陘縣和滏陽河獐犭厶，1996年暴雨中心在野溝門水庫，2016年暴雨中心分別發(fā)生在邢臺(tái)臨城一帶、石家莊贊皇、石家莊井陘一帶，可以看出暴雨中心基本分布在太行山迎風(fēng)坡，滹沱河主要以冶河中上游為中心，滏陽河主要以臨城水庫、朱莊水庫等為中心。

綜上分析，若將“7·20”鄭州暴雨移植至子牙河流域，那么暴雨中心應(yīng)在滏陽河朱莊、臨城水庫上游或者滹沱河冶河中上游，現(xiàn)以暴雨中心在冶河微水上游和泜河臨城水庫上游進(jìn)行移植分析。

3.移植技術(shù)應(yīng)用

基于數(shù)字高程，提取流域邊界，使用反距離權(quán)重內(nèi)插將平移后點(diǎn)雨量值內(nèi)插至網(wǎng)格，再計(jì)算流域面雨量。

3.1 填洼處理

在依據(jù)數(shù)字高程圈畫流域中，流域高低起伏。計(jì)算水流方向時(shí)，存在有些水不能流出洼地，使得水系產(chǎn)生誤差或不連續(xù)。為了使水系在遇到洼地流向依然連續(xù)，需對(duì)DEM數(shù)據(jù)進(jìn)行洼地處理。填洼前基于數(shù)字化流域電子水系，采用AGREE算法對(duì)原始DEM數(shù)據(jù)進(jìn)行預(yù)處理，使DEM和輸入的矢量圖層相一致。先設(shè)置緩沖區(qū)、平滑、增益3個(gè)參數(shù)。河道及周圍區(qū)域地形呈“V”形，以確保水最終流入河流網(wǎng)格。根據(jù)設(shè)置的增益參數(shù)，河道所在網(wǎng)格高程將再次降低相應(yīng)值，以便根據(jù)矢量河道校正DEM。AGREE算法修正后流域數(shù)字高程模型（DEM）GTOPO30見圖1。

圖1 AGREE算法修正后流域數(shù)字高程模型（DEM）GTOPO30

填洼處理的目的是將洼地內(nèi)部高程填至與出口高程一致，按照增加中心柵格高程方法填平洼地，使處理后的DEM能生成收斂連續(xù)的河道。

3.2 網(wǎng)格流向計(jì)算

在填洼后DEM基礎(chǔ)上計(jì)算網(wǎng)格流向，使用D8單流向算法：根據(jù)公式（1）計(jì)算每個(gè)網(wǎng)格單元與周圍網(wǎng)格的坡度，其中為中心網(wǎng)格與相鄰網(wǎng)格中心之間的高程差。按最陡坡度原則設(shè)定該單元水流流向，流向用一系列數(shù)字代表。D8算法提取流向見圖2。

圖2 D8算法提取流向

3.3 集水面積確定

在確定流向后，可根據(jù)流向分析集水面積。若上游柵格的水量都匯集到某一柵格點(diǎn)，便可據(jù)此計(jì)算該柵格的累計(jì)水量。流域出口處，集水面積最大，流域內(nèi)的水最終都流向它。根據(jù)圖2中DEM確定的集水面積矩陣見圖3。

圖3 根據(jù)圖2中DEM確定的集水面積矩陣

3.4 河道閾值確定及流域提取

劃清河道及水系需搞清楚形成永久性水道所必須的集水面積，集水面積大于等于集水面積閾值的柵格被定義為河道，反之則為坡地。有些研究認(rèn)為絕大多數(shù)侵蝕性斜坡上凸下凹，并且在兩部分間有一個(gè)轉(zhuǎn)折，轉(zhuǎn)折出現(xiàn)的臨界區(qū)域確定了河網(wǎng)的界限。在具有足夠侵蝕介質(zhì)的區(qū)域（例如具有良好土壤覆蓋的流域），侵蝕能力（輸沙率）f可表示為：

式中：

s—坡度；

q—流量；

k—侵蝕擴(kuò)散率；

m、n—系數(shù)。

坡地的侵蝕主要由坡度決定，河道的侵蝕除坡度外也與流量關(guān)系密切。對(duì)于坡地n=0、m=1，對(duì)于河道n通常大于1、m大致為2。Smith和Bretherton進(jìn)一步研究認(rèn)為形成凹形坡的條件為：

反之，則容易形成一個(gè)凸斜率。對(duì)于流量Q和流域A，通常滿足Q∝Ax，指標(biāo)x與流量回歸周期有關(guān)。當(dāng)考慮較大的回歸周期(如年流量)時(shí)，x=1。Krikby即以單位等高線的集水面積a取代式（2）和（3）中q，得到：

式中：

F—輸沙量；

β、u、w—系數(shù)。

(5)式建立時(shí)，容易形成凹坡，反之則容易形成凸坡。經(jīng)過進(jìn)一步的研究，Tarboton得出結(jié)論，在當(dāng)?shù)貤l件下有一個(gè)固定的輸沙速率。在此假設(shè)下，輸沙量只與流域面積有關(guān):

U為常數(shù)，（6）式兩邊分別對(duì)a求導(dǎo)并都乘以a得到公式（7）進(jìn)而得到公式（8）：

（8）式右邊與（5）式形式相近，因而公式（9）成立時(shí)易形成凹形坡，a(?F/?s)為正值，故凹形坡滿足ds/da<0；反之，凸形坡則為ds/da>0。

基于上述理論，閾值通常由集水區(qū)面積和河道平均坡度之間的關(guān)系得出。首先給出集水區(qū)的取值范圍，計(jì)算每個(gè)數(shù)字排水系統(tǒng)平均坡度，繪制平均坡度與流域面積關(guān)系曲線。平均坡度與流域面積關(guān)系曲線拐點(diǎn)對(duì)應(yīng)的流域面積，可視為流域內(nèi)河流地貌發(fā)育的關(guān)鍵支撐區(qū)域，即流域面積閾值。

通過計(jì)算分析，結(jié)合數(shù)字化矢量水系，推求出流域集水面積閾值為50柵格。確定集水面積閾值后，集水面積大于等于50的柵格被定義為河道，其余均為坡地。海河流域數(shù)字水系見圖4。

圖4 海河流域數(shù)字水系

3.5 流域確定

確定分水嶺邊界的算法如下：確定閾值后，搜索DEM矩陣以找到分水嶺面積大于或等于閾值的網(wǎng)格，網(wǎng)格是流域的出口（DEM矩陣可以包含多個(gè)流域），按順序編號(hào)（n=1,2,3，…）。對(duì)于任何編號(hào)的出口網(wǎng)格，需根據(jù)流量矩陣表向上搜索，任何流向它的網(wǎng)格都分配給這個(gè)一編號(hào)的網(wǎng)格，并繼續(xù)循環(huán)搜索，直到再也沒有流向該網(wǎng)格為止。水系矩陣（a）與流域矩陣（b）見圖5。

圖5 水系矩陣（a）與流域矩陣（b）

3.6 控制站分區(qū)

控制站分區(qū)一般按水文站控制，先劃分最下游站控制的集水區(qū)，再依次劃分上游的集水區(qū)，這樣上游區(qū)會(huì)覆蓋在原來的下游區(qū)中。各區(qū)域在生成過程中，以流量控制站的柵格為該集水區(qū)成員，標(biāo)上區(qū)號(hào)，不斷循環(huán)直到所有柵格被標(biāo)上區(qū)號(hào)。

3.7 降水空間分布及控制站分區(qū)面雨量計(jì)算

空間分布，也稱為空間插值。根據(jù)是否考慮高程影響，降水的空間分布可分為僅考慮平面位置影響的二維空間插值、僅考慮高程影響的降雨—高程線性回歸和同時(shí)考慮平面位置和高程影響的三維空間插值?？紤]到數(shù)據(jù)的局限性，此文采用二維空間插值，只考慮插值。

泰森多邊形網(wǎng)格和泰森多邊形雨量測(cè)量站之間的距離是流域內(nèi)所有雨量測(cè)量站中最短的。然而，根據(jù)泰森多邊形劃分降水量將導(dǎo)致降水量不均勻。在泰森多邊形的內(nèi)部，降雨是一致的，泰森多邊形之間的降雨可能會(huì)存在較大差異。

與泰森多邊形法不同，反距離加權(quán)平均法在降雨連續(xù)性方面更好。思路就是網(wǎng)格降雨與雨量站所在的網(wǎng)格降雨成正比，與從網(wǎng)格到每個(gè)雨量測(cè)量站所在網(wǎng)格的距離成反比。計(jì)算公式如下：

式(10)（11）中：

P—期望格網(wǎng)的降雨值；

N—雨測(cè)站個(gè)數(shù)；

P（si）—第i個(gè)雨測(cè)站的降雨值

λ1—第i個(gè)雨測(cè)站的權(quán)重

d1—第i個(gè)雨測(cè)站與期望格網(wǎng)的距離；

b—權(quán) 重指標(biāo)，b=0時(shí)，式(10)演化為算術(shù)平均算法;b=1時(shí)，式(10)為倒數(shù)距離比法;b=2時(shí)，式(10)為距離平方的倒數(shù)比法，此文使用了廣泛使用的距離平方倒數(shù)比法。

待控制站分區(qū)內(nèi)逐柵格雨量插值后，計(jì)算控制站分區(qū)內(nèi)所有柵格雨量時(shí)間序列的算數(shù)平均值，可得到控制站分區(qū)的面雨量時(shí)間序列。距離反比權(quán)重插值見圖6。

圖6 距離反比權(quán)重插值

4.移植結(jié)果分析

尖崗水庫7月18日—22日5天累積雨量950.8mm，其中7月20日15～16時(shí)最大1h雨量147mm。暴雨中心移植至不同位置時(shí)主要區(qū)間面雨量統(tǒng)計(jì)見表1。移植降雨量如下：方案一移植至冶河微水位置，子牙河流域面雨量184.3mm，其中滹沱河崗南水庫以上241.6mm，崗黃區(qū)間448mm，滏陽河山區(qū)產(chǎn)流區(qū)214.5mm，暴雨中心位于冶河中游區(qū)域。方案二移植至滏陽河臨城水庫，子牙河流域面雨量202.9mm，其中滹沱河崗南水庫以上110.8mm，崗黃區(qū)間382.5mm，滏陽河山區(qū)產(chǎn)流區(qū)339.9mm，暴雨中心位于滏陽河山區(qū)水庫上游。通過對(duì)比，兩種方案僅滹沱河降雨較“63·8”小，方案一冶河降水略大于“96·8”和“16·7”降水，滏陽河山區(qū)降水小于“96·8”和“16·7”降水。方案二冶河降水小于“96·8”和“16·7”降水，滏陽河山區(qū)降水接近“96·8”降水，大于“16·7”降水?？傮w上，移植降水與“96·8”降水量級(jí)大致相當(dāng)。□

表1 暴雨中心移植至不同位置時(shí)主要區(qū)間面雨量統(tǒng)計(jì)