模型簡介

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

源碼下載

該項目已在Github開源，您可以點擊此處下載。

個人編寫，難免有理解不到位和疏漏的地方，敬請諒解，僅供參考！
感謝導師、大李老師、小李老師、姚老師、小蕾在模型原理、程序設計、率定驗證等方面的指導！

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

設計思路

嚴格按照“流域分單元，蒸散發(fā)分層次，產流分水源，匯流分階段”的建模流程，以面向對象的方式編寫三水源新安江模型C++程序。

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

程序驗證

程序模塊

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

流域分塊

Watershed.h

#pragma once #include <string> #include <vector>class IO //從文本中導入降雨和蒸發(fā)數(shù)據(jù)，輸出流量過程到文本中 { public:void ReadFromFile(std::string StrPath = ""); //從文本中導入流域降雨或蒸發(fā)數(shù)據(jù)void WriteToFile(std::string StrPath = ""); //輸出流域出口斷面流量過程到文本中IO();~IO();public:std::vector< std::vector<double> > m_P; //各雨量站逐時段降雨量，mmstd::vector< std::vector<double> > m_EM; //各蒸發(fā)站逐時段水面蒸發(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); //計算各單元流域降雨量和水面蒸發(fā)量double GetP(int nt, int nw); //得到nt時刻，第nw單元流域的降雨量，mmdouble GetEM(int nt, int nw); //得到nt時刻，第nw單元流域的水面蒸發(fā)量，mmdouble GetF(int nw); //得到第nw單元流域的面積，km2int GetnW(); //得到單元流域個數(shù)Watershed();~Watershed();protected: private:std::string m_name; //流域名稱double m_area; //流域控制面積，km2int m_numRainfallStation; //雨量站個數(shù)int m_numEvaporationStation; //蒸發(fā)站個數(shù)int m_numSubWatershed; //單元流域個數(shù)double* m_areaSubWatershed; //單元流域面積，km2std::vector< std::vector<double> > m_rateRainfallStation; //各單元流域對應各雨量站比例，按泰森多邊形分配std::vector< std::vector<double> > m_rateEvaporationStation; //各單元流域對應各蒸發(fā)站比例，按泰森多邊形分配std::string* m_nameRainfallStation; //雨量站點名稱std::string* m_nameEvaporationStation; //蒸發(fā)站點名稱std::vector< std::vector<double> > m_P; //各單元流域逐時段降雨量，mmstd::vector< std::vector<double> > m_EM; //各單元流域逐時段水面蒸發(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)建子流域行*雨量站列的二維動態(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ā)站列的二維動態(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];}}//讀取雨量站點名稱m_nameRainfallStation = new std::string[m_numRainfallStation];for (int j = 0; j < m_numRainfallStation; j++){fin >> m_nameRainfallStation[j];}//讀取蒸發(fā)站點名稱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ù)int ncols = m_numSubWatershed; //單元流域個數(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]; //按比例計算單元流域降雨量}}}//計算各單元流域逐時段水面蒸發(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]; //按比例計算單元流域水面蒸發(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 = "默認流域";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ù)fin >> nrows >> ncols;//創(chuàng)建動態(tài)二維數(shù)據(jù)，用于存儲各雨量站降雨量，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ù)表示蒸發(fā)站數(shù)fin >> nrows >> ncols;//創(chuàng)建動態(tài)二維數(shù)據(jù)，用于存儲各蒸發(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ù)void SetState(const State* state); //設置狀態(tài)void UpdateState(State* state); //更新狀態(tài)void calculate(); //計算三層蒸散發(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); //構造函數(shù)~Evapotranspiration(); //析構函數(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ā)計算EP = KC * EM; //計算流域蒸發(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;}}}//計算總的蒸散發(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() {}

產流

Source.h

#pragma once #include "Data.h"//水源劃分 class Source { public:void SetParmameter(const Parameter* parameter); //設置參數(shù)void SetState(const State* state); //設置狀態(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; //表層自由水蓄水庫對地下水的日出流系數(shù)，敏感double KI; //表層自由水蓄水庫對壤中流的日出流系數(shù)，敏感//========模型狀態(tài)========//double R; //總徑流量，mmdouble RS; //地面徑流，mmdouble RI; //壤中流，mmdouble RG; //地下徑流，mmdouble PE; //凈雨量，mmdouble FR; //本時段產流面積比例double S0; //本時段初的自由水蓄量，mmdouble S; //本時段的自由水蓄量，mmdouble M; //一天劃分的計算時段數(shù)double KID; //表層自由水蓄水庫對壤中流的計算時段出流系數(shù)，敏感double KGD; //表層自由水蓄水庫對地下水的計算時段出流系數(shù)，敏感double SMM; //全流域單點最大的自由水蓄水容量，mmdouble SMMF; //產流面積上最大一點的自由水蓄水容量，mmdouble SMF; //產流面積上的平均自由水蓄水容量深，mmdouble AU; //相應平均蓄水深的最大蓄水深，S0值對應的縱坐標，mmdouble RSD; //計算步長地面徑流，mmdouble RID; //計算步長壤中流，mmdouble RGD; //計算步長地下徑流，mmdouble FR0; //上一時段產流面積比例int N; //N 為計算時段分段數(shù)，每一段為計算步長double Q; //Q 是每個計算步長內的凈雨量，mmdouble KIDD; //表層自由水蓄水庫對壤中流的計算步長出流系數(shù)，敏感double KGDD; //表層自由水蓄水庫對地下水的計算步長出流系數(shù)，敏感double dt; //模型計算時段長,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; //一天劃分的計算時段數(shù)KID = (1 - pow(1 - (KI + KG), 1.0 / M)) / (1 + KG / KI); //表層自由水蓄水庫對壤中流的計算時段出流系數(shù)，敏感KGD = KID * KG / KI; //表層自由水蓄水庫對地下水的計算時段出流系數(shù)，敏感//三分水源[4]if (PE <= 1e-5){//凈雨量小于等于0時,這里認為凈雨量小于1e-5時即為小于等于0RS = 0;RI = KID * S0 * FR; /*當凈雨量小于等于0時，消耗自由水蓄水庫中的水產流面積比例仍為上一時段的比例不變*/RG = KGD * S0 * FR;S = S0 * (1 - KID - KGD); //更新下一時段初的自由水蓄量}else{//凈雨量大于0時SMM = SM * (1 + EX); //全流域單點最大的自由水蓄水容量，mmFR0 = FR; //上一時段產流面積比例FR = R / PE; //計算本時段產流面積比例if (FR > 1) //如果FR由于小數(shù)誤差而計算出大于1的情況，則強制置為1{FR = 1;}S = S0 * FR0 / FR;N = int(PE / 5.0) + 1; //N 為計算時段分段數(shù)，每一段為計算步長Q = PE / N; //Q 是每個計算步長內的凈雨量，mm//R 是總徑流量，PE是單元流域降雨深度。//把R按5mm分，實際操作就是把PE按5mm分，是為了減小產流面積變化的差分誤差。//自由水蓄水庫只發(fā)生在產流面積上，其底寬為產流面積FR，顯然FR是隨時間變化的。//產流量R進入水庫即在產流面積上產生PE的徑流深。//FR = f / F = R / PE KIDD = (1 - pow(1 - (KID + KGD), 1.0 / N)) / (1 + KGD / KID); //表層自由水蓄水庫對壤中流的計算步長出流系數(shù)，敏感KGDD = KIDD * KGD / KID; //表層自由水蓄水庫對地下水的計算步長出流系數(shù)，敏感//把該時段的RS、RI、RG置0，用于后續(xù)累加計算步長內的RSD、RID、RGDRS = 0.0;RI = 0.0;RG = 0.0;//計算產流面積上最大一點的自由水蓄水容量 SMMFif (EX == 0.0){SMMF = SMM; //EX等于0時，流域自由水蓄水容量分布均勻}else{//假定SMMF與產流面積FR及全流域上最大點的自由水蓄水容量SMM仍為拋物線分布//則SMMF應該用下式計算SMMF = (1 - pow(1 - FR, 1.0 / EX)) * SMM;}SMF = SMMF / (1.0 + EX);//將每個計算時段的入流量R分成N段，計算各個計算步長內的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){//計算步長內的凈雨量進入自由水蓄水庫后，使得自由水蓄水深超過產流面積上最大單點自由水蓄水深RSD = (Q + S - SMF) * FR;RID = SMF * KIDD * FR;RGD = SMF * KGDD * FR;S = SMF * (1 - KIDD - KGDD);}else{//自由水蓄水深未超過產流面積上最大單點自由水蓄水深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;}}//累加計算時段內的地面徑流、壤中流和地下徑流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ù)void SetState(const State* state); //設置狀態(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)總入流//（進入單元面積的地面徑流、壤中流和地下徑流之和），m3/sdouble QU; //單元流域出口流量，m3/sdouble RS; //地面徑流量，mmdouble RI; //壤中流徑流量，mmdouble RG; //地下徑流量，mmdouble RIM; //不透水面積上的產流量，mmdouble QI0; //QI(t-1)，前一時刻壤中流，m3/sdouble QG0; //QG(t-1)，前一時刻地下徑流，m3/sdouble QU0; //QU(t-1)，前一時刻單元流域出口流量，m3/sdouble F; //單元流域面積，km2double U; //單位轉換系數(shù)double M; //一天劃分的計算時段數(shù)double CSD; //計算時段內地面徑流蓄水庫的消退系數(shù)double CID; //計算時段內壤中流蓄水庫的消退系數(shù)double CGD; //計算時段內地下水蓄水庫的消退系數(shù)double CRD; //計算時段內河網(wǎng)蓄水消退系數(shù)double dt; //模型計算時段長,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() {//把透水面積上的產流量均攤到單元流域上RS = RS * (1 - IM);RI = RI * (1 - IM);RG = RG * (1 - IM);//出流系數(shù)換算M = 24.0 / dt; //一天劃分的計算時段數(shù)CSD = pow(CS, 1.0 / M); //計算時段內地面徑流蓄水庫的消退系數(shù)CID = pow(CI, 1.0 / M); //計算時段內壤中流蓄水庫的消退系數(shù)CGD = pow(CG, 1.0 / M); //計算時段內地下水蓄水庫的消退系數(shù)CRD = pow(CR, 1.0 / M); //計算時段內河網(wǎng)蓄水消退系數(shù)//坡地匯流U = F / 3.6 / 24; //單位轉換系數(shù)QS = CSD * QS + (1 - CSD) * (RS + RIM) * U; //地面徑流流入地面徑流水庫，//經(jīng)過消退(CSD)，成為地面徑流對河網(wǎng)的總入流QSQI = CID * QI + (1 - CID) * RI * U; //壤中流流入壤中流水庫，//經(jīng)過消退(CID)，成為壤中流對河網(wǎng)的總入流QIQG = CGD * QG + (1 - CGD) * RG * U; //地下徑流進入地下水蓄水庫，經(jīng)過地下水//蓄水庫的消退(CGD)，成為地下水對河網(wǎng)的總入流QGQT = QS + QI + QG;//河網(wǎng)匯流，采用線性水庫法，且僅當單元流域面積大于200km2時才計算河網(wǎng)匯流if (F < 200){QU = QT; //單元流域面積不大且河道較短，對水流運動的調蓄作用通常較小//將這種調蓄作用合并在地面徑流和地下徑流中一起考慮所帶來的誤差通常可以忽略}else{QU = CRD * QU + (1 - CRD) * QT; //線性水庫法//只有在單元流域面積較大或流域坡面匯流極其復雜時//才考慮單元面積內的河網(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ù)void SetState(const State* state); //設置狀態(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; //時段初的河段入流量double I2; //時段末的河段入流量double O1; //時段初的河段出流量double O2; //時段末的河段出流量double * O; //各子河段出流量double dt; //模型計算時段長,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; //為了保證馬斯京根法的兩個線性條件，每個單元河取 KL = ?tN = int(KE / KL); //單元河段數(shù)XL = 0.5 - N * (1 - 2 * XE) / 2; //計算單元河段XLdouble denominator = 0.5 * dt + KL - KL * XL; //計算分母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)建存儲單元流域在子河段出口斷面的出流量的動態(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]; //子河段時段初出流量O2 = C0 * I2 + C1 * I1 + C2 * O1; //計算時段末單元流域在子河段出口斷面的出流量，m3/sO[n] = O2; //更新子河段時段初出流量I1 = O1; //上一河段時段初出流為下一河段時段初入流I2 = O2; //上一河段時段末出流為下一河段時段末入流} }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() {}

參考文獻

---

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

包為民.水文預報，第5版[M]，北京：中國水利水電出版社，2017 ??

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

翟家瑞.常用水文預報算法和計算程序[M]，鄭州：黃河水利出版社，1995 ??

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

總結

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

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