泰森多邊形距離反比加權(quán)法在BTOPMC模型中的應(yīng)用分析

馬亞麗，敖天其，王漢濤

(1.甘肅農(nóng)業(yè)大學(xué) 水利水電工程學(xué)院，甘肅 蘭州 730070；2.四川大學(xué) 水力學(xué)與山區(qū)河流開(kāi)發(fā)保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610065；3.中國(guó)長(zhǎng)江三峽集團(tuán)公司，湖北 宜昌 443133)

我國(guó)監(jiān)測(cè)的雨量以點(diǎn)雨量為主，產(chǎn)流計(jì)算過(guò)程要求將點(diǎn)雨量轉(zhuǎn)化為面雨量，為此，通過(guò)對(duì)雨量數(shù)據(jù)進(jìn)行插值得到每一個(gè)網(wǎng)格計(jì)算單元的雨量值。目前，并不存在某一個(gè)降水插值模型能夠適應(yīng)諸多分布式水文模型，當(dāng)然這也是水文學(xué)界的一個(gè)重要難題，因此針對(duì)某一個(gè)特定的水文模型，探討某種雨量插值方法對(duì)其模擬精度以及參數(shù)的影響變化規(guī)律是有一定意義且十分必要的，不僅可以提高模型模擬的效果，而且更有利于我們深入地認(rèn)識(shí)該水文模型。朱會(huì)義等[1]采用反距離加權(quán)插值法、克里金插值法、樣條函數(shù)插值法等，根據(jù)潮白河流域58個(gè)雨量站1990年的降雨觀測(cè)數(shù)據(jù)，得到插值結(jié)果的不確定性與雨量站個(gè)數(shù)成反比、與時(shí)間尺度成反比。陳艷英等[2]用樣條函數(shù)法、克里金內(nèi)插法及距離平方反比法，對(duì)溫州地區(qū)1991—2001年90個(gè)觀測(cè)站的月平均降雨進(jìn)行雨量插值，結(jié)果表明各種方法效果均較好，距離平方反比法精度誤差較小。覃建明等[3]提出了基于泰森多邊形法法估算流溪河模型單元網(wǎng)格面雨量的方法，通過(guò)分析六場(chǎng)實(shí)測(cè)洪水過(guò)程的模擬效果，說(shuō)明泰森多邊形法適用于流溪河模型洪水預(yù)報(bào)降雨估算的需要。李金潔等[4]以西南地區(qū)1996—2000年93個(gè)氣象臺(tái)站觀測(cè)的月均降雨量為基礎(chǔ)，采用反距離加權(quán)法(IDW)和普通克里金(O-Kriging)兩種方法進(jìn)行空間插值，通過(guò)交叉驗(yàn)證結(jié)果對(duì)兩種方法進(jìn)行分析比對(duì)，結(jié)果表明：在降雨量的空間插值中，由于研究區(qū)域和時(shí)間尺度的不同，并不存在絕對(duì)的最優(yōu)方法，應(yīng)根據(jù)實(shí)際應(yīng)用效果選擇最適合的方法。汪青靜等[5]采用4種空間插值算法對(duì)研究區(qū)2006—2014年的日降雨進(jìn)行插值，并將面均雨量作為新安江模型的輸入，分析和比較其降雨徑流響應(yīng)，建議選用計(jì)算簡(jiǎn)便的插值算法計(jì)算面雨量，比如泰森多邊形、反距離權(quán)重法。本文為了探究泰森多邊形的距離反比加權(quán)法與BTOPMC水文模型自帶的最鄰近插值法對(duì)洪水模擬效果的影響差異，利用泰森多邊形的距離反比加權(quán)法對(duì)BTOPMC水文模型降雨空間插值方法進(jìn)行改進(jìn)，并對(duì)兩者的模擬效果進(jìn)行比較分析，以期提高BTOPMC模型洪水模擬的可靠性，從而進(jìn)一步改進(jìn)和完善BTOPMC模型，為將來(lái)深入研究該模型和廣泛應(yīng)用該模型奠定一定的基礎(chǔ)。

1 BTOPMC模型簡(jiǎn)介

BTOPMC(Block-Wise use of TOPMODEL with Muskingum-Cunge Method)模型是一個(gè)基于物理機(jī)制的分布式流域水文模型[6- 7]，具體需要率定的參數(shù)個(gè)數(shù)少、且全部具有物理意義等顯著優(yōu)點(diǎn)。BTOPMC模型已在國(guó)內(nèi)外得到一定程度的推廣和應(yīng)用，在日本的富士川流域，黃河流域[8]、印尼的 kali Brantas 流域、尼泊爾的8個(gè)流域[9]以及美國(guó)的十幾個(gè)中小型流域取得了較好的驗(yàn)證。模型采用TOPMODEL作為基本的產(chǎn)流子模型、Muskingum-Cunge(馬斯京根—康奇法) 作為匯流子模型，并加入了由敖天其等開(kāi)發(fā)的洼地處理、河網(wǎng)生成和自然子流域劃分等有關(guān)地形解析的子模型[10]。BTOPMC主要子模型組成結(jié)構(gòu)如圖1所示。目前模型主要有5個(gè)參數(shù)：飽和土壤導(dǎo)水率T0(m2/s)、飽和土壤導(dǎo)水率的衰減因子m(m)、根層最大儲(chǔ)水量Srmax(m)、平均土壤飽和差初始值Sbar0(m)、曼寧糙率系數(shù)n0。其中Srmax主要反映植被/土地利用的影響，T0和m反映流域土壤種類的影響，Sbar0主要反映流域地形的影響，n0同時(shí)反映土壤種類和植被/土地利用的影響。其中，前4個(gè)為產(chǎn)流子模型中的參數(shù)，n0為匯流子模型參數(shù)，這些參數(shù)都具有明確的物理意義。Srmax、T0、m、n0是以每一個(gè)DEM網(wǎng)格來(lái)率定的，而Sbar0(k)是以每一個(gè)子流域來(lái)率定的。

圖1 BTOPMC主要子模型組成結(jié)構(gòu)

BTOPMC模型采用了兩個(gè)精度評(píng)價(jià)指標(biāo)進(jìn)行模型評(píng)價(jià)：Vr和E。

Vr是總模擬流量(TVsim)與實(shí)際資料觀測(cè)值(TVobs)之比，公式如下：

Vr=TVsim/TVobs

(1)

Nash效率E[11- 12]，計(jì)算公式如下：

(2)

式中，N—時(shí)間步長(zhǎng)的總數(shù)；Qsim—時(shí)間步長(zhǎng)的模擬流量；Qobs—時(shí)間步長(zhǎng)的觀測(cè)流量；Qav—是整個(gè)時(shí)期的觀測(cè)流量平均值。

Nash值越高則表明模擬效果越好，100%代表模擬結(jié)果與實(shí)測(cè)資料完全一樣。

除了以上兩個(gè)評(píng)價(jià)指標(biāo)，可以根據(jù)具體需要，采用模擬和觀測(cè)洪峰時(shí)間及平均流量等其他參數(shù)作為輔助的評(píng)價(jià)指標(biāo)。

2 泰森多邊形的距離反比加權(quán)法與BTOPMC模型的耦合

2.1 泰森多邊形的距離反比加權(quán)法

距離反比加權(quán)法(IDW)也叫反距離加權(quán)法，是一種常用的加權(quán)平均插值方法，由Shepard[13]最早提出。該法將計(jì)算區(qū)域劃分成若干矩形網(wǎng)格，每個(gè)網(wǎng)格的寬度和長(zhǎng)度分別為Δx和Δy，各網(wǎng)格點(diǎn)的雨量Zj由周圍鄰近的雨量站實(shí)測(cè)資料按距離平方的倒數(shù)插值求得[14]。每個(gè)DEM網(wǎng)格點(diǎn)與所有已知的雨量站點(diǎn)的權(quán)重之和為1，各雨量站的貢獻(xiàn)值(即權(quán)重)隨該DEM網(wǎng)格點(diǎn)與各雨量站的距離增大而減少，公式如下：

(3)

式中：Zj—第j個(gè)待插值DEM網(wǎng)格點(diǎn)的雨量值，(j=1，2，3，…，n)；Zi—第i個(gè)已知雨量站點(diǎn)的雨量值，(i=1，2，3，…，n)；n—流域內(nèi)已知雨量站總數(shù)；Di—待插值DEM網(wǎng)格點(diǎn)與第i個(gè)雨量站點(diǎn)的距離；p—距離的冪，距離反比加權(quán)法中p=2。

泰森多邊形法是一種極端的邊界內(nèi)插方法，將各雨量站兩兩相連并作連線的中垂線，中垂線相交所形成得若干個(gè)多邊形即為泰森多邊形。以泰森多邊形內(nèi)僅有的一個(gè)氣象站的雨量值來(lái)代表這個(gè)多邊形區(qū)域內(nèi)的降雨量。利用泰森多邊形法將流域劃分為多個(gè)子流域，對(duì)每個(gè)泰森多邊形內(nèi)的所有已經(jīng)利用距離反比加權(quán)法插值的DEM網(wǎng)格點(diǎn)的雨量求平均值，即為泰森多邊形的距離反比加權(quán)法。以每個(gè)泰森多邊形為單位，對(duì)經(jīng)過(guò)距離反比加權(quán)法插值的DEM網(wǎng)格點(diǎn)雨量值進(jìn)行雨量平均，計(jì)算公式如下：

(4)

式中，pk—第k個(gè)泰森多邊形內(nèi)的的平均雨量，(k=1，2，3，…，n)；Zj，k—第k個(gè)泰森多邊形內(nèi)第j個(gè)DEM網(wǎng)格的雨量值，j=1，2，3，…，n；nk—第k個(gè)泰森多邊形內(nèi)的DEM網(wǎng)格個(gè)數(shù)，k=1，2，3，…，n；k—泰森多邊形的編號(hào)；j—DEM網(wǎng)格的編號(hào)。

泰森多邊形的距離反比加權(quán)法結(jié)合了距離反比加權(quán)法和泰森多邊形法兩種降雨插值方法，使各雨量站對(duì)流域內(nèi)所有DEM網(wǎng)格點(diǎn)的雨量值都有貢獻(xiàn)，卻又不會(huì)影響B(tài)TOPMC模型的結(jié)構(gòu)，便于與BTOPMC模型的耦合。

2.2 與BTOPMC模型的耦合

在分布式水文模型BTOPMC中，敖天其等采用的降雨插值方法為最鄰近插值法，即每個(gè)DEM網(wǎng)格點(diǎn)的雨量值采用距離該DEM網(wǎng)格點(diǎn)最近的那個(gè)雨量站的雨量值，并且將每個(gè)DEM網(wǎng)格點(diǎn)最近的雨量站編號(hào)存儲(chǔ)在以后綴“*.sit”的文件中，結(jié)合存儲(chǔ)在“*.pri”文件中對(duì)應(yīng)雨量站的雨量值，從而可以得到每個(gè)DEM網(wǎng)格點(diǎn)的雨量值，值得注意的是在“*.pri”文件中每一列對(duì)應(yīng)一個(gè)已知雨量站點(diǎn)各時(shí)刻的雨量，每列的編號(hào)與雨量站的編號(hào)一致。

本文為實(shí)現(xiàn)泰森多邊形的距離反比加權(quán)法與BTOPMC模型的耦合，首先利用泰森多邊形法將流域劃分為多個(gè)子流域并依次編號(hào)，每個(gè)泰森多邊形內(nèi)的DEM網(wǎng)格點(diǎn)均采用該泰森多邊形的編號(hào)，DEM網(wǎng)格點(diǎn)編碼信息存儲(chǔ)在后綴為“*.sit”文件，通過(guò)距離反比加權(quán)法對(duì)各DEM網(wǎng)格點(diǎn)的雨量值進(jìn)行插值，每個(gè)泰森多邊形內(nèi)內(nèi)所有DEM網(wǎng)格的雨量值逐時(shí)刻求平均值代表該子流域的雨量，并存儲(chǔ)在后綴“*.pri”文件中，泰森多邊形的編號(hào)與“*.pri”文件的列編碼相對(duì)應(yīng)。

3 改進(jìn)的BTOPMC模型在明江流域的應(yīng)用

3.1 明江流域及資料概況

渠江(渠河、南江河)是嘉陵江的一條支流，發(fā)源于川、陜兩省交界處的米倉(cāng)山，流經(jīng)南江、渠縣、廣安、鄰水等縣境，在合川市匯入嘉陵江，并最終流入長(zhǎng)江。總流域面積3.92萬(wàn)km2，多年平均流量約為663m3/s，多年平均徑流量約為235億m3，河道主干流長(zhǎng)約305.5km，天然落差532m，干流平均比降0.17‰。明江為渠江的二級(jí)支流，河長(zhǎng)93km，流域面積1204km2，河口流量25.6m3/s，總落差1155m。渠江流域大部分屬亞熱帶濕潤(rùn)季風(fēng)氣候，年降水量在1014～1500mm以上，雨季主要集中在7—9月份，占全年降雨量的50%以上。本文研究區(qū)域位于明江上游趕場(chǎng)(二)水文站以上，控制流域面積約273km2，多年平均流量7.54m3/s，水位變幅5.9m。明江水系圖如圖2所示。

圖2 明江流域水系圖

BTOPMC模型徑流模擬所需數(shù)據(jù)資料見(jiàn)表1。該研究范圍內(nèi)收集了3個(gè)雨量站和1個(gè)水文站資料，由中國(guó)氣象局成都高原氣象研究所提供：趕場(chǎng)(二)、貴民、沙壩雨量站；趕場(chǎng)(二)水文站。水文站與雨量站基本情況見(jiàn)表2。

表1 模型所需數(shù)據(jù)資料

表2 水文站及雨量站基本情況表

3.2 參數(shù)率定和驗(yàn)證的結(jié)果與分析

本文采用趕場(chǎng)(二)(2007/7/2 1∶00～2007/7/8 0∶00)的洪水?dāng)?shù)據(jù)進(jìn)行BTOPMC模型參數(shù)的率定；利用趕場(chǎng)(二)(2004/9/29 1∶00～2004/10/3 0∶00)的洪水?dāng)?shù)據(jù)進(jìn)行BTOPMC模型參數(shù)的驗(yàn)證。根據(jù)BTOPMC模型精度評(píng)價(jià)指標(biāo)要求，結(jié)合本文洪水模擬的特點(diǎn)，選用以下4個(gè)評(píng)價(jià)指標(biāo)對(duì)模擬結(jié)果進(jìn)行評(píng)價(jià)，包括：Nash效率(E)，模擬總徑流量與實(shí)測(cè)總徑流量之比(Vr)、模擬洪峰時(shí)刻與實(shí)測(cè)洪峰時(shí)刻之差(ΔQt)、模擬洪峰洪量與實(shí)測(cè)洪峰洪量的相對(duì)誤差(ΔQmax)，率定期和驗(yàn)證期的的評(píng)價(jià)指標(biāo)值見(jiàn)表3。

通過(guò)模擬結(jié)果比較分析，可以看出，在參數(shù)率定過(guò)程中，模型自帶的最鄰近點(diǎn)插值法Nash效率高于泰森多邊形的距離反比加權(quán)法Nash效率，但是在驗(yàn)證過(guò)程中，泰森多邊形的距離反比加權(quán)法Nash效率高于最鄰近點(diǎn)插值法Nash效率，同時(shí)，泰森多邊形的距離反比加權(quán)法模擬洪峰時(shí)刻與實(shí)測(cè)洪峰時(shí)刻相吻合，對(duì)洪峰時(shí)刻的模擬效果更準(zhǔn)確，模擬洪峰流量與實(shí)測(cè)洪峰流量相對(duì)誤差更小，可見(jiàn)，在更為注重對(duì)洪峰時(shí)刻及洪峰流量的模擬效果的洪水模擬中，泰森多邊形的距離反比加權(quán)法優(yōu)于模型自帶的雨量插值方法，對(duì)BTOPMC模型降雨插值方法的改進(jìn)，有利于提高模型對(duì)于洪水模擬的精度。

表3 明江流域水文斷面參數(shù)率定和驗(yàn)證評(píng)價(jià)指標(biāo)值

模型自帶的最鄰近點(diǎn)插值法與泰森多邊形的距離反比加權(quán)法兩種雨量插值方法對(duì)洪水過(guò)程的形態(tài)變化影響不同，相比之下，泰森多邊形的距離反比加權(quán)法對(duì)洪峰流量、洪峰時(shí)刻的模擬效果更佳。率定期洪水過(guò)程線如圖3—4所示，驗(yàn)證期洪水過(guò)程線如圖5—6所示。

圖3 率定期洪水過(guò)程線(模型自帶最鄰近點(diǎn)插值法)

圖4 率定期洪水過(guò)程線(泰森多邊形的距離反比加權(quán)法)

圖5 驗(yàn)證期洪水過(guò)程線(模型自帶最鄰近點(diǎn)插值法)

圖6 驗(yàn)證期洪水過(guò)程線(泰森多邊形的距離反比加權(quán)法)

BTOPMC模型主要有5個(gè)參數(shù)，具有需要率定的參數(shù)個(gè)數(shù)少、且全部具有物理意義等顯著優(yōu)點(diǎn)，并且模型采用SCE-UA優(yōu)化算法對(duì)參數(shù)進(jìn)行率定并

表4 明江流域水文斷面模擬參數(shù)自動(dòng)優(yōu)化值

確定最終參數(shù)值，模擬參數(shù)自動(dòng)優(yōu)化值具體見(jiàn)表4。因目前模型并不具有參數(shù)敏感性分析模塊，本文通過(guò)增加模型參數(shù)的初始值和上下限，并以此為參照，與兩種降雨插值方法自動(dòng)優(yōu)化確定的參數(shù)進(jìn)行比較，依據(jù)兩種方法的參數(shù)變化幅度，找出變化最大的參數(shù)值，以相對(duì)誤差η的平方ε表示，ε保留整數(shù)部分，公式3如下：

(5)

式中，εi—相對(duì)于初始值對(duì)應(yīng)參數(shù)的相對(duì)誤差平方，ε1最鄰近點(diǎn)插值法；ε2泰森多邊形的距離反比加權(quán)法；x0—模擬參數(shù)初始值；xi—兩種降雨插值法模擬參數(shù)自動(dòng)優(yōu)化值。

采用最鄰近點(diǎn)插值法，模型參數(shù)曼寧糙率系數(shù)(建筑區(qū))受雨量插值方法的影響最大，Sbar0(21≤λ≤30、31≤λ≤40)、m值(粘土)受雨量插值方法的影響較大；采用泰森多邊形的距離反比加權(quán)法，模型參數(shù)m值(粉土)、Sbar0(31≤λ≤40)受雨量插值方法的影響最大，m值(粘土)、Sbar0(21≤λ≤30)受雨量插值方法的影響較大，這說(shuō)明了BTOPMC模型的各參數(shù)是具有一定物理意義的。

4 結(jié)論

本文采用泰森多邊形的距離反比加權(quán)法對(duì)BTOPMC模型的降雨空間插值方法進(jìn)行改進(jìn)，改變模型自帶的最鄰近點(diǎn)插值法僅僅考慮最近雨量站的雨量代表DEM網(wǎng)格點(diǎn)雨量值的局限性，而采用泰森多邊形的距離反比加權(quán)法的雨量插值方法，使得所有已知雨量站均對(duì)某DEM網(wǎng)格點(diǎn)雨量值有貢獻(xiàn)，通過(guò)比較兩種降雨量插值方法在明江流域洪水過(guò)程的模擬結(jié)果，可以看出，在更注重對(duì)洪峰時(shí)刻及洪峰流量的模擬效果的洪水模擬中，泰森多邊形的距離反比加權(quán)法優(yōu)于模型自帶的雨量差值方法，通過(guò)對(duì)BTOPMC模型降雨插值方法的改進(jìn)，有利于提高BTOPMC模型洪水模擬精度，對(duì)洪峰流量、洪峰時(shí)刻有較好的模擬；采用最鄰近點(diǎn)插值法，模型參數(shù)n0(建筑區(qū))受雨量插值方法的影響最大，Sbar0(21≤λ≤30、31≤λ≤40)、m值(粘土)受雨量插值方法的影響較大；采用泰森多邊形的距離反比加權(quán)法，模型參數(shù)m值(粉土)、Sbar0(31≤λ≤40)受雨量插值方法的影響最大，m值(粘土)、Sbar0(21≤λ≤30)受雨量插值方法的影響較大，這說(shuō)明BTOPMC模型的各參數(shù)是具有一定物理意義的。

本文利用泰森多邊形的距離反比加權(quán)法對(duì)BTOPMC模型的降雨空間插值方法進(jìn)行改進(jìn)，并取得一定成果，但是雨量插值方法在不同時(shí)間尺度上對(duì)模型的影響、不同的DEM分辨率下雨量插值方法對(duì)BTOPMC模型的影響等還需要進(jìn)一步探討。