前言
马斯京根法是由G.T.麦卡锡于1938年提出，在我国河段流量演算中得到了广泛地运用，经过了几十年的发展，马斯京根法已经相当的成熟，但是实际使用中，马斯京根法依然存在许多限制，在二参数马斯京根法中，存在无区间旁侧入流的假定，若需要用在有旁侧入流的河段中，需要根据具体情况对马斯京根法进行修改，如使用三参数马斯京根法。
基本马斯京根方法
马斯京根法（以下简称马法）中河段内入流、出流及槽蓄量之间的关系以连续方程和槽蓄方程表示：
I-O=dW/dt （1）
W=K[xI+(1-x)O] （2）
式中，I为入流；O为出流 ；W为槽蓄量； K 为槽蓄系数 ( 量纲为时间) ； X 为权重系数 ,，反映入流和出流在决定河段蓄量中的相对重 要性 。
在基本马法中，假定不存在区间旁侧入流，将式（1）以差分形式表示，并将式（2）代入式中即可求得出流量公式：
Q2=C0I1+C1I2+C2Q1
式中，
C0=(0.5 * Δt – K * x) / (0.5 * Δt + K – K * x)，
C1=(0.5 * Δt + K * x) / (0.5 * Δt + K – K * x)，
C2=(-0.5 * Δt + K – K * x)/(0.5 * Δt + K – K * x)。
参数率定
可以看到，对于马法而言，系数K和x的确定十分重要，直接影响了演算结果的精度，对于参数估计，常用的方法是利用图解法试算估计K、x两个参数即根据已 知的实测洪水资料(I , O 值) , 点绘 W 和〔xI + ( 1-x ) O 〕值 , 试错 x 值。若点据接近直线分布 , 则 x 值即可采用 , 直线的坡度即为 K 值 。在计算机中，个人认为采用相关系数判断是否点距时候接近直线分布，实现代码如下：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 | using System; using System.Collections.Generic; using System.IO; using System.Linq; using System.Text; using System.Threading.Tasks; namespace Maskinen { class Program { //程序运行目录 static string ProgrameDir = System.IO.Directory.GetCurrentDirectory(); //演算方程系数 static double[] C = new double[3]; //Δt static double dt = 0; //试算参数 static double K, x; static void Main(string[] args) { dt = double.Parse(Console.ReadLine()); //出入流数据 double[] IData = File.ReadAllLines(ProgrameDir + "\\I.txt").Select(x => double.Parse(x)).ToArray(); double[] OData = File.ReadAllLines(ProgrameDir + "\\O.txt").Select(x => double.Parse(x)).ToArray(); double[] dIQ = new double[IData.Length]; double[] dQ = new double[IData.Length]; double[] S = new double[IData.Length]; for (int i = 1; i < IData.Length - 1; i++) { //ΔQ=I-O dIQ[i] = IData[i] - OData[i]; //ΔQ蓄量变化量 dQ[i] = (IData[i + 1] + IData[i]) / 2 - (OData[i + 1] + OData[i]) / 2; } for (int i = 2; i < IData.Length - 1; i++) { S[i] = S[i - 1] + dQ[i - 1]; } callTrial(OData, dIQ, S); K = Math.Round(K * dt, 0); Console.WriteLine("初步试算结果："); Console.WriteLine("x={0} \t K={1}", x, K); } static void callTrial(double[] OData,double[] dIQ,double[] S) { Console.Write("相关系数阈值R="); double RThreshold = double.Parse(Console.ReadLine()); double[] result = Trial(OData, dIQ, S, 0.10, RThreshold); x = result[0]; K = result[1]; if (x >= 1.0) { Console.WriteLine("计算无结果，请重新指定相关系数阈值"); callTrial(OData, dIQ, S); } } static double[] Trial(double[] O, double[] dIQ, double[] S, double x, double RThreshold) { double[] _Q = new double[S.Length]; for (int i = 1; i < S.Length - 1; i++) { _Q[i] = O[i] + x * dIQ[i]; } //回归前数据处理 List<double> _QList = new List<double>(_Q); List<double> SList = new List<double>(S); _QList.RemoveAt(0); _QList.RemoveAt(_QList.Count - 1); SList.RemoveAt(0); SList.RemoveAt(SList.Count - 1); double[] _Qtmp = _QList.ToArray<double>(); double[] Stmp = SList.ToArray<double>(); //线性回归系数 double[] coff = LeastSquare.MultiLine(Stmp, _Qtmp, Stmp.Length, 1); //回归方程函数值 double[] _QEqua = new double[_Qtmp.Length]; for (int i = 0; i < _Qtmp.Length; i++) { _QEqua[i] = coff[0] + coff[1] * Stmp[i]; } //计算相关系数 double correl = Statistic.getCorrel(_Qtmp, _QEqua); double[] result = new double[2]; result[0] = x; result[1] = coff[1]; if (correl >= RThreshold || x>=1.0) return result; else return Trial(O, dIQ, S, x + 0.05, RThreshold); } } } |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 | using System; public static class Statistic { /// <summary> /// 平均数 /// </summary> /// <param name="Data"></param> /// <returns></returns> public static double getAverage(double[] Data) { double a = 0; foreach(double d in Data) { a += d; } a /= Data.Length; return a; } /// <summary> /// 标准差 /// </summary> /// <param name="Data"></param> /// <returns></returns> public static double getSTD(double[] Data) { double avg = getAverage(Data); double sumstdev = 0; foreach(double d in Data) { sumstdev += (d - avg) * (d - avg); } return Math.Sqrt(sumstdev); } /// <summary> /// 相关系数 /// </summary> /// <param name="Data1"></param> /// <param name="Data2"></param> /// <returns></returns> public static double getCorrel(double[] Data1,double[] Data2) { double avg1 = getAverage(Data1); double std1 = getSTD(Data1); double avg2 = getAverage(Data2); double std2 = getSTD(Data2); double sum = 0; for(int i = 0; i < Data1.Length && i < Data2.Length; i++) { sum += ((Data1[i] - avg1)*(Data2[i] - avg2) / (std2*std1)); } return sum; } /// <summary> /// 剩余标准差 /// </summary> /// <param name="Data1"></param> /// <param name="Data2"></param> /// <returns></returns> public static double getRSD(double[] DataC,double[] DataO) { double sum = 0; for(int i = 0; i < DataC.Length; i++) { sum += Math.Pow(DataC[i] - DataO[i], 2) / (DataC.Length - 2); } sum = Math.Sqrt(sum); return sum; } } |
