模型簡(jiǎn)介

新安江模型是河海大學(xué)水文院趙人俊教授于上世紀(jì)80年代提出的一個(gè)具有世界影響力的水文模型¹²³⁴。新安江模型是分布式模型，可用于濕潤地區(qū)與半濕潤地區(qū)的濕潤季節(jié)。當(dāng)流域面積較小時(shí)，新安江模型采用集總模型，當(dāng)面積較大時(shí)，采用分塊模型。它把全流域分為許多塊單元流域，對(duì)每個(gè)單元流域作產(chǎn)匯流計(jì)算，得出單元流域的出口流量過程。再進(jìn)行出口以下的河道洪水演算，求得流域出口的流量過程。把每個(gè)單元流域的出流過程相加，就求得了流域的總出流過程。
該模型按照三層蒸散發(fā)模式計(jì)算流域蒸散發(fā)，按蓄滿產(chǎn)流概念計(jì)算降雨產(chǎn)生的總徑流量，采用流域蓄水曲線考慮下墊面不均勻?qū)Ξa(chǎn)流面積變化的影響。在徑流成分劃分方面，對(duì)三水源情況，按“山坡水文學(xué)”產(chǎn)流理論用一個(gè)具有有限容積和測(cè)孔、底孔的自由水蓄水庫把總徑流劃分成飽和地面徑流、壤中水徑流和地下水徑流。在匯流計(jì)算方面，單元面積的地面徑流匯流一般采用單位線法，壤中水徑流和地下水徑流的匯流則采用線性水庫法。河網(wǎng)匯流一般采用分段連續(xù)演算的Muskingum法或滯時(shí)演算法⁵。

源碼下載

該項(xiàng)目已在Github開源，您可以點(diǎn)擊此處下載。

個(gè)人編寫，難免有理解不到位和疏漏的地方，敬請(qǐng)諒解，僅供參考！
感謝導(dǎo)師、大李老師、小李老師、姚老師、小蕾在模型原理、程序設(shè)計(jì)、率定驗(yàn)證等方面的指導(dǎo)！

如有問題，歡迎交流探討！ 郵箱：lujiabo@hhu.edu.cn 盧家波
來信請(qǐng)說明博客標(biāo)題及鏈接，謝謝。

設(shè)計(jì)思路

嚴(yán)格按照“流域分單元，蒸散發(fā)分層次，產(chǎn)流分水源，匯流分階段”的建模流程，以面向?qū)ο蟮姆绞骄帉懭葱掳步Ｐ虲++程序。

將流域分塊、蒸散發(fā)、產(chǎn)流、匯流等模塊封裝成單獨(dú)的C++類，便于移植和復(fù)用。

程序驗(yàn)證

程序模塊

以下僅為主要模塊代碼，完整程序請(qǐng)至Github下載。

流域分塊

Watershed.h

#pragma once #include <string> #include <vector>class IO //從文本中導(dǎo)入降雨和蒸發(fā)數(shù)據(jù)，輸出流量過程到文本中 { public:void ReadFromFile(std::string StrPath = ""); //從文本中導(dǎo)入流域降雨或蒸發(fā)數(shù)據(jù)void WriteToFile(std::string StrPath = ""); //輸出流域出口斷面流量過程到文本中IO();~IO();public:std::vector< std::vector<double> > m_P; //各雨量站逐時(shí)段降雨量，mmstd::vector< std::vector<double> > m_EM; //各蒸發(fā)站逐時(shí)段水面蒸發(fā)量，mmdouble* m_Q; //流域出口斷面流量過程，m3/sint nrows; //數(shù)據(jù)行數(shù)int ncols; //數(shù)據(jù)列數(shù)};//流域分塊 class Watershed { public:void ReadFromFile(std::string StrPath = "");void SetValues(std::string name, double area, int numRainfallStation, int numEvaporationStation, int numSubWatershed,double* areaSubWatershed, std::vector< std::vector<double> > rateRainfallStation,std::vector< std::vector<double> > rateEvaporationStation,std::string* nameRainfallStation, std::string* nameEvaporationStation,std::vector< std::vector<double> > P, std::vector< std::vector<double> > EM);void calculate(const IO * io); //計(jì)算各單元流域降雨量和水面蒸發(fā)量double GetP(int nt, int nw); //得到nt時(shí)刻，第nw單元流域的降雨量，mmdouble GetEM(int nt, int nw); //得到nt時(shí)刻，第nw單元流域的水面蒸發(fā)量，mmdouble GetF(int nw); //得到第nw單元流域的面積，km2int GetnW(); //得到單元流域個(gè)數(shù)Watershed();~Watershed();protected: private:std::string m_name; //流域名稱double m_area; //流域控制面積，km2int m_numRainfallStation; //雨量站個(gè)數(shù)int m_numEvaporationStation; //蒸發(fā)站個(gè)數(shù)int m_numSubWatershed; //單元流域個(gè)數(shù)double* m_areaSubWatershed; //單元流域面積，km2std::vector< std::vector<double> > m_rateRainfallStation; //各單元流域?qū)?yīng)各雨量站比例，按泰森多邊形分配std::vector< std::vector<double> > m_rateEvaporationStation; //各單元流域?qū)?yīng)各蒸發(fā)站比例，按泰森多邊形分配std::string* m_nameRainfallStation; //雨量站點(diǎn)名稱std::string* m_nameEvaporationStation; //蒸發(fā)站點(diǎn)名稱std::vector< std::vector<double> > m_P; //各單元流域逐時(shí)段降雨量，mmstd::vector< std::vector<double> > m_EM; //各單元流域逐時(shí)段水面蒸發(fā)量，mm};

Watershed.cpp

#include "Watershed.h" #include <fstream>void Watershed::ReadFromFile(std::string StrPath) {std::ifstream fin(StrPath + "watershed.txt");fin >> m_name>> m_area>> m_numRainfallStation>> m_numEvaporationStation>> m_numSubWatershed;//讀取單元流域面積m_areaSubWatershed = new double[m_numSubWatershed];for (int i = 0; i < m_numSubWatershed; i++){fin >> m_areaSubWatershed[i];}//讀取雨量站分配比例//創(chuàng)建子流域行*雨量站列的二維動(dòng)態(tài)向量，并初始化為0m_rateRainfallStation = std::vector< std::vector<double> >(m_numSubWatershed, std::vector<double>(m_numRainfallStation, 0.0));for (int i = 0; i < m_numSubWatershed; i++){for (int j = 0; j < m_numRainfallStation; j++){fin >> m_rateRainfallStation[i][j];}}//讀取蒸發(fā)站分配比例//創(chuàng)建子流域行*蒸發(fā)站列的二維動(dòng)態(tài)向量，并初始化為0m_rateEvaporationStation = std::vector< std::vector<double> >(m_numSubWatershed, std::vector<double>(m_numEvaporationStation, 0.0));for (int i = 0; i < m_numSubWatershed; i++){for (int j = 0; j < m_numEvaporationStation; j++){fin >> m_rateEvaporationStation[i][j];}}//讀取雨量站點(diǎn)名稱m_nameRainfallStation = new std::string[m_numRainfallStation];for (int j = 0; j < m_numRainfallStation; j++){fin >> m_nameRainfallStation[j];}//讀取蒸發(fā)站點(diǎn)名稱m_nameEvaporationStation = new std::string[m_numEvaporationStation];for (int j = 0; j < m_numEvaporationStation; j++){fin >> m_nameEvaporationStation[j];}fin.close(); }void Watershed::SetValues(std::string name, double area,int numRainfallStation, int numEvaporationStation, int numSubWatershed,double* areaSubWatershed, std::vector< std::vector<double> > rateRainfallStation,std::vector< std::vector<double> > rateEvaporationStation,std::string* nameRainfallStation, std::string* nameEvaporationStation,std::vector< std::vector<double> > P, std::vector< std::vector<double> > EM) {m_name = name;m_area = area;m_numRainfallStation = numRainfallStation;m_numEvaporationStation = numEvaporationStation;m_numSubWatershed = numSubWatershed;m_areaSubWatershed = areaSubWatershed;m_rateRainfallStation = rateRainfallStation;m_rateEvaporationStation = rateEvaporationStation;m_nameRainfallStation = nameRainfallStation;m_nameEvaporationStation = nameEvaporationStation;m_P = P;m_EM = EM; }void Watershed::calculate(const IO * io) {int nrows = io->nrows; //記錄條數(shù)，即時(shí)段數(shù)int ncols = m_numSubWatershed; //單元流域個(gè)數(shù)//計(jì)算各單元流域逐時(shí)段降雨量，mmm_P = std::vector< std::vector<double> >(nrows, std::vector<double>(ncols, 0.0)); for (int r = 0; r < nrows; r++){for (int c = 0; c < ncols; c++){for (int i = 0; i < m_numRainfallStation; i++){m_P[r][c] += io->m_P[r][i] * m_rateRainfallStation[c][i]; //按比例計(jì)算單元流域降雨量}}}//計(jì)算各單元流域逐時(shí)段水面蒸發(fā)量，mmm_EM = std::vector< std::vector<double> >(nrows, std::vector<double>(ncols, 0.0));for (int r = 0; r < nrows; r++){for (int c = 0; c < ncols; c++){for (int i = 0; i < m_numEvaporationStation; i++){m_EM[r][c] += io->m_EM[r][i] * m_rateEvaporationStation[c][i]; //按比例計(jì)算單元流域水面蒸發(fā)量}}} }double Watershed::GetP(int nt, int nw) {return m_P[nt][nw]; }double Watershed::GetEM(int nt, int nw) {return m_EM[nt][nw]; }double Watershed::GetF(int nw) {return m_areaSubWatershed[nw]; }int Watershed::GetnW() {return m_numSubWatershed; }Watershed::Watershed() {m_name = "默認(rèn)流域";m_area = 0.0;m_numRainfallStation = 0;m_numEvaporationStation = 0;m_numSubWatershed = 0;m_areaSubWatershed = nullptr;m_rateRainfallStation = std::vector< std::vector<double> >(m_numSubWatershed, std::vector<double>(m_numRainfallStation, 0.0));m_rateEvaporationStation = std::vector< std::vector<double> >(m_numSubWatershed, std::vector<double>(m_numEvaporationStation, 0.0));m_nameRainfallStation = nullptr;m_nameEvaporationStation = nullptr;m_P = std::vector< std::vector<double> >(0, std::vector<double>(m_numSubWatershed, 0.0));m_EM = std::vector< std::vector<double> >(0, std::vector<double>(m_numSubWatershed, 0.0)); }Watershed::~Watershed() {delete[] m_areaSubWatershed;}void IO::ReadFromFile(std::string StrPath) {std::ifstream fin(StrPath + "P.txt");//讀入降雨數(shù)據(jù)的行數(shù)和列數(shù)，行數(shù)表示時(shí)段數(shù)，列數(shù)表示雨量站數(shù)fin >> nrows >> ncols;//創(chuàng)建動(dòng)態(tài)二維數(shù)據(jù)，用于存儲(chǔ)各雨量站降雨量，mmm_P = std::vector< std::vector<double> >(nrows, std::vector<double>(ncols, 0.0));//讀取各雨量站降雨量數(shù)據(jù)，mmfor (int r = 0; r < nrows; r++){for (int c = 0; c < ncols; c++){fin >> m_P[r][c];}}fin.close();fin.open(StrPath + "EM.txt");//讀入蒸發(fā)數(shù)據(jù)的行數(shù)和列數(shù)，行數(shù)表示時(shí)段數(shù)，列數(shù)表示蒸發(fā)站數(shù)fin >> nrows >> ncols;//創(chuàng)建動(dòng)態(tài)二維數(shù)據(jù)，用于存儲(chǔ)各蒸發(fā)站蒸發(fā)量，mmm_EM = std::vector< std::vector<double> >(nrows, std::vector<double>(ncols, 0.0));//讀取各蒸發(fā)站蒸發(fā)量數(shù)據(jù)，mmfor (int r = 0; r < nrows; r++){for (int c = 0; c < ncols; c++){fin >> m_EM[r][c];}}fin.close(); }void IO::WriteToFile(std::string StrPath) {//打開Q.txt輸出流域出口斷面流量過程，沒有該文件則新建std::ofstream fout(StrPath + "Q.txt");//輸出流量過程for (int i = 0; i < nrows; i++){fout << m_Q[i] << std::endl;}fout.close(); }IO::IO() {m_Q = nullptr;nrows = 0;ncols = 0;m_P = std::vector< std::vector<double> >(nrows, std::vector<double>(ncols, 0.0));m_EM = std::vector< std::vector<double> >(nrows, std::vector<double>(ncols, 0.0)); }IO::~IO() {delete[] m_Q; }

蒸散發(fā)

Evapotranspiration.h

#pragma once #include "Data.h"//三層蒸散發(fā)模型 class Evapotranspiration { public:void SetParmameter(const Parameter * parameter); //設(shè)置參數(shù)void SetState(const State* state); //設(shè)置狀態(tài)void UpdateState(State* state); //更新狀態(tài)void calculate(); //計(jì)算三層蒸散發(fā)量Evapotranspiration(double kc = 0.8, double lm = 80.0, double c = 0.15,double wu = 10, double wl = 30, double ep = 0.0,double e = 0.0, double eu = 0.0, double el = 0.0, double ed = 0.0, double p = 0.0, double em = 0.0); //構(gòu)造函數(shù)~Evapotranspiration(); //析構(gòu)函數(shù)protected: private://========模型參數(shù)========//double KC; //流域蒸散發(fā)折算系數(shù)，敏感double LM; //下層張力水容量/mm，敏感，60~90double C; //深層蒸散發(fā)折算系數(shù)，不敏感，0.10~0.20//========模型狀態(tài)========//double WU; //上層張力水蓄量，mmdouble WL; //下層張力水蓄量，mmdouble EP; //單元流域蒸發(fā)能力，mmdouble E; //總的蒸散發(fā)量，mmdouble EU; //上層蒸散發(fā)量，mmdouble EL; //下層蒸散發(fā)量，mmdouble ED; //深層蒸散發(fā)量，mm//========外部輸入========//double P; //降雨量，mmdouble EM; //水面蒸發(fā)量，mm};

Evapotranspiration.cpp

#include "Evapotranspiration.h"void Evapotranspiration::SetParmameter(const Parameter* parameter) {KC = parameter->m_KC;LM = parameter->m_LM;C = parameter->m_C; }void Evapotranspiration::SetState(const State* state) {WU = state->m_WU;WL = state->m_WL;P = state->m_P;EM = state->m_EM;}void Evapotranspiration::UpdateState(State* state) {state->m_EP = EP;state->m_E = E;state->m_EU = EU;state->m_EL = EL;state->m_ED = ED; }void Evapotranspiration::calculate() {//三層蒸散發(fā)計(jì)算EP = KC * EM; //計(jì)算流域蒸發(fā)能力if (P + WU >= EP){EU = EP;EL = 0;ED = 0;}else{EU = P + WU;if (WL >= C * LM){EL = (EP - EU) * WL / LM;ED = 0;}else{if (WL >= C * (EP - EU)){EL = C * (EP - EU);ED = 0;}else{EL = WL;ED = C * (EP - EU) - EL;}}}//計(jì)算總的蒸散發(fā)量E = EU + EL + ED; }Evapotranspiration::Evapotranspiration(double kc, double lm, double c, double wu, double wl, double ep,double e, double eu, double el, double ed, double p, double em) {KC = kc;LM = lm;C = c;WU = wu;WL = wl;EP = ep;E = e;EU = eu;EL = el;ED = ed;P = p;EM = em; }Evapotranspiration::~Evapotranspiration() {}

產(chǎn)流

Source.h

#pragma once #include "Data.h"//水源劃分 class Source { public:void SetParmameter(const Parameter* parameter); //設(shè)置參數(shù)void SetState(const State* state); //設(shè)置狀態(tài)void UpdateState(State* state); //更新狀態(tài)void calculate(); //劃分三水源Source(double sm = 20.0, double ex = 1.5, double kg = 0.35, double ki = 0.35,double r = 0.0, double rs = 0.0, double ri = 0.0, double rg = 0.0,double pe = 0.0, double fr = 0.0, double s0 = 0.0, double s = 0.0,double m = 0.0, double kid = 0.0, double kgd = 0.0,double smm = 0.0, double smmf = 0.0, double smf = 0.0, double au = 0.0,double rsd = 0.0, double rid = 0.0, double rgd = 0.0, double fr0 = 0.0,int n = 1, double q = 0.0, double kidd = 0.0, double kgdd = 0.0, double dt_ = 0.0);~Source();protected: private://========模型參數(shù)========//double SM; //表層自由水蓄水容量/mm ，敏感double EX; //表層自由水蓄水容量方次，不敏感，1.0~1.5double KG; //表層自由水蓄水庫對(duì)地下水的日出流系數(shù)，敏感double KI; //表層自由水蓄水庫對(duì)壤中流的日出流系數(shù)，敏感//========模型狀態(tài)========//double R; //總徑流量，mmdouble RS; //地面徑流，mmdouble RI; //壤中流，mmdouble RG; //地下徑流，mmdouble PE; //凈雨量，mmdouble FR; //本時(shí)段產(chǎn)流面積比例double S0; //本時(shí)段初的自由水蓄量，mmdouble S; //本時(shí)段的自由水蓄量，mmdouble M; //一天劃分的計(jì)算時(shí)段數(shù)double KID; //表層自由水蓄水庫對(duì)壤中流的計(jì)算時(shí)段出流系數(shù)，敏感double KGD; //表層自由水蓄水庫對(duì)地下水的計(jì)算時(shí)段出流系數(shù)，敏感double SMM; //全流域單點(diǎn)最大的自由水蓄水容量，mmdouble SMMF; //產(chǎn)流面積上最大一點(diǎn)的自由水蓄水容量，mmdouble SMF; //產(chǎn)流面積上的平均自由水蓄水容量深，mmdouble AU; //相應(yīng)平均蓄水深的最大蓄水深，S0值對(duì)應(yīng)的縱坐標(biāo)，mmdouble RSD; //計(jì)算步長地面徑流，mmdouble RID; //計(jì)算步長壤中流，mmdouble RGD; //計(jì)算步長地下徑流，mmdouble FR0; //上一時(shí)段產(chǎn)流面積比例int N; //N 為計(jì)算時(shí)段分段數(shù)，每一段為計(jì)算步長double Q; //Q 是每個(gè)計(jì)算步長內(nèi)的凈雨量，mmdouble KIDD; //表層自由水蓄水庫對(duì)壤中流的計(jì)算步長出流系數(shù)，敏感double KGDD; //表層自由水蓄水庫對(duì)地下水的計(jì)算步長出流系數(shù)，敏感double dt; //模型計(jì)算時(shí)段長,h };

Source.cpp

#include "Source.h"void Source::SetParmameter(const Parameter* parameter) {SM = parameter->m_SM;EX = parameter->m_EX;KG = parameter->m_KG;KI = parameter->m_KI; }void Source::SetState(const State* state) {R = state->m_R;PE = state->m_PE;FR = state->m_FR;S0 = state->m_S0;dt = state->m_dt; }void Source::UpdateState(State* state) {state->m_RS = RS;state->m_RI = RI;state->m_RG = RG;state->m_S0 = S;state->m_FR = FR;}void Source::calculate() {//出流系數(shù)換算M = 24.0 / dt; //一天劃分的計(jì)算時(shí)段數(shù)KID = (1 - pow(1 - (KI + KG), 1.0 / M)) / (1 + KG / KI); //表層自由水蓄水庫對(duì)壤中流的計(jì)算時(shí)段出流系數(shù)，敏感KGD = KID * KG / KI; //表層自由水蓄水庫對(duì)地下水的計(jì)算時(shí)段出流系數(shù)，敏感//三分水源[4]if (PE <= 1e-5){//凈雨量小于等于0時(shí),這里認(rèn)為凈雨量小于1e-5時(shí)即為小于等于0RS = 0;RI = KID * S0 * FR; /*當(dāng)凈雨量小于等于0時(shí)，消耗自由水蓄水庫中的水產(chǎn)流面積比例仍為上一時(shí)段的比例不變*/RG = KGD * S0 * FR;S = S0 * (1 - KID - KGD); //更新下一時(shí)段初的自由水蓄量}else{//凈雨量大于0時(shí)SMM = SM * (1 + EX); //全流域單點(diǎn)最大的自由水蓄水容量，mmFR0 = FR; //上一時(shí)段產(chǎn)流面積比例FR = R / PE; //計(jì)算本時(shí)段產(chǎn)流面積比例if (FR > 1) //如果FR由于小數(shù)誤差而計(jì)算出大于1的情況，則強(qiáng)制置為1{FR = 1;}S = S0 * FR0 / FR;N = int(PE / 5.0) + 1; //N 為計(jì)算時(shí)段分段數(shù)，每一段為計(jì)算步長Q = PE / N; //Q 是每個(gè)計(jì)算步長內(nèi)的凈雨量，mm//R 是總徑流量，PE是單元流域降雨深度。//把R按5mm分，實(shí)際操作就是把PE按5mm分，是為了減小產(chǎn)流面積變化的差分誤差。//自由水蓄水庫只發(fā)生在產(chǎn)流面積上，其底寬為產(chǎn)流面積FR，顯然FR是隨時(shí)間變化的。//產(chǎn)流量R進(jìn)入水庫即在產(chǎn)流面積上產(chǎn)生PE的徑流深。//FR = f / F = R / PE KIDD = (1 - pow(1 - (KID + KGD), 1.0 / N)) / (1 + KGD / KID); //表層自由水蓄水庫對(duì)壤中流的計(jì)算步長出流系數(shù)，敏感KGDD = KIDD * KGD / KID; //表層自由水蓄水庫對(duì)地下水的計(jì)算步長出流系數(shù)，敏感//把該時(shí)段的RS、RI、RG置0，用于后續(xù)累加計(jì)算步長內(nèi)的RSD、RID、RGDRS = 0.0;RI = 0.0;RG = 0.0;//計(jì)算產(chǎn)流面積上最大一點(diǎn)的自由水蓄水容量 SMMFif (EX == 0.0){SMMF = SMM; //EX等于0時(shí)，流域自由水蓄水容量分布均勻}else{//假定SMMF與產(chǎn)流面積FR及全流域上最大點(diǎn)的自由水蓄水容量SMM仍為拋物線分布//則SMMF應(yīng)該用下式計(jì)算SMMF = (1 - pow(1 - FR, 1.0 / EX)) * SMM;}SMF = SMMF / (1.0 + EX);//將每個(gè)計(jì)算時(shí)段的入流量R分成N段，計(jì)算各個(gè)計(jì)算步長內(nèi)的RSD、RID、RGD，再累加得到RS、RI、RGfor (int i = 1; i <= N; i++){if (S > SMF){S = SMF;}AU = SMMF * (1 - pow(1 - S / SMF, 1.0 / (1 + EX)));if (Q + AU <= 0){RSD = 0;RID = 0;RGD = 0;S = 0;}else{if (Q + AU >= SMMF){//計(jì)算步長內(nèi)的凈雨量進(jìn)入自由水蓄水庫后，使得自由水蓄水深超過產(chǎn)流面積上最大單點(diǎn)自由水蓄水深RSD = (Q + S - SMF) * FR;RID = SMF * KIDD * FR;RGD = SMF * KGDD * FR;S = SMF * (1 - KIDD - KGDD);}else{//自由水蓄水深未超過產(chǎn)流面積上最大單點(diǎn)自由水蓄水深RSD = (S + Q - SMF + SMF * pow(1 - (Q + AU) / SMMF, 1 + EX)) * FR;RID = (S * FR + Q * FR - RSD) * KIDD;RGD = (S * FR + Q * FR - RSD) * KGDD;S = S + Q - (RSD + RID + RGD) / FR;}}//累加計(jì)算時(shí)段內(nèi)的地面徑流、壤中流和地下徑流RS = RS + RSD;RI = RI + RID;RG = RG + RGD;}} }Source::Source(double sm, double ex, double kg, double ki,double r, double rs, double ri, double rg,double pe, double fr, double s0, double s,double m, double kid, double kgd,double smm, double smmf, double smf, double au,double rsd, double rid, double rgd, double fr0,int n, double q, double kidd, double kgdd, double dt_) {SM = sm;EX = ex;KG = kg;KI = ki;R = r;RS = rs;RI = ri;RG = rg;PE = pe;FR = fr;S0 = s0;S = s;M = m;KID = kid;KGD = kgd;SMM = smm;SMMF = smmf;SMF = smf;AU = au;RSD = rsd;RID = rid;RGD = rgd;FR0 = fr0;N = n;Q = q;KIDD = kidd;KGDD = kgdd;dt = dt_; }Source::~Source() {}

單元流域匯流

Confluence.h

#pragma once #include "Data.h"//單元流域匯流，包括坡地匯流和河網(wǎng)匯流，都采用線性水庫法 class Confluence { public:void SetParmameter(const Parameter* parameter); //設(shè)置參數(shù)void SetState(const State* state); //設(shè)置狀態(tài)void UpdateState(State* state); //更新狀態(tài)void calculate(); //坡地匯流和河網(wǎng)匯流Confluence(double cs = 0.1, double ci = 0.6, double cg = 0.95, double cr = 0.2, double im = 0.01,double qs = 0.0, double qi = 0.0, double qg = 0.0, double qt = 0.0, double qu = 0.0,double rs = 0.0, double ri = 0.0, double rg = 0.0, double rim = 0.0,double qi0 = 0.0, double qg0 = 0.0, double qu0 = 0.0, double f = 0.0,double u = 0.0, double m = 0.0, double csd = 0.0, double cid = 0.0, double cgd = 0.0, double crd = 0.0, double dt_ = 24);~Confluence(); protected: private://========模型參數(shù)========//double CS; //地面徑流消退系數(shù)，敏感double CI; //壤中流消退系數(shù)，敏感double CG; //地下水消退系數(shù)，敏感double CR; //河網(wǎng)蓄水消退系數(shù)，敏感double IM; //不透水面積占全流域面積的比例，不敏感//========模型狀態(tài)========//double QS; //單元流域地面徑流，m3/sdouble QI; //單元流域壤中流，m3/sdouble QG; //單元流域地下徑流，m3/sdouble QT; //單元流域河網(wǎng)總?cè)肓?/span>//（進(jìn)入單元面積的地面徑流、壤中流和地下徑流之和），m3/sdouble QU; //單元流域出口流量，m3/sdouble RS; //地面徑流量，mmdouble RI; //壤中流徑流量，mmdouble RG; //地下徑流量，mmdouble RIM; //不透水面積上的產(chǎn)流量，mmdouble QI0; //QI(t-1)，前一時(shí)刻壤中流，m3/sdouble QG0; //QG(t-1)，前一時(shí)刻地下徑流，m3/sdouble QU0; //QU(t-1)，前一時(shí)刻單元流域出口流量，m3/sdouble F; //單元流域面積，km2double U; //單位轉(zhuǎn)換系數(shù)double M; //一天劃分的計(jì)算時(shí)段數(shù)double CSD; //計(jì)算時(shí)段內(nèi)地面徑流蓄水庫的消退系數(shù)double CID; //計(jì)算時(shí)段內(nèi)壤中流蓄水庫的消退系數(shù)double CGD; //計(jì)算時(shí)段內(nèi)地下水蓄水庫的消退系數(shù)double CRD; //計(jì)算時(shí)段內(nèi)河網(wǎng)蓄水消退系數(shù)double dt; //模型計(jì)算時(shí)段長,h };

Confluence.cpp

#include "Confluence.h"void Confluence::SetParmameter(const Parameter* parameter) {CS = parameter->m_CS;CI = parameter->m_CI;CG = parameter->m_CG;CR = parameter->m_CR;IM = parameter->m_IM; }void Confluence::SetState(const State* state) {RIM = state->m_RIM;RS = state->m_RS;RI = state->m_RI;RG = state->m_RG;QS = state->m_QS;QI = state->m_QI;QG = state->m_QG;QU0 = state->m_QU;F = state->m_F;dt = state->m_dt; }void Confluence::UpdateState(State* state) {state->m_QS = QS;state->m_QI = QI;state->m_QG = QG;state->m_QU = QU;state->m_QU0 = QU0; }void Confluence::calculate() {//把透水面積上的產(chǎn)流量均攤到單元流域上RS = RS * (1 - IM);RI = RI * (1 - IM);RG = RG * (1 - IM);//出流系數(shù)換算M = 24.0 / dt; //一天劃分的計(jì)算時(shí)段數(shù)CSD = pow(CS, 1.0 / M); //計(jì)算時(shí)段內(nèi)地面徑流蓄水庫的消退系數(shù)CID = pow(CI, 1.0 / M); //計(jì)算時(shí)段內(nèi)壤中流蓄水庫的消退系數(shù)CGD = pow(CG, 1.0 / M); //計(jì)算時(shí)段內(nèi)地下水蓄水庫的消退系數(shù)CRD = pow(CR, 1.0 / M); //計(jì)算時(shí)段內(nèi)河網(wǎng)蓄水消退系數(shù)//坡地匯流U = F / 3.6 / 24; //單位轉(zhuǎn)換系數(shù)QS = CSD * QS + (1 - CSD) * (RS + RIM) * U; //地面徑流流入地面徑流水庫，//經(jīng)過消退(CSD)，成為地面徑流對(duì)河網(wǎng)的總?cè)肓鱍SQI = CID * QI + (1 - CID) * RI * U; //壤中流流入壤中流水庫，//經(jīng)過消退(CID)，成為壤中流對(duì)河網(wǎng)的總?cè)肓鱍IQG = CGD * QG + (1 - CGD) * RG * U; //地下徑流進(jìn)入地下水蓄水庫，經(jīng)過地下水//蓄水庫的消退(CGD)，成為地下水對(duì)河網(wǎng)的總?cè)肓鱍GQT = QS + QI + QG;//河網(wǎng)匯流，采用線性水庫法，且僅當(dāng)單元流域面積大于200km2時(shí)才計(jì)算河網(wǎng)匯流if (F < 200){QU = QT; //單元流域面積不大且河道較短，對(duì)水流運(yùn)動(dòng)的調(diào)蓄作用通常較小//將這種調(diào)蓄作用合并在地面徑流和地下徑流中一起考慮所帶來的誤差通常可以忽略}else{QU = CRD * QU + (1 - CRD) * QT; //線性水庫法//只有在單元流域面積較大或流域坡面匯流極其復(fù)雜時(shí)//才考慮單元面積內(nèi)的河網(wǎng)匯流} }Confluence::Confluence(double cs, double ci, double cg, double cr, double im,double qs, double qi, double qg, double qt, double qu,double rs, double ri, double rg, double rim,double qi0, double qg0, double qu0, double f,double u, double m, double csd, double cid, double cgd, double crd, double dt_) {CS = cs;CI = ci;CG = cg;CR = cr;IM = im;QS = qs;QI = qi;QG = qg;QT = qt;QU = qu;RS = rs;RI = ri;RG = rg;RIM = rim;QI0 = qi0;QG0 = qg0;QU0 = qu0;F = f;U = u;M = m;CSD = csd;CID = cid;CGD = cgd;CRD = crd;dt = dt_; }Confluence::~Confluence() {}

河道匯流

Muskingum.h

#pragma once #include "Data.h"//河道匯流：馬斯京根分段連續(xù)演算法（分段馬法） class Muskingum { public:void SetParmameter(const Parameter* parameter); //設(shè)置參數(shù)void SetState(const State* state); //設(shè)置狀態(tài)void UpdateState(State* state); //更新狀態(tài)void calculate();Muskingum(double ke = 0.0, double xe = 0.0, double kl = 0.0, double xl = 0.0, int n = 1, double c0 = 0.0, double c1 = 0.0, double c2 = 0.0,double i1 = 0.0, double i2 = 0.0, double o1 = 0.0, double o2 = 0.0, double* o = nullptr, double dt_ = 24);~Muskingum();protected: private://========模型參數(shù)========//double KE; //馬斯京根法演算參數(shù)，h，敏感，KE = N * ?tdouble XE; //馬斯京根法演算參數(shù)，敏感，0.0~0.5//========模型狀態(tài)========//double KL; //子河段的馬斯京根法演算參數(shù)/hdouble XL; //子河段的馬斯京根法演算參數(shù)int N; //單元河段數(shù)，即分段數(shù)double C0; //馬斯京根流量演算公式I2系數(shù)double C1; //馬斯京根流量演算公式I1系數(shù)double C2; //馬斯京根流量演算公式O1系數(shù)double I1; //時(shí)段初的河段入流量double I2; //時(shí)段末的河段入流量double O1; //時(shí)段初的河段出流量double O2; //時(shí)段末的河段出流量double * O; //各子河段出流量double dt; //模型計(jì)算時(shí)段長,h};

Muskingum.cpp

#include "Muskingum.h"void Muskingum::SetParmameter(const Parameter* parameter) {KE = parameter->m_KE;XE = parameter->m_XE; }void Muskingum::SetState(const State* state) {I1 = state->m_QU0;I2 = state->m_QU;O = state->m_O;dt = state->m_dt; }void Muskingum::UpdateState(State* state) {state->m_O = O;state->m_O2 = O2; }void Muskingum::calculate() {KL = dt; //為了保證馬斯京根法的兩個(gè)線性條件，每個(gè)單元河取 KL = ?tN = int(KE / KL); //單元河段數(shù)XL = 0.5 - N * (1 - 2 * XE) / 2; //計(jì)算單元河段XLdouble denominator = 0.5 * dt + KL - KL * XL; //計(jì)算分母C0 = (0.5 * dt - KL * XL) / denominator;C1 = (0.5 * dt + KL * XL) / denominator;C2 = (-0.5 * dt + KL - KL * XL) / denominator;if (!O){O = new double[N]; //創(chuàng)建存儲(chǔ)單元流域在子河段出口斷面的出流量的動(dòng)態(tài)數(shù)組for (int n = 0; n < N; n++){O[n] = 0.0; //單元流域在子河段出口斷面的出流量為0}}for (int n = 0; n < N; n++){O1 = O[n]; //子河段時(shí)段初出流量O2 = C0 * I2 + C1 * I1 + C2 * O1; //計(jì)算時(shí)段末單元流域在子河段出口斷面的出流量，m3/sO[n] = O2; //更新子河段時(shí)段初出流量I1 = O1; //上一河段時(shí)段初出流為下一河段時(shí)段初入流I2 = O2; //上一河段時(shí)段末出流為下一河段時(shí)段末入流} }Muskingum::Muskingum(double ke, double xe, double kl, double xl, int n, double c0, double c1, double c2,double i1, double i2, double o1, double o2, double* o, double dt_) {KE = ke;XE = xe;KL = kl;XL = xl;N = n;C0 = c0;C1 = c1;C2 = c2;I1 = i1;I2 = i2;O1 = o1;O2 = o2;O = o;dt = dt_; }Muskingum::~Muskingum() {}

參考文獻(xiàn)

---

趙人俊.流域水文模型——新安江模型與陜北模型[M]. 北京：水利電力出版社，1983. ??

包為民.水文預(yù)報(bào)，第5版[M]，北京：中國水利水電出版社，2017 ??

李致家.現(xiàn)代水文模擬與預(yù)報(bào)技術(shù)[M]，南京：河海大學(xué)出版社，2010 ??

翟家瑞.常用水文預(yù)報(bào)算法和計(jì)算程序[M]，鄭州：黃河水利出版社，1995 ??

新安江模型，百度百科 ??

總結(jié)

以上是生活随笔為你收集整理的【水文模型】10 新安江模型C++实现的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。