安航永，趙文龍，張雅萍

（1.廣州禺山水務(wù)勘測(cè)設(shè)計(jì)股份有限公司，廣東 廣州 511400；2.黃河水利水電開(kāi)發(fā)總公司，河南 鄭州 450000；3.黃河建工集團(tuán)有限公司，河南 鄭州 450000）

調(diào)洪演算是水庫(kù)規(guī)劃設(shè)計(jì)中的重要內(nèi)容之一，其目的是根據(jù)已知的水位～庫(kù)容關(guān)系、洪水過(guò)程、調(diào)洪調(diào)度方案和起調(diào)水位，演算并輸出各計(jì)算時(shí)段末的庫(kù)水位、泄流量和庫(kù)容變化量等[1]。調(diào)洪演算計(jì)算方法主要有試算法、半圖解法、數(shù)值解法等，其中半圖解法是最為簡(jiǎn)單、有效的演算方法[2]。但是，隨著計(jì)算機(jī)的普及，適合編程的數(shù)值解法得到快速發(fā)展，半圖解法逐漸在實(shí)用中淡出[3]。筆者借助MATLAB強(qiáng)大的數(shù)據(jù)處理和數(shù)據(jù)插值能力以及Excel友好的人機(jī)交互能力，程序化半圖解法，使其能夠更加高效、便捷的解決調(diào)洪演算問(wèn)題。

1 MATLAB在數(shù)據(jù)插值方面的應(yīng)用

MATLAB內(nèi)置多種數(shù)據(jù)插值函數(shù)，可以進(jìn)行一維數(shù)據(jù)插值（interp1）、二維數(shù)據(jù)插值（interp2）和三維數(shù)據(jù)插值（interp3）等。在工程計(jì)算中，較常用的插值函數(shù)為interp1，其運(yùn)用的一般格式為：y=interp1（X，Y，x，method），其中 method 為插值方法[4]。常用的插值方法有：nearest（最近鄰點(diǎn)插值）、linear（線性插值）、spline（三次樣條函數(shù)插值）、pchip（分段三次Hermite插值）等，函數(shù)默認(rèn)插值方法為linear。

2 半圖解法基本原理

水庫(kù)調(diào)洪演算的基本原理是求解水庫(kù)的水量平衡方程[5]：

半圖解法在求解過(guò)程中將水量平衡方程改寫(xiě)為[6]：

式中：Q1、Q2分別為計(jì)算時(shí)段初、末的入庫(kù)流量，m3/s；q1、q2分別為計(jì)算時(shí)段初、末的下泄流量，m3/s；V1、V2分別為計(jì)算時(shí)段初、末水庫(kù)的蓄水量，m3；Δt計(jì)算時(shí)段，s。

半圖解法計(jì)算步驟如下：

（1）根據(jù)水位～庫(kù)容關(guān)系、水位～泄流關(guān)系以及計(jì)算時(shí)段等繪制輔助曲線，見(jiàn)圖1。

（2）確定起調(diào)水位 Z1，及相應(yīng)的q1、V1，計(jì)算各時(shí)段平均入庫(kù)流量Qp。

（3）在圖1水位坐標(biāo)軸上確定Z1位置，記為a點(diǎn)。作水平線ac交 A線于b點(diǎn)，使 bc=Qp。因A 線為（V/Δt-q/2）=f1(Z)，則 ab 等 于 （V1/Δt-q1/2），ac 等 于 Qp+（V1/Δt-q1/2）=（V1/Δt+q1/2）。

（4）在c點(diǎn)做垂線交B線于d點(diǎn)，由d點(diǎn)作水平線de交Z坐標(biāo)軸于 e點(diǎn)，可見(jiàn) de=ac=（V2/Δt+q2/2）。因 B線為（V/Δt+q/2）=f2(Z)，d 點(diǎn)位于 B 線上，則 e點(diǎn)為 Z2值。

（5）過(guò)de與C線交點(diǎn)f作垂線交q坐標(biāo)軸于g點(diǎn)，則g點(diǎn)為q2值。

（6）根據(jù)Z2值，查水位～庫(kù)容曲線可得V2。

（7）將 e點(diǎn)的 Z2值作為第二時(shí)段的 Z1，重復(fù)（2）～（6）即可得下一時(shí)段的特征值。由此逐時(shí)段進(jìn)行計(jì)算，即可完成全部計(jì)算。

圖1 半圖解法輔助曲線

3 計(jì)算程序及使用方法

本方法由Excel完成基礎(chǔ)數(shù)據(jù)的錄入工作，然后通過(guò)MATLAB程序調(diào)用基礎(chǔ)數(shù)據(jù)并完成調(diào)洪演算。

（1）數(shù)據(jù)錄入

首先，根據(jù)程序要求創(chuàng)建名稱(chēng)為data的Excel文件。在data.xls內(nèi)建立5個(gè)工作表，分別用來(lái)存儲(chǔ)洪水過(guò)程（RuLiu）、水位 ～ 庫(kù)容關(guān)系（KuRong）、水位 ～ 泄流關(guān)系（XieLiu）、起算條件（Set）以及演算結(jié)果（OutPut）。各工作表數(shù)據(jù)錄入格式見(jiàn)圖2～圖5。

圖2 洪水過(guò)程錄入格式

圖3 水位～庫(kù)容關(guān)系錄入格式

圖4 水位～泄流關(guān)系錄入格式

該程序通過(guò)MATLAB中xlsread函數(shù)自動(dòng)調(diào)用data.xls中基礎(chǔ)數(shù)據(jù)，采用interp1完成半圖解法中的插值運(yùn)算，然后通過(guò)xlswrite函數(shù)將演算結(jié)果輸出至“OutPut”工作表，整個(gè)計(jì)算過(guò)程均由計(jì)算機(jī)完成，即提高了演算速度又避免了人工繪圖和讀數(shù)的誤差。

4 應(yīng)用實(shí)例

以官涌水庫(kù)調(diào)洪演算為例，驗(yàn)證該方法的方便適用性。官涌水庫(kù)位于廣州市番禺區(qū)，為?。?）型水庫(kù)，水庫(kù)大壩為均質(zhì)土壩，壩高9.53 m，壩頂高程19.97 m，設(shè)計(jì)洪峰流量為2.93 m3/s，校核洪峰流量為4.39 m3/s。官涌水庫(kù)泄流建筑物為一座1.5 m×1.55 m（寬×高）開(kāi)敞式溢洪道，溢洪道底高程16.82 m。

官涌水庫(kù)設(shè)計(jì)洪水過(guò)程、水位～庫(kù)容關(guān)系、水位～泄流關(guān)系見(jiàn)圖6～圖8，水庫(kù)調(diào)度原則為敞泄，即起調(diào)水位為16.82 m。

圖6 設(shè)計(jì)洪水過(guò)程線

圖7 水位～庫(kù)容關(guān)系曲線

圖8 水位～泄流關(guān)系曲線

根據(jù)官涌水庫(kù)設(shè)計(jì)洪水過(guò)程、水位～庫(kù)容關(guān)系、水位～泄流關(guān)系及水庫(kù)調(diào)度原則，在錄入基礎(chǔ)數(shù)據(jù)后運(yùn)行MATLAB程序完成計(jì)算。程序輸出結(jié)果包括入庫(kù)流量（m3/s）、出庫(kù)流量（m3/s）、時(shí)段蓄量（萬(wàn) m3）、庫(kù)容（萬(wàn) m3）和水位（m），演算結(jié)果見(jiàn)圖9～圖10，結(jié)果輸出格式如表1。

圖9 調(diào)洪演算結(jié)果（出/入庫(kù)流量）

圖10 調(diào)洪演算結(jié)果（庫(kù)水位）

表1 演算結(jié)果輸出格式

5 結(jié)語(yǔ)

本文以官涌水庫(kù)調(diào)洪演算為例，聯(lián)合運(yùn)用MATLAB程序和Excel，方便、快速的完成了官涌水庫(kù)調(diào)洪演算工作。該方法不僅操作方便、計(jì)算迅速，而且避免了人工繪圖和讀數(shù)的誤差，計(jì)算精度也得以提高?？梢?jiàn)，以半圖解法為理論基礎(chǔ)，聯(lián)合運(yùn)用MATLAB程序和Excel解決水庫(kù)調(diào)洪演算問(wèn)題是一次成功的嘗試，不僅解決了半圖解法的程序化問(wèn)題，而且為廣大水利工程技術(shù)人員提供了一個(gè)解決水庫(kù)調(diào)洪演算工作的方便、快捷的方法。