前言
在水文预报中，利用1小时单位线计算而来的S曲线可以很方便的完成单位线的时段转换。实际运用中，当目标单位线时段长Δtk是原始单位线时段长Δt0的整数倍，即当Δtk=KΔt0 (K=1,2,…N)时，可以不需要计算出Δt0=1h的S曲线即可完成单位线的时段转换，其方法为由原始单位线计算得S(t) (Δt=Δt0)，接着计算得S(t-KΔt)，两个S曲线相减再除以K即可得到目标单位线的q(t)曲线，但当Δtk!=KΔt0 (K=1,2,…N)时，就必须得通过三次样条插值得到Δt=1h的S曲线后再进行接下来的时段转换计算。
插值实现
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 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 | using System; using System.Collections.Generic; using System.Linq; using System.Text; using System.Threading.Tasks; namespace UHTools { public static class LeastSquare { ///<summary> ///用最小二乘法拟合二元多次曲线 ///</summary> ///<param name="arrX">已知点的x坐标集合</param> ///<param name="arrY">已知点的y坐标集合</param> ///<param name="length">已知点的个数</param> ///<param name="dimension">方程的最高次数</param> public static double[] MultiLine(double[] arrX, double[] arrY, int length, int dimension)//二元多次线性方程拟合曲线 { int n = dimension + 1; //dimension次方程需要求 dimension+1个 系数 double[,] Guass = new double[n, n + 1]; //高斯矩阵 例如：y=a0+a1*x+a2*x*x for (int i = 0; i < n; i++) { int j; for (j = 0; j < n; j++) { Guass[i, j] = SumArr(arrX, j + i, length); } Guass[i, j] = SumArr(arrX, i, arrY, 1, length); } return ComputGauss(Guass, n); } public static double SumArr(double[] arr, int n, int length) //求数组的元素的n次方的和 { double s = 0; for (int i = 0; i < length; i++) { if (arr[i] != 0 || n != 0) s = s + Math.Pow(arr[i], n); else s = s + 1; } return s; } public static double SumArr(double[] arr1, int n1, double[] arr2, int n2, int length) { double s = 0; for (int i = 0; i < length; i++) { if ((arr1[i] != 0 || n1 != 0) && (arr2[i] != 0 || n2 != 0)) s = s + Math.Pow(arr1[i], n1) * Math.Pow(arr2[i], n2); else s = s + 1; } return s; } public static double[] ComputGauss(double[,] Guass, int n) { int i, j; int k, m; double temp; double max; double s; double[] x = new double[n]; for (i = 0; i < n; i++) x[i] = 0.0;//初始化 for (j = 0; j < n; j++) { max = 0; k = j; for (i = j; i < n; i++) { if (Math.Abs(Guass[i, j]) > max) { max = Guass[i, j]; k = i; } } if (k != j) { for (m = j; m < n + 1; m++) { temp = Guass[j, m]; Guass[j, m] = Guass[k, m]; Guass[k, m] = temp; } } if (0 == max) { // "此线性方程为奇异线性方程" return x; } for (i = j + 1; i < n; i++) { s = Guass[i, j]; for (m = j; m < n + 1; m++) { Guass[i, m] = Guass[i, m] - Guass[j, m] * s / (Guass[j, j]); } } }//结束for (j=0;j<n;j++) for (i = n - 1; i >= 0; i--) { s = 0; for (j = i + 1; j < n; j++) { s = s + Guass[i, j] * x[j]; } x[i] = (Guass[i, n] - s) / Guass[i, i]; } return x; }//返回值是函数的系数 } public static class SplineMath { /// <summary> /// 三次样条插值 /// </summary> /// <param name="points">排序好的数</param> /// <param name="xs">需要计算的插值点</param> /// <param name="chf">写1</param> /// <returns>返回计算好的数值</returns> public static double[] SplineInsertPoint(PointClass[] points, double[] xs, int chf) { int plength = points.Length; double[] h = new double[plength]; double[] f = new double[plength]; double[] l = new double[plength]; double[] v = new double[plength]; double[] g = new double[plength]; for (int i = 0; i < plength - 1; i++) { h[i] = points[i + 1].x - points[i].x; f[i] = (points[i + 1].y - points[i].y) / h[i]; } for (int i = 1; i < plength - 1; i++) { l[i] = h[i] / (h[i - 1] + h[i]); v[i] = h[i - 1] / (h[i - 1] + h[i]); g[i] = 3 * (l[i] * f[i - 1] + v[i] * f[i]); } double[] b = new double[plength]; double[] tem = new double[plength]; double[] m = new double[plength]; double f0 = (points[0].y - points[1].y) / (points[0].x - points[1].x); double fn = (points[plength - 1].y - points[plength - 2].y) / (points[plength - 1].x - points[plength - 2].x); b[1] = v[1] / 2; for (int i = 2; i < plength - 2; i++) { // Console.Write(" " + i); b[i] = v[i] / (2 - b[i - 1] * l[i]); } tem[1] = g[1] / 2; for (int i = 2; i < plength - 1; i++) { //Console.Write(" " + i); tem[i] = (g[i] - l[i] * tem[i - 1]) / (2 - l[i] * b[i - 1]); } m[plength - 2] = tem[plength - 2]; for (int i = plength - 3; i > 0; i--) { //Console.Write(" " + i); m[i] = tem[i] - b[i] * m[i + 1]; } m[0] = 3 * f[0] / 2.0; m[plength - 1] = fn; int xlength = xs.Length; double[] insertRes = new double[xlength]; for (int i = 0; i < xlength; i++) { int j = 0; for (j = 0; j < plength; j++) { if (xs[i] < points[j].x) break; } j = j - 1; Console.WriteLine(j); if (j == -1 || j == points.Length - 1) { if (j == -1) throw new Exception("插值下边界超出"); if (j == points.Length - 1 && xs[i] == points[j].x) insertRes[i] = points[j].y; else throw new Exception("插值下边界超出"); } else { double p1; p1 = (xs[i] - points[j + 1].x) / (points[j].x - points[j + 1].x); p1 = p1 * p1; double p2; p2 = (xs[i] - points[j].x) / (points[j + 1].x - points[j].x); p2 = p2 * p2; double p3; p3 = p1 * (1 + 2 * (xs[i] - points[j].x) / (points[j + 1].x - points[j].x)) * points[j].y + p2 * (1 + 2 * (xs[i] - points[j + 1].x) / (points[j].x - points[j + 1].x)) * points[j + 1].y; double p4; p4 = p1 * (xs[i] - points[j].x) * m[j] + p2 * (xs[i] - points[j + 1].x) * m[j + 1]; // Console.WriteLine(m[j] + " " + m[j + 1] + " " + j); p4 = p4 + p3; insertRes[i] = p4; //Console.WriteLine("f(" + xs[i] + ")= " + p4); } } //Console.ReadLine(); return insertRes; } } } |
三次样条插值需要使用的PointClass如下：
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 | using System; using System.Collections.Generic; using System.Linq; using System.Text; using System.Threading.Tasks; namespace UHTools { public class PointClass { public double x = 0; public double y = 0; public PointClass() { x = 0; y = 0; } //-------写一个排序函数，使得输入的点按顺序排列，是因为插值算法的要求是，x轴递增有序的--------- public static PointClass[] DeSortX(PointClass[] points) { int length = points.Length; double temx, temy; for (int i = 0; i < length - 1; i++) { for (int j = 0; j < length - i - 1; j++) if (points[j].x > points[j + 1].x) { temx = points[j + 1].x; points[j + 1].x = points[j].x; points[j].x = temx; temy = points[j + 1].y; points[j + 1].y = points[j].y; points[j].y = temy; } } return points; } } } |
使用方法
三次样条插值
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 | double[] t = tSeries.ToArray();//tSeries为时间序列 double[] s = sSeries.ToArray();//sSeries为初始S曲线系列 PointClass[] points = new PointClass[t.Length]; for (int i = 0; i < t.Length; i++) { points[i] = new PointClass(); points[i].x = t[i]; points[i].y = s[i]; } PointClass.DeSortX(points); //排序 List<double> time = new List<double>(); for (int x = 0; x < rowCount; x++) time.Add(x); //求出需要插值的时段 List<double> tInterpolation = time.Except(tSeries).ToList(); // 差集 double[] Y = SplineMath.SplineInsertPoint(points, tInterpolation.ToArray(), 1); |
最小二乘插值
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 | double[] t = tSeries.ToArray(); double[] s = sSeries.ToArray(); int partLength = t.Length / 2; double[] tpart1 = SplitArray(t, 0, partLength - 1); double[] tpart2 = SplitArray(t, partLength, t.Length - 1); double[] spart1 = SplitArray(s, 0, partLength - 1); double[] spart2 = SplitArray(s, partLength, s.Length - 1); double[] efficient1 = LeastSquare.MultiLine(tpart1, spart1, tpart1.Length, 4); double[] efficient2 = LeastSquare.MultiLine(tpart2, spart2, tpart2.Length, 4); //计算St //分部分计算 for (int x = 1; x < tpart1.Length * dtOrigin - 1; x++) { val = efficient1[1] * x + efficient1[2] * Math.Pow(x, 2) + efficient1[3] * Math.Pow(x, 3) + efficient1[4] * Math.Pow(x, 4) + efficient1[0]; Console.WriteLine(val); } for (int x = tpart1.Length * dtOrigin - 1; x < rowCount; x++) { val = efficient2[1] * x + efficient2[2] * Math.Pow(x, 2) + efficient2[3] * Math.Pow(x, 3) + efficient2[4] * Math.Pow(x, 4) + efficient2[0]; Console.WriteLine(val); } |
