1 public class XinAnJiangModel { 2 // FORCING 驱动数据 3 private double[] m_pP; // 降水数据 4 private double[] m_pEm; // 水面蒸发数据 5 // 6 private int m_nSteps; // 模型要运行的步长 7 8 // OUTPUT 输出数据 9 private double[] m_pR; // 流域内每一步长的产流量(径流深度) 10 private double[] m_pRg; // 每一步长的地表径流深(毫米) 11 private double[] m_pRs; // 每一步长的基流径流深(毫米) 12 private double[] m_pE; // 每一步长的蒸发(毫米) 13 private double[] m_pQrs; // 流域出口地表径流量 14 private double[] m_pQrg; // 流域出口地下径流量 15 private double[] m_pQ; // 流域出口的总流量 16 // 17 private double m_U; // for 24h. U=A(km^2)/3.6/delta_t 18 // SOIL 土壤数据 19 private double[] m_pW; // 流域内土壤湿度 20 private double[] m_pWu; // 流域内上层土壤湿度 21 private double[] m_pWl; // 流域内下层土壤湿度 22 private double[] m_pWd; // 流域内深层土壤湿度 23 24 private double m_Wum; // 流域内上层土壤蓄水容量，植被良好的流域，约为20mm,差的流域,2~10mm 25 private double m_Wlm; // 流域内下层土壤蓄水容量，可取60~90mm 26 private double m_Wdm; // 流域内深层土壤蓄水容量，WDM=WM-WUM-WLM 27 28 // EVAPORATION 蒸发 29 private double[] m_pEu; // 上层土壤蒸发量（毫米） 30 private double[] m_pEl; // 下层土壤蒸发量（毫米） 31 private double[] m_pEd; // 深层土壤蒸发量（毫米） 32 33 // PARAMETER 模型参数 34 private double m_K; // 流域蒸散发能力与实测蒸散发值的比 35 private double m_IMP; // 不透水面积占全流域面积之比 36 private double m_B; // 蓄水容量曲线的方次，小流域（几平方公里）B为0.1左右， 37 // 中等面积（300平方公里以内）0.2~0.3，较大面积0.3~0.4 38 private double m_WM; // 流域平均蓄水容量（毫米）(WM=WUM+WLM+WDM) 39 40 private double m_C; // 流域内深层土壤蒸发系数，江南湿润地区：0.15-0.2，华北半湿润地区：0.09-0.12 41 private double m_FC; // 稳定入渗率，毫米/小时 42 private double m_KKG; // 地下径流消退系数 43 //double m_UH; // 单元流域上地面径流的单位线 44 private double m_Kstor; // 脉冲汇流计算的参数,Liang 45 private double m_WMM; // 流域内最大蓄水容量 46 private double m_Area; // 流域面积 47 private int m_DeltaT; // 每一步长的小时数 48 //double m_U; // 49 50 public void XinanjiangModel () 51 { 52 this.m_pE = null; 53 this.m_pP = null; 54 55 this.m_pE = null; 56 this.m_pEd = null; 57 this.m_pEl = null; 58 this.m_pEu = null; 59 60 this.m_pQrg = null; 61 this.m_pQrs = null; 62 63 this.m_pR = null; 64 this.m_pRg = null; 65 this.m_pRs = null; 66 }; 67 68 // 模型析构函数 69 protected void finalize() 70 { 71 this.m_pP = null; 72 this.m_pEm = null; 73 74 this.m_pE = null; 75 this.m_pEd = null; 76 this.m_pEl = null; 77 this.m_pEu = null; 78 79 this.m_pQrg = null; 80 this.m_pQrs = null; 81 this.m_pQ = null; 82 83 this.m_pR = null; 84 this.m_pRg = null; 85 this.m_pRs = null; 86 87 this.m_pW = null; 88 this.m_pWd = null; 89 this.m_pWl = null; 90 this.m_pWu = null; 91 }; 92 93 94 // 通过文件初始化模型设置，包括步长，步数，驱动数据 95 public void XinanjiangModel (File FileName) 96 { 97 98 }; 99 100 // 初始化模型 101 public void InitModel (int nSteps, double Area, int DeltaT, File ForcingFile) 102 { 103 File fp; 104 int i; 105 this.m_nSteps = nSteps; 106 // 驱动数据 107 this.m_pP = new double[this.m_nSteps]; 108 this.m_pEm = new double[this.m_nSteps]; 109 // 模型输出，蒸散发项 110 this.m_pE = new double[this.m_nSteps]; 111 this.m_pEd = new double[this.m_nSteps]; 112 this.m_pEl = new double[this.m_nSteps]; 113 this.m_pEu = new double[this.m_nSteps]; 114 // 模型输出，出流项，经过汇流的产流 115 this.m_pQrg = new double[this.m_nSteps]; 116 this.m_pQrs = new double[this.m_nSteps]; 117 this.m_pQ = new double[this.m_nSteps]; 118 // 模型输出，产流项 119 this.m_pR = new double[this.m_nSteps]; 120 this.m_pRg = new double[this.m_nSteps]; 121 this.m_pRs = new double[this.m_nSteps]; 122 // 模型状态量，土壤湿度 123 this.m_pW = new double[this.m_nSteps]; 124 this.m_pWd = new double[this.m_nSteps]; 125 this.m_pWl = new double[this.m_nSteps]; 126 this.m_pWu = new double[this.m_nSteps]; 127 128 this.m_Area = Area; 129 this.m_DeltaT = DeltaT; 130 this.m_U = 1.0; //this.m_Area/(3.6*this.m_DeltaT); 131 // Forcing文件的格式：第一列：降水（单位毫米）空格第二列水面蒸发（毫米） 132 fp = ForcingFile; 133 134 try 135 { 136 if (!fp.exists()||fp.isDirectory()) 137 throw new FileNotFoundException(); 138 BufferedReader in = new BufferedReader(new FileReader(fp)); 139 String line; //一行数据 140 // int row=0; 141 //逐行读取，并将每个数组放入到数组中 142 // while((line = in.readLine()) != null){ 143 // String[] temp = line.split("t"); 144 // 145 // for(int j=0;j<temp.length;j++){ 146 // arr2[row][j] = Double.parseDouble(temp[j]); 147 // this.m_pP[i] = 148 // this.m_pEm[i]= 149 // } 150 // row++; 151 // } 152 for (i = 0; i < this.m_nSteps; i++) //读取驱动数据 153 { 154 // fscanf (fp, "%lf%lf", &(this.m_pP[i]), &(this.m_pEm[i])); 155 line = in.readLine(); 156 if (line != null){ 157 String[] temp = line.split("\t"); 158 this.m_pP[i] = Double.parseDouble(temp[0]); 159 this.m_pEm[i]= Double.parseDouble(temp[1]); 160 } 161 } 162 in.close(); 163 } 164 catch (Exception e) 165 { 166 System.out.println("找不到指定的驱动数据！"+ e.getMessage()); 167 return; 168 } 169 170 }; 171 172 // 初始化模型 173 public void InitModel (int nSteps, double Area, int DeltaT, double[][] forcingData) 174 { 175 File fp; 176 int i; 177 this.m_nSteps = nSteps; 178 // 驱动数据 179 this.m_pP = new double[this.m_nSteps]; 180 this.m_pEm = new double[this.m_nSteps]; 181 // 模型输出，蒸散发项 182 this.m_pE = new double[this.m_nSteps]; 183 this.m_pEd = new double[this.m_nSteps]; 184 this.m_pEl = new double[this.m_nSteps]; 185 this.m_pEu = new double[this.m_nSteps]; 186 // 模型输出，出流项，经过汇流的产流 187 this.m_pQrg = new double[this.m_nSteps]; 188 this.m_pQrs = new double[this.m_nSteps]; 189 this.m_pQ = new double[this.m_nSteps]; 190 // 模型输出，产流项 191 this.m_pR = new double[this.m_nSteps]; 192 this.m_pRg = new double[this.m_nSteps]; 193 this.m_pRs = new double[this.m_nSteps]; 194 // 模型状态量，土壤湿度 195 this.m_pW = new double[this.m_nSteps]; 196 this.m_pWd = new double[this.m_nSteps]; 197 this.m_pWl = new double[this.m_nSteps]; 198 this.m_pWu = new double[this.m_nSteps]; 199 200 this.m_Area = Area; 201 this.m_DeltaT = DeltaT; 202 this.m_U = 1.0; //this.m_Area/(3.6*this.m_DeltaT); 203 // Forcing数据格式：第一维：降水（单位毫米）；第二维水面蒸发（毫米） 204 205 for (i = 0; i < this.m_nSteps; i++) //读取驱动数据 206 { 207 // fscanf (fp, "%lf%lf", &(this.m_pP[i]), &(this.m_pEm[i])); 208 this.m_pP[i] = forcingData[i][0]; 209 this.m_pEm[i]= forcingData[i][1]; 210 } 211 212 }; 213 214 // 设置模型参数 215 public void SetParameters (double[] Params) 216 { 217 this.m_K = Params[0]; // (1) 流域蒸散发能力与实测水面蒸发之比 218 this.m_IMP = Params[1]; // (2) 流域不透水面积占全流域面积之比 219 this.m_B = Params[2]; // (3) 蓄水容量曲线的方次 220 this.m_Wum = Params[3]; // (4) 上层蓄水容量 221 this.m_Wlm = Params[4]; // (5) 下层蓄水容量 222 this.m_Wdm = Params[5]; // (6) 深层蓄水容量 223 this.m_C = Params[6]; // (7) 深层蒸散发系数 224 this.m_FC = Params[7]; // (8) 稳定入渗率（毫米／小时） 225 this.m_KKG = Params[8]; // (9) 地下径流消退系数 226 this.m_Kstor = Params[9]; // (10)汇流计算参数 227 228 this.m_WM = this.m_Wum + this.m_Wlm + this.m_Wdm; 229 this.m_WMM = this.m_WM * (1.0 + this.m_B) / (1.0 - this.m_IMP); 230 }; 231 // 运行新安江模型 232 public void RunModel () 233 { 234 int i; 235 236 // 模型的状态变量 237 double PE; // 大于零时为净雨量，小于零时为蒸发不足量，单位（毫米） 238 239 double R; // 产流深度，包括地表径流深度和地下径流深度两部分（毫米） 240 double RG; // 地下径流深度，单位（毫米） 241 double RS; // 地表径流深度，单位（毫米） 242 243 double A; // 当流域内的土壤湿度为W是,土壤含水量折算成的径流深度,单位（毫米） 244 245 double E = 0.0; // 蒸散发 246 double EU = 0.0; // 上层土壤蒸散发量（毫米） 247 double EL = 0.0; // 下层土壤蒸散发量（毫米） 248 double ED = 0.0; // 深层土壤蒸散发量（毫米） 249 250 // 假设流域经历了长时间降水，各层土壤含水量均为该层土壤的蓄水能力 251 double W = this.m_WM; // 流域内土壤湿度 252 double WU = this.m_Wum; // 流域内上层土壤湿度 253 double WL = this.m_Wlm; // 流域内下层土壤适度 254 double WD = this.m_Wdm; // 流域内深层土壤湿度 255 256 for (i = 0; i < this.m_nSteps; i++) 257 { 258 PE = m_pP[i] - m_K * m_pEm[i]; 259 // 如果降水量小于足蒸发需求 260 if (PE < 0) 261 { 262 R = 0.0; // 产流总量为零 263 RG = 0.0; // 地下径流量为零 264 RS = 0.0; // 地表径流量为零 265 266 // 如果上层土壤含水量大于蒸发不足量 267 if ((WU + PE) > 0.0) 268 { 269 // 上层土壤为流域蒸散发提供水量 270 EU = m_K * m_pEm[i]; 271 // 没有降水量用于增加土壤湿度 272 EL = 0.0; /* 降水用来增加土壤湿度的部分 */ 273 // 274 ED = 0.0; 275 // 更新上层土壤含水量 276 WU = WU + PE; 277 } 278 // 上层土壤含水量小于蒸发不足量 279 else 280 { 281 EU = WU + m_pP[i]; // 上层土壤蒸发,降水全部用于蒸发 282 WU = 0.0; // 上层含水量为0，全部水分被蒸发 283 // 如果下层含水量大于下层土壤的蒸散发潜力 284 if (WL > (m_C * m_Wlm)) 285 { 286 EL = (m_K * m_pEm[i] - EU) * (WL / m_Wlm); 287 WL = WL - EL; 288 ED = 0; 289 } 290 // 如果下层土壤含水量小于下层土壤的蒸散发潜力 291 else 292 { 293 // 如果下层土壤的含水量蒸发之后还有剩余 294 if (WL > m_C * (m_K * m_pEm[i] - EU)) 295 { 296 EL = m_C * (m_K * m_pEm[i] - EU); 297 WL = WL - EL; ///////////////////////////////// 298 ED = 0.0; 299 } 300 // 如果下层土壤含水量全部蒸发之后尚不能满足蒸发需求 301 else 302 { 303 EL = WL; /* 下层土壤含水量全部用于蒸散发 */ 304 WL = 0.0; /* 下层土剩余壤含水量为0 */ 305 ED = m_C * (m_K * m_pEm[i] - EU) - EL; /* 深层土壤含水量参与蒸发 */ 306 WD = WD - ED; /* 深层土壤含水量更新 */ 307 } 308 } 309 } 310 } 311 // 如果降水量大于或者等于蒸散发需求,即降水满足蒸发后还有剩余 312 else 313 { 314 /*************** 以下代码负责径流划分计算 **************/ 315 // 初始化变量 316 R = 0.0; // 产流总量为零 317 RG = 0.0; // 地下径流产流深度为零 318 RS = 0.0; // 地表径流产流深度为零 319 320 // 计算流域当天土壤含水量折算成的径流深度Ａ 321 // m_WM:流域平均蓄水容量(一个参数), 322 // m_W:流域内土壤湿度(一个状态变量) 323 // m_B:蓄水容量曲线的方次(一个参数) 324 A = m_WMM * (1 - Math.pow((1.0 - W / m_WM), 1.0 / (1 + m_B))); 325 // 土壤湿度折算净雨量加上降水后蒸发剩余雨量小于流域内最大含水容量 326 if ((A + PE) < this.m_WMM) 327 { 328 // 流域内的产流深度计算 329 R = PE /* 降水蒸发后的剩余量(PE=P-E:状态变量) */ 330 + W /* 流域内土壤湿度 (W:状态变量) */ 331 + m_WM * Math.pow ((1 - (PE + A) / m_WMM), (1 + m_B)) - m_WM; /* 减去流域平均蓄水容量（m_WM:参数） */ 332 } 333 // 土壤湿度折算净雨量加上降水后蒸发剩余雨量大于流域内最大含水容量 334 else 335 { 336 // 流域内的产流深度计算 337 R = PE /* 降水蒸发后的剩余量(PE=P-E:状态变量) */ 338 + W /* 流域内土壤湿度 (W:状态变量) */ 339 - m_WM; /* 减去流域平均蓄水容量（m_WM:参数） */ 340 } 341 // 如果降水经过蒸散发后的剩余量大于等于土壤稳定入渗率// 342 if (PE > m_FC) 343 { 344 // 计算地下径流的产流深度 345 RG = (R - this.m_IMP * PE) * (m_FC/ PE); 346 // 计算地表径流的产流深度 347 RS = R - RG; 348 } 349 // 如果降水蒸发后的剩余量小于土壤的稳定入渗率(m_FC:参数) 350 //除了不透水面积上的地表产流外，全部下渗，形成地下径流 351 else 352 { 353 // 计算地下径流的产流深度 354 RG = R - /* 总产流深度 */ 355 m_IMP * PE; /* 不透水面积上的产流深度，m_IMP:参数 */ 356 // 计算地表径流的产流深度 357 RS = R - RG; 358 } 359 /*************** 径流划分计算结束 **************/ 360 361 // 计算上层土壤蒸散发量 362 EU = m_K * /* 流域蒸散发能力与实测蒸散发值的比 */ 363 m_pEm[i]; /* 当前时段的水面蒸发 */ 364 ED = 0.0; 365 EL = 0.0; /* 降水用来增加土壤湿度的部分 */ 366 367 /*************** 以下代码负责土壤含水量的更新计算 **************/ 368 // 如果上层土壤含水量与降水蒸散发剩余量之和减去产流量之后 369 // 大于上层土壤的蓄水能力 370 if ((WU + PE - R) >= m_Wum) 371 { 372 // 上层含水量+下层含水量+降水剩余量-产流量-上层土壤蓄水需求 373 // 后的水量大于下层土壤蓄水需求，多余水量分配到深层土壤 374 if ((WU + WL + PE - R - m_Wum) > m_Wlm) 375 { 376 WU = m_Wum; /* 上层土壤含水量=上层土壤蓄水容量 */ 377 WL = m_Wlm; /* 下层土壤含水量=下层土壤蓄水容量 */ 378 WD = W + PE - R - WU - WL; /* 绝对降水剩余量补充到深层土壤中 */ 379 } 380 // 上层含水量+下层含水量+降水剩余量-产流量-上层土壤蓄水需求 381 // 后的水量小于下层土壤蓄水需求，剩余水量补充到下层土壤中 382 else 383 { 384 WL = WU + WL + PE - R - m_Wum; /* 下层土壤含水量 */ 385 WU = m_Wum; /* 上层土壤蓄水容量得到满足 */ 386 } 387 } 388 // 如果上层土壤含水量与降水蒸散发剩余量之和减去产流量之后 389 // 小于上层土壤的蓄水能力 390 else 391 { 392 WU = WU + PE - R; 393 // WU 可能小于零，应该加一些控制代码.......... 394 } 395 /*************** 土壤含水量的更新计算结束 **************/ 396 } 397 398 E = EU + EL + ED; 399 W = WU + WL + WD; 400 /* 以下部分是状态量：总蒸发量、上、下和深层土壤的蒸发的保存 */ 401 /* 1 */ this.m_pE[i] = E; 402 // 当前步长的蒸发 （模型重要输出） 403 /* 2 */ this.m_pEu[i] = EU; 404 // 当前步长上层土壤蒸发 405 /* 3 */ this.m_pEl[i] = EL; 406 // 当前步长下层土壤蒸发 407 /* 4 */ this.m_pEd[i] = ED; 408 // 当前步长深层土壤蒸发 409 /* 5 */ this.m_pW[i] = W; 410 // 当前步长流域平均土壤含水量 411 /* 6 */ this.m_pWu[i] = WU; 412 // 当前步长流域上层土壤含水量 413 /* 7 */ this.m_pWl[i] = WL; 414 // 当前步长流域下层土壤含水量 415 /* 8 */ this.m_pWd[i] = WD; 416 // 当前步长流域深层土壤含水量 417 /* 9 */ this.m_pRg[i] = RG; 418 // 当前步长流域基流径流深度 419 /* 10*/ this.m_pRs[i] = RS; 420 // 当前步长流域地表径流径流深度 421 /* 11*/ this.m_pR[i] = R; 422 // 当前步长的总产流径流深度 423 } 424 this.Routing (); 425 }; 426 427 // 保存模拟结果到文件 428 public void SaveResults (String FileName) 429 { 430 DataWriter outWriter; 431 432 int i; 433 try 434 { 435 // if(!FileName.exists()) 436 // FileName.createNewFile(); 437 438 // fprintf (fp," E(mm) EU(mm) EL(mm) ED(mm) W(mm) WU(mm) WL(mm) WD(mm) R(mm) RG(mm) RS(mm) Q(m3/d) QS(m3/d) QG(m3/d) "); 439 outWriter = new DataWriter(FileName); 440 System.out.println("模型运算结果输出至文件开始"); 441 outWriter.writeLine(" E(mm) EU(mm) EL(mm) ED(mm) W(mm) WU(mm) WL(mm) WD(mm) R(mm) RG(mm) RS(mm) Q(m3/d) QS(m3/d) QG(m3/d)"); 442 443 for (i = 0; i < this.m_nSteps; i++) 444 { 445 // fprintf (fp," %11.5lf %11.5lf %11.5lf %11.5lf %11.5lf %11.5lf %11.5lf %11.5lf %11.5lf %11.5lf %11.5lf %11.5lf %11.5lf %11.5lf ", 446 // this.m_pE[i], this.m_pEu[i], this.m_pEl[i], this.m_pEd[i], 447 // this.m_pW[i], this.m_pWu[i], this.m_pWl[i], this.m_pWd[i], 448 // this.m_pR[i], this.m_pRs[i], this.m_pRg[i], this.m_pQ[i], 449 // this.m_pQrs[i], this.m_pQrg[i]); 450 outWriter.writeLine(String.format("%1$16.3f %2$16.3f %3$16.3f %4$16.3f %5$16.3f %6$16.3f %7$16.3f %8$16.3f %9$16.3f %10$16.3f %11$16.3f %12$16.3f %13$16.3f %14$16.3f", 451 this.m_pE[i], this.m_pEu[i], this.m_pEl[i], this.m_pEd[i], 452 this.m_pW[i], this.m_pWu[i], this.m_pWl[i], this.m_pWd[i], 453 this.m_pR[i], this.m_pRs[i], this.m_pRg[i], this.m_pQ[i], 454 this.m_pQrs[i], this.m_pQrg[i])); 455 456 } 457 System.out.println("模型运算结果输出至文件结束"); 458 // out.close(); 459 } 460 catch(Exception e) 461 { 462 System.out.println("无法创建模型参数输出文件！"); 463 return; 464 } 465 }; 466 public void Runoff (double[] runoff) 467 { 468 /*从1990年1月1日到1996年12月31日为模型的标定期，共有2557天，其总从 469 1990年1月1日到1990年12月31日为模型运行的预热期，不参与标定 */ 470 int i; 471 472 for (i = 0; i < 2556; i++) 473 { 474 runoff[i] = this.m_pQ[i]; 475 // if(runoff[i]<0) 476 //{ 477 //// printf("runoff=%lf ",this.m_pQ[i]); 478 //// getch(); 479 // runoff[i]=0.0; 480 //} 481 } 482 }; 483 484 // 进行汇流计算，将径流深度转换为流域出口的流量 485 private void Routing () 486 { 487 double[] UH = new double[100]; // 单位线,假定最长的汇流时间为100天 488 int N ; // 汇流天数 489 N = 0; 490 double K; // 汇流参数 491 double sum; 492 int i, j; 493 494 K = this.m_Kstor; 495 // 单位线推导 496 for (i = 0; i < 100; i++) 497 { 498 UH[i] = (1.0 / K) * Math.exp((-1.0 * i) / K); 499 } 500 UH[0]=(UH[1]+UH[2])*0.5; 501 sum = 0.0; 502 for (i = 0; i < 100; i++) 503 { 504 sum += UH[i]; 505 if (sum > 1.0) 506 { 507 UH[i] = 1.0 - (sum - UH[i]); 508 N = i; 509 break; 510 } 511 } 512 // 单位线汇流计算 513 for (i = 0; i < this.m_nSteps; i++) 514 { 515 this.m_pQrs[i] = 0.0; 516 for (j = 0; j <= N; j++) 517 { 518 if ((i - j) < 0) 519 { 520 continue; 521 } 522 this.m_pQrs[i] += this.m_pRs[i - j] * UH[j] * this.m_U; 523 } 524 } 525 //地下水汇流计算 526 this.m_pQrg[0] = 0.0; 527 for (i = 1; i < this.m_nSteps; i++) 528 { 529 this.m_pQrg[i] = this.m_pQrg[i - 1] * this.m_KKG + 530 this.m_pRg[i] * (1.0 - this.m_KKG) * this.m_U; 531 } 532 for (i = 0; i < this.m_nSteps; i++) 533 { 534 this.m_pQ[i] = this.m_pQrs[i] + this.m_pQrg[i]; 535 } 536 }; 537 }