徐 晗，徐建剛

（1.南京大學(xué) 建筑與城市規(guī)劃學(xué)院，江蘇 南京 210093；2.上海營邑城市規(guī)劃設(shè)計(jì)股份有限公司，上海 200030）

通過近十年來的不斷探索，中國海綿城市規(guī)劃工作開始從流域入手，從流域整體排洪防澇來進(jìn)行規(guī)劃設(shè)計(jì)。為了更科學(xué)地確定海綿城市大中型海綿體的選址、滯水容量與淹沒范圍，需要引入水文模型對河道斷面在暴雨情景下通過的流量進(jìn)行模擬與預(yù)測。伴隨著空間信息科學(xué)與技術(shù)的進(jìn)步，基于流域柵格化技術(shù)的雨洪時(shí)空模擬方法逐漸發(fā)展起來。運(yùn)用數(shù)學(xué)與物理方程對水文循環(huán)中的各個(gè)過程進(jìn)行解釋與描述，通過參數(shù)與變量來表達(dá)其空間分布上的差異性[1]。這種模擬現(xiàn)實(shí)物理法則的物理模型能夠更準(zhǔn)確地描述水文循環(huán)過程，并且具有很強(qiáng)的適應(yīng)性[2]。因而更適用于對空間精度要求較高，空間差異性敏感的城市規(guī)劃領(lǐng)域。本文在林蔚[3]等提出的暴雨匯流水文模型集成分析的基礎(chǔ)上，結(jié)合海綿城市規(guī)劃的實(shí)際設(shè)計(jì)工作需求，面向海綿城市的應(yīng)用進(jìn)行優(yōu)化。通過對水文學(xué)、水力學(xué)相關(guān)模型的研究，基于多層次GIS空間分析模型技術(shù)，系統(tǒng)地建立了從流域雨洪過程相關(guān)影響因子?xùn)鸥窕幚怼⒌疆a(chǎn)匯流過程模型的集成化模擬方法，首次實(shí)現(xiàn)了對流域任意柵格單元內(nèi)河道斷面的徑流過程線的可視化測算。同時(shí)，還對于模型中單元水流長度與河道單元流速計(jì)算方法進(jìn)行了優(yōu)化改進(jìn)。并以閩西國家歷史文化名城長汀縣城汀江觀音橋水文站上游流域?yàn)閷?shí)驗(yàn)區(qū)，使用水文站實(shí)測數(shù)據(jù)對模型進(jìn)行對比驗(yàn)證，實(shí)現(xiàn)了通過對多個(gè)情景的模型參數(shù)可視化調(diào)整來提升模擬精度的技術(shù)應(yīng)用。

1 分時(shí)段產(chǎn)流過程模擬

在流域徑流形成過程中，流域內(nèi)的降水經(jīng)過植被截留、下滲、蒸發(fā)等損失而產(chǎn)生凈雨過程即為產(chǎn)流過程。目前國內(nèi)主要采用的產(chǎn)流模擬大多基于GIS平臺(tái)對流域空間數(shù)據(jù)進(jìn)行柵格化表達(dá)，進(jìn)而構(gòu)建了美國農(nóng)業(yè)部水土保持局在上世紀(jì)50年代提出的SCS-CN流域水文模型，通過對暴雨過程中各個(gè)時(shí)段流域空間產(chǎn)流量的累積計(jì)算，以獲得流域分時(shí)段產(chǎn)流量柵格數(shù)據(jù)。SCS-CN模型綜合考慮了降雨量、土壤透水性、土地利用類型、前期土壤濕度等因子與流域產(chǎn)流的關(guān)系，

通過徑流曲線數(shù)（CN）來推求最大可能滯留量（S），CN是一個(gè)經(jīng)驗(yàn)性無因次參數(shù)，是對降雨前流域下墊面土壤類型、土地利用方式、前期土壤含水量所造成影響的一個(gè)綜合反映。

公式如下：

式中，S為最大可能滯留量（mm）；CN為徑流曲線數(shù)。

美國農(nóng)業(yè)部土壤保持局在分析了大量長期的實(shí)驗(yàn)結(jié)果基礎(chǔ)上，經(jīng)過一系列推導(dǎo)得到以下計(jì)算產(chǎn)流量公式：

式中，R為實(shí)際產(chǎn)流量（mm）；P為降雨量（mm）；S為最大可能滯留量（mm）。

基于上述模型原理與國內(nèi)相關(guān)研究成果，本研究結(jié)合實(shí)際海綿城市規(guī)劃項(xiàng)目中所能獲得的相關(guān)資料，對產(chǎn)流量的影響因子進(jìn)行梳理，并運(yùn)用空間插值、地圖代數(shù)等GIS技術(shù)實(shí)現(xiàn)了參數(shù)的全部柵格化，最終形成產(chǎn)流過程模擬計(jì)算的系統(tǒng)性多層次概念模型（圖1）。

圖1 產(chǎn)流過程模擬概念模型

1.1 降雨量（P）柵格化

根據(jù)圖1流程所示，降雨量指標(biāo)值的空間分布式計(jì)算是模擬的基礎(chǔ)。降雨資料通常為散布于流域內(nèi)的多個(gè)雨量監(jiān)測站點(diǎn)位數(shù)據(jù)，可利用IDW插值、泰森多邊形等方法柵格化為覆蓋全流域的降雨量柵格數(shù)據(jù)。

在實(shí)際應(yīng)用中，雨量監(jiān)測站通常記錄的是各單位時(shí)段內(nèi)的降雨量（mm）。由于SCS-CN模型中對產(chǎn)流量計(jì)算公式是基于一場完整降雨，公式（2）中的初損量0.2 s不會(huì)隨著降雨時(shí)長變化。這使得在計(jì)算分時(shí)段降雨所形成的徑流時(shí)需要將每個(gè)小時(shí)降雨都看作一場完整降雨，并分別運(yùn)用公式（2）進(jìn)行計(jì)算。則有，對于歷時(shí)m小時(shí)的降雨場次α，其中第n小時(shí)（1≤n≤m）的累計(jì)降雨量（Pn）為n小時(shí)降雨量之和，公式如下：式中，Pn為第n小時(shí)降雨的累計(jì)降雨量（mm）；pi為第i小時(shí)降雨的降雨量（mm）。

1.2 徑流曲線數(shù)（CN）柵格化

徑流曲線數(shù)CN值在流域空間的柵格化推算是模擬的關(guān)鍵，本研究從流域下墊面三大因子特征量化、坡降影響修正和土壤前期濕度分級三方面探析，建立了適用于實(shí)驗(yàn)區(qū)的參數(shù)及其科學(xué)計(jì)算方法。首先，通過查閱美國農(nóng)業(yè)部在網(wǎng)絡(luò)發(fā)布的說明文件《National Engineering Handbook Hydrology Chapters》中的表格，可以獲取到在不同土地利用類型、土壤透水度與植被覆蓋度組合下的CN值。土地利用類型因子主要為流域土地利用現(xiàn)狀圖；土壤透水性因子則為流域土壤的HSG分級的空間分布數(shù)據(jù)，其包括了A/B/C/D 4個(gè)級別的指標(biāo)，透水性由A至D依次降低；植被覆蓋率基于多波段衛(wèi)星影像計(jì)算獲得的歸一化植被指數(shù)（NDVI）柵格數(shù)據(jù)，根據(jù)楊勝天[5]等的研究內(nèi)容計(jì)算獲得，并依據(jù)流域所在地的植被特征，重分類為好、中、差3個(gè)等級的柵格數(shù)據(jù)。最后將三大因子的空間數(shù)據(jù)運(yùn)用疊加標(biāo)識(shí)、柵格重分類等GIS技術(shù)，進(jìn)而生成平均土壤濕度狀態(tài)下的流域CN值柵格數(shù)據(jù)。

其次，坡降作為水文學(xué)中影響徑流的重要因子，在地形變化較大的山地丘陵地區(qū)應(yīng)用SCS-CN模型時(shí)還應(yīng)考慮坡度影響并對徑流的影響。本研究采用陳正維[6]的研究成果，對CN值進(jìn)行修正：

式中，CN為徑流曲線數(shù)，slp為坡降值。

最后，由于土壤前期濕度會(huì)顯著影響土壤對降雨的初損吸收量，因此，SCS-CN模型以當(dāng)次降雨前5 d的雨量為基準(zhǔn)來判斷土壤前期濕度，其將土壤前期濕度分為AMCI，AMCII，AMCIII三級指標(biāo)，分別代表干、平均、濕3種狀態(tài)。本研究采用羅鵬[7]的研究成果對AMC進(jìn)行修正，以上述經(jīng)過坡度修正的CN值作為平均狀態(tài)，其他2種狀態(tài)的情況則根據(jù)下面2個(gè)公式進(jìn)行轉(zhuǎn)換：

式中，CNI為干狀態(tài)下的CN值；CNII為平均狀態(tài)下的CN值；CNIII為濕狀態(tài)下的CN值。

最終獲得3種狀態(tài)下的流域CN值柵格數(shù)據(jù)，根據(jù)暴雨情景前5 d雨量選擇使用。

1.3 分時(shí)段產(chǎn)流量計(jì)算

基于地圖代數(shù)GIS技術(shù)，將前述所獲得的Pn柵格化成果及CN值柵格代入SCS-CN模型公式（1）、（2）中作為參數(shù)，計(jì)算得出降雨n小時(shí)后流域的累計(jì)產(chǎn)流量（Qn），而降雨在第n小時(shí)的產(chǎn)流量（qn）則為第n場降雨累計(jì)產(chǎn)流量（Qn）與第n-1場累計(jì)產(chǎn)流量（Qn-1）之差：

式中，qn為第n小時(shí)降雨的分時(shí)段產(chǎn)流量（mm）；Qn為第n小時(shí)降雨的累計(jì)產(chǎn)流量（mm）。計(jì)算獲得流域一場降雨中的分時(shí)段產(chǎn)流量柵格。

2 匯流過程模擬

匯流過程，即流域凈雨流入河道中的聚集過程。流域匯流過程模擬采用空間分布式旅行法模型[8]，模型通過追蹤各單元徑流匯流至流域出口所經(jīng)過路徑，并應(yīng)用曼寧公式、連續(xù)方程等水力學(xué)公式計(jì)算獲得徑流通過各單元所花費(fèi)的時(shí)間成本，進(jìn)而累加計(jì)算各流域單元到達(dá)流域出口所要花費(fèi)的匯流時(shí)間，生成流域等流時(shí)面柵格，即所產(chǎn)生徑流同時(shí)達(dá)到出口斷面的柵格單元數(shù)目，定義為該時(shí)刻的無因次流量。最終結(jié)合產(chǎn)流模型所生成的分時(shí)段流域產(chǎn)流量柵格，匯總統(tǒng)計(jì)生成流域出口斷面流量過程線。這一匯流過程計(jì)算模式的概念層次模型如下（圖2）。

圖2 匯流過程模擬概念模型

2.1 流域等流時(shí)面柵格計(jì)算

降雨徑流在流域進(jìn)行匯流過程時(shí)存在坡面漫流與河道流2種不同的匯流模式，需要將流域劃分為坡面單元與河道單元，計(jì)算流域單元水流長度并分別計(jì)算各單元的匯流成本，進(jìn)而生成等流時(shí)面。因此，計(jì)算流域各單元的徑流通過時(shí)間成本需要?jiǎng)澐譃?個(gè)計(jì)算步驟。

2.1.1 坡面單元與河道單元?jiǎng)澐钟?jì)算

本研究引入河道提取閾值來界定坡面與河道兩類單元。該域值被定義為形成河道匯流所需要的最小集水面積占流域總面積的比值?；谟闪饔駾EM生成的流域流向柵格，可以通過向流動(dòng)方向累加的方法生成的無因次流量柵格（累計(jì)數(shù)為柵格像元個(gè)數(shù)，圖3），則河道單元即是像元值大于“河道提取閾值*流域柵格總數(shù)”的柵格像元，其余像元便是坡面單元。由于降雨量的不同會(huì)導(dǎo)致最小集水面積的變化，這使得河道提取閾值需要通過實(shí)驗(yàn)對比實(shí)現(xiàn)河流情況獲得。

2.1.2 單元水流長度計(jì)算

傳統(tǒng)單元水流長度的確定，通常依據(jù)柵格像元大小計(jì)算，當(dāng)流向?yàn)檎媳毕蚧蛘龞|西向時(shí)，取柵格像元邊長大小；當(dāng)流向?yàn)樾睂蔷€方向時(shí)，取柵格像元斜對角線長度（柵格像元邊長×20.5）。這種方法存在2個(gè)問題：①在現(xiàn)實(shí)中水流通常沿地形成曲線匯流，這種方法會(huì)導(dǎo)致分析用柵格像元越大，所獲取值的徑流長度相比現(xiàn)實(shí)徑流長度越小。②這種方法所獲取到的僅是路徑像元的平面長度，而實(shí)際水流是在路徑單元的坡面進(jìn)行匯流，忽略了路徑坡度的影響，坡度越大計(jì)算結(jié)果相較現(xiàn)實(shí)徑流長度越小。為解決以上2個(gè)問題，研究提出以下優(yōu)化方法。

1）地圖精度糾正。如圖4所示，不同比例尺地圖下同一條匯水線的形態(tài)與長度并不相同，受比例尺精度影響，比例尺越小，誤差越大。根據(jù)國內(nèi)學(xué)者對數(shù)字地圖曲線長度精度的相關(guān)研究[9-10]，本研究采用公式（8）進(jìn)行徑流長度糾正。

圖4 基于不同比例尺地形圖所生成的匯水線差異

式中，ld為修正后徑流長度；li為糾正前徑流長度；d為地圖精度糾正系數(shù)。

通過對河流進(jìn)行不同柵格單元大小（設(shè)單元邊長為r）的柵格化處理，并統(tǒng)計(jì)獲得各個(gè)柵格的柵格數(shù)N，將測量結(jié)果做log(N)～log(r)雙對數(shù)散點(diǎn)圖，擬合得一負(fù)斜率直線，該斜率的絕對值為分維值D。作者通過對研究流域內(nèi)多條主要河流進(jìn)行上述計(jì)算，獲得的分維值D約為1.05～1.11，綜合分析實(shí)驗(yàn)與前人研究的結(jié)果，本研究取1.1為實(shí)證流域單位水流長度的地圖精度糾正系數(shù)。

2）坡度糾正。在計(jì)算路徑像元長度時(shí)，可乘以該像元坡度的正割值，進(jìn)行坡度糾正。可轉(zhuǎn)化為與坡降值（即像元坡度的正切值）相關(guān)的公式：

式中，lp為坡度修正后路徑像元徑流長度；li為路徑像元徑流長度；slp為路徑像元坡降值

2.1.3 坡面單元通過成本計(jì)算

坡面單元通過成本時(shí)間由穩(wěn)態(tài)運(yùn)動(dòng)波近似值及曼寧公式計(jì)算獲得,計(jì)算公式如下[11]：

式中，to為坡面單元時(shí)間成本（s）；L為單元水流長度（m）；n為曼寧糙率；ie為平均凈雨強(qiáng)度（m/s）；slp為坡降值。

其中平均凈雨強(qiáng)度柵格ie表示由累計(jì)產(chǎn)流量柵格除以降雨時(shí)長獲得：

式中，ie為平均凈雨強(qiáng)度（m/s）；Rn為第n時(shí)段降雨的累計(jì)產(chǎn)流量（m）；T為降雨時(shí)長（s）。

坡降值柵格由流域DEM柵格，為每個(gè)像元計(jì)算從該像元到與其相鄰的像元方向上的最大變化率，并輸出像元值為坡降值的坡降柵格。這樣計(jì)算獲得的非常小的甚至為0的坡降值，這使得部分坡面單元的流速將接近于0，與現(xiàn)實(shí)徑流情況不符。研究通過設(shè)定最小坡降值，去除柵格中過小的坡降值，通常設(shè)值為0.01。

曼寧糙率柵格基于流域土地利用柵格生成，通過查閱相關(guān)文獻(xiàn)標(biāo)準(zhǔn)，賦予對應(yīng)土地利用糙率數(shù)值。

2.1.4 河道單元通過成本計(jì)算

河道單元流速公式由寬河道的穩(wěn)態(tài)連續(xù)方程結(jié)合曼寧公式轉(zhuǎn)化而來，計(jì)算公式如下[11]

式中，tc為河道單元時(shí)間成本（m/s）；B為河道的有效寬度（m）；n為曼寧糙率；L為單元水流長度（m）；Qe為平均累計(jì)匯流量（m3/s）；sl為坡降值。

其中對于河道有效寬度的計(jì)算，由于現(xiàn)實(shí)中河道寬度與其通過流量成正相關(guān)，則對于無因次流量柵格中的河道像元而言，其匯積值越大，該像元的河道有效寬度越大。從而就可以通過計(jì)算流域出口斷面的河道有效寬度，運(yùn)用河道的無因次流量柵格推求河道有效寬度柵格。公式如下：

式中，Bi為當(dāng)前像元河道有效寬度（m）；Bm為流域出口斷面河道有效寬度（m）；Qi為當(dāng)前像元無因次流量；Qm為流域出口無因次流量。

河道平均累計(jì)匯流量柵格計(jì)算類似無因次流量，采用累計(jì)流量法，并采用平均凈雨強(qiáng)度柵格作為權(quán)重柵格計(jì)算獲得。

對于河道曼寧糙率值，由于地形圖及土地利用圖等資料難以真實(shí)表達(dá)河道底部情況，故在計(jì)算時(shí)均采用固定數(shù)值。河道坡降取值與坡面坡降一致，取坡降柵格值。

2.1.5 匯流時(shí)間柵格生成

基于流域流向柵格可追蹤獲得各柵格像元流向流域出口的唯一匯流路徑，結(jié)合匯流成本柵格，就可以通過累加計(jì)算得到該像元產(chǎn)生徑流流至流域出口所花費(fèi)的總時(shí)間，將時(shí)間作為該像元的新柵格值，便生成了匯流時(shí)間柵格。公式如下：

式中，Tm為流域內(nèi)某一柵格像元匯流至流域出口的時(shí)間；m為匯流路徑上的路徑像元數(shù)量；τi為徑流通過某一路徑像元所花費(fèi)的時(shí)間。

在計(jì)算生成流域匯流時(shí)間柵格后，進(jìn)而通過對時(shí)間柵格值按單位時(shí)段長度進(jìn)行分類，通常與降雨數(shù)據(jù)的分時(shí)段長度一致，從而獲得流域等流時(shí)面柵格（圖5），表達(dá)在同一時(shí)段內(nèi)匯流至流域出口流域柵格像元。

圖5 左：匯流時(shí)間柵格；右：等流時(shí)面柵格

2.2 流域出口斷面流量過程線計(jì)算

以產(chǎn)流模型所計(jì)算獲得的分時(shí)段產(chǎn)流量柵格數(shù)據(jù)作為權(quán)重，對流域等流時(shí)面柵格進(jìn)行分區(qū)統(tǒng)計(jì)，統(tǒng)計(jì)各時(shí)段等流時(shí)面內(nèi)對應(yīng)產(chǎn)流量像元值的總和，進(jìn)而可生成該時(shí)段產(chǎn)流在流域出口所形成的流量過程線。最后應(yīng)用傳統(tǒng)單位線法的疊加假定，即出口斷面的流量過程線等于各個(gè)時(shí)段流量過程錯(cuò)開時(shí)段疊加之和，對分時(shí)段流量過程線錯(cuò)開時(shí)段進(jìn)行疊加，最終生成流域出口斷面流量過程線。

3 模擬實(shí)證

本研究以長汀縣觀音橋水文站控制斷面以上的汀江上游流域?yàn)槟M實(shí)證應(yīng)用示范區(qū)，流域面積約為382 km2。以林蔚[3]研究成果為對照組進(jìn)行方法優(yōu)化設(shè)計(jì)與結(jié)果精度驗(yàn)證，選擇觀音橋水文站20 150 519場次洪水，對應(yīng)重現(xiàn)期十年一遇。該場次為2012－2015年間最大洪水場次，且林蔚[3]研究中的模型效率系數(shù)最低、峰現(xiàn)誤差時(shí)間最高。結(jié)果如表1所示,實(shí)驗(yàn)設(shè)定在其他參數(shù)及數(shù)據(jù)一致的條件下進(jìn)行。其中對照組不進(jìn)行任何優(yōu)化，實(shí)驗(yàn)組01僅對單元水流長度進(jìn)行優(yōu)化計(jì)算；實(shí)驗(yàn)組02僅對河道有效寬度進(jìn)行優(yōu)化計(jì)算；實(shí)驗(yàn)組03對二者均進(jìn)行優(yōu)化計(jì)算，各實(shí)驗(yàn)組形成流域出口斷面過程線對比圖見圖6。

表1 模型精度對照表

圖6 實(shí)驗(yàn)流域出口斷面流量過程線對比

根據(jù)實(shí)驗(yàn)結(jié)果可以看出，2種優(yōu)化措施效果均較為明顯，采用優(yōu)化后的模擬方法得到的結(jié)果在模型精度驗(yàn)證的4項(xiàng)指標(biāo)中均有提升，模型等級由丙級（0.70≥E＞0.50）提升為乙級（0.90≥E＞0.70）。按照城市總體規(guī)劃一般設(shè)定的20～30年期限，該結(jié)果可以作為海綿城市規(guī)劃中的大型海綿體（濕地公園等）控制性指標(biāo)依據(jù)。

4 結(jié)論與展望

隨著我國海綿城市建設(shè)向縱深方向發(fā)展,海綿城市規(guī)劃建設(shè)方案的實(shí)施效果評價(jià)成為關(guān)鍵，規(guī)劃建設(shè)方案是否科學(xué)、可行和高效將對海綿城市建設(shè)能否真正提升城市韌性產(chǎn)生深遠(yuǎn)影響。

本研究旨在為缺少水文觀測數(shù)據(jù)的小微流域的海綿體規(guī)劃提供科學(xué)的水文模擬數(shù)據(jù)，在優(yōu)化改進(jìn)部分因子的計(jì)算獲取方法的同時(shí)，模擬方法中還存在沒有明確計(jì)算公式取值需要通過率定的參數(shù)。未來可進(jìn)一步深入研究，借助人工智能技術(shù)，采用神經(jīng)網(wǎng)絡(luò)、深度學(xué)習(xí)等算法技術(shù)，實(shí)現(xiàn)計(jì)算機(jī)自動(dòng)求得最優(yōu)參數(shù)。