Files
SzmediTools/GeologyToolkit/KrigingMethod.cs
2025-09-16 16:06:41 +08:00

653 lines
22 KiB
C#

using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
using MathNet.Numerics;
using MathNet.Numerics.LinearAlgebra.Complex;
using MathNet.Numerics.LinearRegression;
namespace GeologyToolkit
{
internal class KrigingMethod
{
/// <summary>
/// 求已知点的距离矩阵
/// </summary>
/// <param name="points"></param>
/// <returns></returns>
public DenseMatrix DistanceDenseMatrix(List<XYPoint> points)
{
int n = points.Count;
var matrix = new DenseMatrix(n, n);
Complex[,] dis = new Complex[n, n];
for (int i = 0; i < n; i++) //行
{
for (int j = 0; j < n; j++) //列
{
Complex distance = points[i].DistanceTo(points[j]);
dis.SetValue(distance, i, j);
}
}
matrix = DenseMatrix.OfArray(dis);
return matrix;
}
/// <summary>
/// 求属性值的半方差矩阵
/// </summary>
/// <param name="points"></param>
public DenseMatrix ElevationVarD(List<XYPoint> points)
{
int n = points.Count;
Complex[,] dis = new Complex[points.Count, points.Count];
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
//属性半方差
Complex var_D = Math.Pow(points[i].Elevation - points[j].Elevation, 2) / 2;
dis.SetValue(var_D, i, j);
}
}
return DenseMatrix.OfArray(dis);
}
/// <summary>
/// 分组后,每组距离的平均值-半方差值平均值(用于拟合半方差函数)
/// </summary>
/// <param name="pairs"></param>
/// <param name="groupCount"></param>
/// <returns></returns>
public List<Pair> GetAveragePair(List<Pair> pairs, int groupCount)
{
var orderedData = pairs.OrderBy(x => x.Distance.Real);
var max = orderedData.Last().Distance;
var min = orderedData.First().Distance;
double step = (max.Real - min.Real + 0.0001) / groupCount;
//分组
List<List<Pair>> ds = new List<List<Pair>>();
for (int i = 0; i < groupCount; i++)
{
List<Pair> d = new List<Pair>();
foreach (var pair in pairs)
{
if (pair.Distance.Real >= step * i + min.Real && pair.Distance.Real <= (i + 1) * step + min.Real)
{
d.Add(pair);
}
}
ds.Add(d);
//IEnumerable<Pair> d =
// orderedData.TakeWhile(n =>
// n.Distance.Real > step * i + min.Real && n.Distance.Real <= (i + 1) * step + min.Real);
//ds.Add(d);
}
List<Pair> averagePairs = new List<Pair>();
foreach (var d in ds)
{
Complex totalDis = 0.0;
Complex totalD = 0.0;
foreach (var complex in d)
{
totalDis += complex.Distance;
totalD += complex.VarD;
}
var pair = new Pair
{
Distance = totalDis / d.Count(),
VarD = totalD / d.Count()
};
averagePairs.Add(pair);
}
return averagePairs;
}
public List<Pair> GroupyDistances(DenseMatrix distances, DenseMatrix elevationD)
{
List<Pair> groups = new List<Pair>();
for (int i = 0; i < distances.RowCount; i++)
{
for (int j = 0; j < elevationD.RowCount; j++)
{
var r = new Pair
{
Distance = distances[i, j],
VarD = elevationD[i, j]
};
groups.Add(r);
}
}
return groups;
}
public double MainMethod(List<XYPoint> points, XYPoint pointToInvestigate, int groupCount = 10)
{
var distancesMatrix = DistanceDenseMatrix(points);
var elevationsVarD = ElevationVarD(points); //高程属性值的半方差
List<Pair> pairs = new List<Pair>();
for (int i = 0; i < distancesMatrix.RowCount; i++)
{
for (int j = 0; j < elevationsVarD.RowCount; j++)
{
for (int k = 0; k < distancesMatrix.ColumnCount; k++)
{
Pair pair = new Pair
{
Distance = distancesMatrix.At(i, k),
VarD = elevationsVarD.At(j, k)
};
pairs.Add(pair);
}
}
}
var fitData = GetAveragePair(pairs, groupCount);
//通过平均距离-半方差,得到拟合函数(距离-半方差关系函数)
double[] x = new double[fitData.Count];
double[] y = new double[fitData.Count];
for (int i = 0; i < x.Length; i++)
{
x[i] = fitData[i].Distance.Real;
y[i] = fitData[i].VarD.Real;
}
//指数模型
var formula_P = Fit.ExponentialFunc(x, y);
//球形模型
var formula_S = Fit.PolynomialFunc(x, y, 3, DirectRegressionMethod.NormalEquations);
//求未知点的距离列矩阵
var unknownDistancesMatrix = UnknownPointDistances(pointToInvestigate, points);
//代入拟合函数,求未知点的半方差矩阵
Complex[] rows = new Complex[unknownDistancesMatrix.Values.Length];
for (var i = 0; i < unknownDistancesMatrix.Values.Length; i++)
{
var d = unknownDistancesMatrix.Values[i];
var distance = d.Real;
//半方差值
var ps = formula_S(distance);
rows[i] = ps;
}
var unknownVarDMatrix = DenseMatrix.OfColumnArrays(rows);
//得到权重(矩阵的逆乘以未知点的距离半方差矩阵)
var inverse = elevationsVarD.Inverse();
var wi = inverse * unknownVarDMatrix;
double n = default;
for (int i = 0; i < wi.RowCount; i++)
{
double na = wi[i, 0].Real;
n += na;
}
//权重乘以对应属性,加权得到
double prop = default;
for (int i = 0; i < points.Count; i++)
{
XYPoint point = points[i];
prop += point.Elevation * wi[i, 0].Real;
}
return prop;
//var data = matrix.Values;
//var orderedData = data.OrderBy(x => x);
//var max = data.Max();
//double step = (max.Real + 10) / groupCount;
//List<IEnumerable<Complex>> ds = new List<IEnumerable<Complex>>();
//for (int i = 0; i < 10; i++)
//{
// IEnumerable<Complex> d = orderedData.TakeWhile(n => n.Real > step * i && n.Real < (i + 1) * step);
// ds.Add(d);
//}
//foreach (var d in ds)
//{
// Complex total = 0.0;
// foreach (var complex in d)
// {
// total += complex;
// }
// var average = total / d.Count();
//}
//matrix.Inverse();
}
/// <summary>
/// 求未知点距离列矩阵(代入拟合半方差函数,可求得半方差)
/// </summary>
/// <param name="pointToInvestigate"></param>
/// <param name="points"></param>
public DenseMatrix UnknownPointDistances(XYPoint pointToInvestigate, List<XYPoint> points)
{
//DenseMatrix matrix = new DenseMatrix(points.Count, 1);
Complex[] rows = new Complex[points.Count];
for (int i = 0; i < points.Count; i++)
{
Complex complex = pointToInvestigate.DistanceTo(points[i]);
rows[i] = complex;
}
return DenseMatrix.OfColumnArrays(rows);
}
private double[,] Deinv(double[,] dMatrix, int Level)
{
double dMatrixValue = MatrixValue(dMatrix, Level);
if (dMatrixValue == 0) return null;
double[,] dReverseMatrix = new double[Level, 2 * Level];
double x, c;
// Init Reverse matrix
for (int i = 0; i < Level; i++)
{
for (int j = 0; j < 2 * Level; j++)
{
if (j < Level)
dReverseMatrix[i, j] = dMatrix[i, j];
else
dReverseMatrix[i, j] = 0;
}
dReverseMatrix[i, Level + i] = 1;
}
for (int i = 0, j = 0; i < Level && j < Level; i++, j++)
{
if (dReverseMatrix[i, j] == 0)
{
int m = i;
for (; dMatrix[m, j] == 0; m++) ;
if (m == Level)
return null;
// Add i-row with m-row
for (int n = j; n < 2 * Level; n++)
dReverseMatrix[i, n] += dReverseMatrix[m, n];
}
// Format the i-row with "1" start
x = dReverseMatrix[i, j];
if (x != 1)
{
for (int n = j; n < 2 * Level; n++)
if (dReverseMatrix[i, n] != 0)
dReverseMatrix[i, n] /= x;
}
// Set 0 to the current column in the rows after current row
for (int s = Level - 1; s > i; s--)
{
x = dReverseMatrix[s, j];
for (int t = j; t < 2 * Level; t++)
dReverseMatrix[s, t] -= dReverseMatrix[i, t] * x;
}
}
// Format the first matrix into unit-matrix
for (int i = Level - 2; i >= 0; i--)
{
for (int j = i + 1; j < Level; j++)
if (dReverseMatrix[i, j] != 0)
{
c = dReverseMatrix[i, j];
for (int n = j; n < 2 * Level; n++)
dReverseMatrix[i, n] -= c * dReverseMatrix[j, n];
}
}
double[,] dReturn = new double[Level, Level];
for (int i = 0; i < Level; i++)
for (int j = 0; j < Level; j++)
dReturn[i, j] = dReverseMatrix[i, j + Level];
return dReturn;
}
//private List<double> pointToInvestigate(XYPoint pointToInvestigate, List<XYPoint> points)
//{
// List<double> dis = new List<double>();
// int n = points.Count;
// for (int i = 0; i < n; i++)
// {
// var distance = pointToInvestigate.DistanceTo(points[i]);
// dis.Add(distance);
// }
// return dis;
//}
/// <summary>
/// 求已知点的距离矩阵数组(无依赖)
/// </summary>
/// <param name="points"></param>
/// <returns></returns>
private double[,] DistanceMatrix(List<XYPoint> points)
{
double[,] dis = new double[points.Count + 1, points.Count];
int n = points.Count;
for (int i = 0; i <= n; i++)
{
for (int j = 0; j < n; j++)
{
double distance = points[i].DistanceTo(points[j]);
dis.SetValue(distance, i, j);
}
}
return dis;
}
private void FitExponential(double[] x, double[] y)
{
var f = Fit.ExponentialFunc(x, y);
double d = 15.1;
var height = f.Invoke(d);
}
/// <summary>
/// 高斯模型
/// </summary>
/// <param name="h">采样点距离</param>
/// <param name="a">变程,与起始采样点距离,</param>
private void FitMethod(double h, double a)
{
if (h == 0)
{
}
else if (h > 0)
{
}
}
private double MatrixValue(double[,] MatrixList, int Level)
{
double[,] dMatrix = new double[Level, Level];
for (int i = 0; i < Level; i++)
for (int j = 0; j < Level; j++)
dMatrix[i, j] = MatrixList[i, j];
double c, x;
int k = 1;
for (int i = 0, j = 0; i < Level && j < Level; i++, j++)
{
if (dMatrix[i, j] == 0)
{
int m = i;
for (; dMatrix[m, j] == 0; m++) ;
if (m == Level)
return 0;
// Row change between i-row and m-row
for (int n = j; n < Level; n++)
{
c = dMatrix[i, n];
dMatrix[i, n] = dMatrix[m, n];
dMatrix[m, n] = c;
}
// Change value pre-value
k *= -1;
}
// Set 0 to the current column in the rows after current row
for (int s = Level - 1; s > i; s--)
{
x = dMatrix[s, j];
for (int t = j; t < Level; t++)
dMatrix[s, t] -= dMatrix[i, t] * (x / dMatrix[i, j]);
}
}
double sn = 1;
for (int i = 0; i < Level; i++)
{
if (dMatrix[i, i] != 0)
sn *= dMatrix[i, i];
else
return 0;
}
return k * sn;
}
/// <summary>
/// 普通克里金
/// </summary>
/// <param name="dpx">数据点坐标x</param>
/// <param name="dpy">数据点坐标y</param>
/// <param name="dpz">数据点坐标z</param>
/// <param name="ndp">数据点个数</param>
/// <param name="ipx">待插入点的坐标x</param>
/// <param name="ipy">待插入点的坐标y</param>
/// <param name="ipz">待插入点的坐标z</param>
/// <param name="nip">待插入点的个数</param>
private void OrdinaryKriging(double[] dpx, double[] dpy, double[] dpz, int ndp, double[] ipx, double[] ipy, double[] ipz, int nip)
{
int LM = 10; //距离分组时的组数
double DR1, DZ;
double[] a = new double[3]; //存储球状模型三个系数
double[] p, A, B, C; //矩阵
int[] LD = new int[LM];
double[] AR = new double[LM];
double[] AZ = new double[LM];
double[] DR = new double[LM + 1];
p = new double[(ndp - 1) * ndp / 2]; //距离矩阵
A = new double[Convert.ToInt16(Math.Pow(ndp + 1, 2))]; //系数矩阵
for (int i = 0; i < Math.Pow(ndp + 1, 2); i++)
{
A[i] = 0;
}
int numOfDis = 0;
for (int i = 0; i < ndp - 1; i++)
{
for (int j = 0; j < ndp; j++)
{
p[numOfDis] = Math.Sqrt(Math.Pow(dpx[i] - dpx[j], 2) + Math.Pow(dpx[i] - dpx[j], 2));
A[i * (ndp + 1) + j] = p[numOfDis]; //记录距离矩阵上半部分
A[j * (ndp + 1) + i] = p[numOfDis]; //记录距离矩阵下半部分
numOfDis++;
}
}
//堆排序,由小到大
p = p.OrderBy(x => x).ToArray();
DR[0] = 0.0;
for (int i = 1; i < LM + 1; i++)
{
var k = i * numOfDis / LM - 1;
DR[i] = p[k];
}
for (int i = 1; i < LM; i++)
{
LD[i] = 0;
AR[i] = 0;
AZ[i] = 0;
}
for (int i = 0; i < ndp - 1; i++)
{
for (int j = 0; j < ndp; j++)
{
DR1 = Math.Sqrt(Math.Pow(dpx[i] - dpx[j], 2) + Math.Pow(dpx[i] - dpx[j], 2)); //距离
DZ = Math.Pow(dpz[i] - dpx[j], 2); //平方增量
for (int k = 0; k <= LM; k++)
{
if (DR1 > DR[k - 1] && DR1 <= DR[k])
{
LD[k - 1] += 1; //点个数
AR[k - 1] += DR1; //每组距离和
AZ[k - 1] += DZ; //每组平方增量和
}
}
}
}
//二元线性回归
//L11*X1+L12*X2=L10
//L21*X1+L22*X2=L20
double X1 = 0.0;
double X2 = 0.0;
double Y = 0.0;
double L11 = 0.0;
double L12 = 0.0;
double L21 = 0.0;
double L22 = 0.0;
double L10 = 0.0;
double L20 = 0.0;
for (int i = 0; i < LM; i++)
{
AR[i] = AR[i] / LD[i]; //距离平均值
AZ[i] = AZ[i] / LD[i] / 2; //平方增量和
X1 += AR[i] / LM;
X2 += Math.Pow(AR[i], 3);
Y += AZ[i] / LM;
}
for (int i = 0; i < LM; i++)
{
L10 += (AR[i] - X1) * (AZ[i] - Y);
L11 += Math.Pow(AR[i] - X1, 2);
L12 += (AR[i] - X1) * (Math.Pow(AR[i], 3) - X2);
L12 += (Math.Pow(AR[i], 3) - X2) * (AZ[i] - Y);
L21 += (Math.Pow(AR[i], 3) - X2) * (AR[i] - X1);
L22 += (Math.Pow(AR[i], 3) - X2) * (Math.Pow(AR[i], 3) - X2);
//球状模型系数
double B1 = (L10 * L22 - L20 * L12) / (L11 * L22 - L12 * L12);
double B2 = (L20 * L11 - L10 * L21) / (L11 * L22 - L12 * L12);
double B0 = Y - B1 * X1 - B2 * X2;
if (B0 > 0)
{
a[0] = B0;
a[1] = Math.Sqrt(-B1 / 3 / B2);
a[2] = 2 * B1 / 3 * Math.Sqrt(-B1 / 3 / B2);
}
}
for (int i = 0; i < ndp; i++)
{
for (int j = 0; j < ndp; j++)
{
if (i == j)
{
a[i * (ndp + 1) + j] = 0.0;
}
else
{
if (A[i * (ndp + 1) + j] < a[1]) //小于变程
{
A[i * (ndp + 1) + j] = a[0] + a[2] * (1.5 * A[i * (ndp + 1) + j] / a[1]) - 0.5 * Math.Pow(A[i * (ndp + 1) + j] / a[1], 3);
}
else
{
A[i * (ndp + 1) + j] = a[0] + a[2]; //大于等于变程
}
}
}
}
for (int i = 0; i < ndp; i++) //对每个数据点加入值为1的项
{
A[(i + 1) * (ndp + 1) - 1] = 1.0;
A[ndp * (ndp + 1) + i] = 1.0;
}
//系数矩阵求逆
//Deinv(A, ndp + 1);
B = new double[ndp + 1]; //存放矩阵方程组等号右边的列向量
C = new double[ndp + 1]; //存放结果
//待插入点循环
for (int i = 0; i < nip; i++)
{
for (int j = 0; j < ndp; j++)
{
DR1 = Math.Sqrt(Math.Pow(ipx[i] - dpx[j], 2) + Math.Pow(ipy[i] - dpy[j], 2));
if (DR1 < 0.00000000001)
{
B[j] = 0.0;
}
else
{
if (DR1 < a[1])
{
B[j] = a[0] + a[2] * (1.5 * DR1 / a[1] - 0.5 * Math.Pow(DR1 / a[1], 3));
}
else
{
B[j] = a[0] + a[2];
}
}
}
B[ndp] = 1.0; //加入值为1的项
for (int k = 0; k < ndp + 1; k++) //插值方程的解
{
C[k] = 0.0;
for (int j = 0; j < ndp + 1; j++)
{
C[k] += A[k * (ndp + 1) + j] * B[j];
}
}
ipz[i] = 0.0;
for (int j = 0; j < ndp; j++)
{
ipz[i] += C[j] * dpz[j]; //普通克里金估值
}
}
}
}
public class XYPoint
{
public XYPoint(double x, double y)
{
X = x;
Y = y;
}
public XYPoint()
{
X = 0.0;
Y = 0.0;
}
public double Elevation { get; set; }
public double X { get; set; }
public double Y { get; set; }
/// <summary>
/// 水平距离
/// </summary>
/// <param name="point"></param>
/// <returns></returns>
public double DistanceTo(XYPoint point)
{
double d = Math.Pow(X - point.X, 2) + Math.Pow(Y - point.Y, 2);
return Math.Sqrt(d);
}
}
public class Pair
{
/// <summary>
/// 距离
/// </summary>
public Complex Distance { get; set; }
/// <summary>
/// 半方差
/// </summary>
public Complex VarD { get; set; }
}
}