From a44a29e8532d0c741675b31ed7c2f3c3c96c2a82 Mon Sep 17 00:00:00 2001 From: biss Date: Tue, 16 Jun 2026 11:43:41 +0800 Subject: [PATCH] =?UTF-8?q?=E6=9C=80=E5=B0=8F=E4=BA=8C=E4=B9=98-=E7=9F=A9?= =?UTF-8?q?=E9=98=B5=E8=BF=90=E7=AE=97?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- Models/LeastSquaresResult.cs | 68 ++++++++ Services/DataCheck.cs | 297 +++++++++++++++++++++++++++++++++++ Services/DataStore.cs | 56 +++++++ Services/LeastSquares.cs | 242 ++++++++++++++++++++++++++++ Services/MatrixOperations.cs | 197 +++++++++++++++++++++++ 5 files changed, 860 insertions(+) create mode 100644 Models/LeastSquaresResult.cs create mode 100644 Services/DataCheck.cs create mode 100644 Services/DataStore.cs create mode 100644 Services/LeastSquares.cs create mode 100644 Services/MatrixOperations.cs diff --git a/Models/LeastSquaresResult.cs b/Models/LeastSquaresResult.cs new file mode 100644 index 0000000..9a880e6 --- /dev/null +++ b/Models/LeastSquaresResult.cs @@ -0,0 +1,68 @@ +using System.Collections.Generic; + +namespace kcsj.Models +{ + public class LeastSquaresResult + { + public int ObservationCount { get; set; } + public int UnknownCount { get; set; } + public int Redundancy { get; set; } + + /// + /// 单位权中误差 + /// + public double Sigma0 { get; set; } + + /// + /// 未知点平差高程 + /// + public Dictionary UnknownElevations { get; } = new(); + + /// + /// 所有点平差后高程，包括已知点和未知点 + /// + public Dictionary AdjustedElevations { get; } = new(); + + /// + /// 未知点高程中误差 + /// + public Dictionary UnknownElevationErrors { get; } = new(); + + /// + /// 每条观测的改正数 v + /// + public List Residuals { get; } = new(); + + /// + /// 每条观测的平差后高差 + /// + public List AdjustedHeightDiffs { get; } = new(); + + /// + /// 每条观测的详细结果 + /// + public List ObservationResults { get; } = new(); + } + + public class ObservationAdjustmentResult + { + public int Index { get; set; } + + public string FromPoint { get; set; } = ""; + public string ToPoint { get; set; } = ""; + + public double ObservedHeightDiff { get; set; } + public double Distance { get; set; } + public double Weight { get; set; } + + /// + /// 改正数 v + /// + public double Residual { get; set; } + + /// + /// 平差后高差 h + v + /// + public double AdjustedHeightDiff { get; set; } + } +} \ No newline at end of file diff --git a/Services/DataCheck.cs b/Services/DataCheck.cs new file mode 100644 index 0000000..2484541 --- /dev/null +++ b/Services/DataCheck.cs @@ -0,0 +1,297 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using kcsj.Models; + +namespace kcsj.Services +{ + public class DataCheckResult + { + public List Errors { get; } = new(); + public List Warnings { get; } = new(); + public List Infos { get; } = new(); + + public bool IsValid => Errors.Count == 0; + + public string ToLogText() + { + var lines = new List(); + + lines.Add("========== 数据检查结果 =========="); + + foreach (string info in Infos) + { + lines.Add("[信息] " + info); + } + + foreach (string warning in Warnings) + { + lines.Add("[警告] " + warning); + } + + foreach (string error in Errors) + { + lines.Add("[错误] " + error); + } + + lines.Add(IsValid ? "数据检查通过。" : "数据检查未通过。"); + lines.Add("================================"); + + return string.Join(Environment.NewLine, lines); + } + } + + public static class DataCheck + { + public static DataCheckResult Check( + List knownPoints, + List observations) + { + var result = new DataCheckResult(); + + if (knownPoints == null || knownPoints.Count == 0) + { + result.Errors.Add("未导入已知点数据。"); + } + else + { + result.Infos.Add($"已知点数量：{knownPoints.Count}"); + } + + if (observations == null || observations.Count == 0) + { + result.Errors.Add("未导入观测数据。"); + } + else + { + result.Infos.Add($"观测数据数量：{observations.Count}"); + } + + if (!result.IsValid) + { + return result; + } + + CheckKnownPoints(knownPoints, result); + CheckObservations(observations, result); + CheckNetwork(knownPoints, observations, result); + CheckRedundancy(knownPoints, observations, result); + + return result; + } + + private static void CheckKnownPoints( + List knownPoints, + DataCheckResult result) + { + var names = new HashSet(StringComparer.OrdinalIgnoreCase); + + for (int i = 0; i < knownPoints.Count; i++) + { + KnownPoint point = knownPoints[i]; + + if (string.IsNullOrWhiteSpace(point.Name)) + { + result.Errors.Add($"第 {i + 1} 个已知点点名为空。"); + continue; + } + + string name = point.Name.Trim(); + + if (!names.Add(name)) + { + result.Errors.Add($"已知点点名重复：{name}"); + } + + if (double.IsNaN(point.Elevation) || double.IsInfinity(point.Elevation)) + { + result.Errors.Add($"已知点 {name} 的高程不是有效数字。"); + } + } + } + + private static void CheckObservations( + List observations, + DataCheckResult result) + { + var edges = new HashSet(StringComparer.OrdinalIgnoreCase); + + for (int i = 0; i < observations.Count; i++) + { + Observation obs = observations[i]; + + if (string.IsNullOrWhiteSpace(obs.FromPoint)) + { + result.Errors.Add($"第 {i + 1} 条观测起点为空。"); + continue; + } + + if (string.IsNullOrWhiteSpace(obs.ToPoint)) + { + result.Errors.Add($"第 {i + 1} 条观测终点为空。"); + continue; + } + + string from = obs.FromPoint.Trim(); + string to = obs.ToPoint.Trim(); + + if (string.Equals(from, to, StringComparison.OrdinalIgnoreCase)) + { + result.Errors.Add($"第 {i + 1} 条观测起点和终点相同：{from}"); + } + + if (double.IsNaN(obs.HeightDiff) || double.IsInfinity(obs.HeightDiff)) + { + result.Errors.Add($"第 {i + 1} 条观测高差不是有效数字。"); + } + + if (double.IsNaN(obs.Distance) || double.IsInfinity(obs.Distance)) + { + result.Errors.Add($"第 {i + 1} 条观测距离不是有效数字。"); + } + else if (obs.Distance <= 0) + { + result.Errors.Add($"第 {i + 1} 条观测距离必须大于 0。"); + } + + string edge1 = from + "-" + to; + string edge2 = to + "-" + from; + + if (edges.Contains(edge1) || edges.Contains(edge2)) + { + result.Warnings.Add($"第 {i + 1} 条观测边可能重复：{from} - {to}"); + } + else + { + edges.Add(edge1); + } + } + } + + private static void CheckNetwork( + List knownPoints, + List observations, + DataCheckResult result) + { + if (!result.IsValid) + { + return; + } + + var knownNames = new HashSet( + knownPoints.Select(p => p.Name.Trim()), + StringComparer.OrdinalIgnoreCase); + + var allNames = new HashSet(StringComparer.OrdinalIgnoreCase); + + foreach (KnownPoint point in knownPoints) + { + allNames.Add(point.Name.Trim()); + } + + foreach (Observation obs in observations) + { + allNames.Add(obs.FromPoint.Trim()); + allNames.Add(obs.ToPoint.Trim()); + } + + var graph = new Dictionary>(StringComparer.OrdinalIgnoreCase); + + foreach (string name in allNames) + { + graph[name] = new List(); + } + + foreach (Observation obs in observations) + { + string from = obs.FromPoint.Trim(); + string to = obs.ToPoint.Trim(); + + graph[from].Add(to); + graph[to].Add(from); + } + + var visited = new HashSet(StringComparer.OrdinalIgnoreCase); + var queue = new Queue(); + + foreach (string knownName in knownNames) + { + if (graph.ContainsKey(knownName)) + { + visited.Add(knownName); + queue.Enqueue(knownName); + } + } + + while (queue.Count > 0) + { + string current = queue.Dequeue(); + + foreach (string next in graph[current]) + { + if (!visited.Contains(next)) + { + visited.Add(next); + queue.Enqueue(next); + } + } + } + + foreach (string name in allNames) + { + if (!visited.Contains(name)) + { + result.Errors.Add($"点 {name} 没有与任何已知点连通。"); + } + } + + int unknownCount = allNames.Count(name => !knownNames.Contains(name)); + + result.Infos.Add($"总点数：{allNames.Count}"); + result.Infos.Add($"未知点数：{unknownCount}"); + } + + private static void CheckRedundancy( + List knownPoints, + List observations, + DataCheckResult result) + { + if (!result.IsValid) + { + return; + } + + var knownNames = new HashSet( + knownPoints.Select(p => p.Name.Trim()), + StringComparer.OrdinalIgnoreCase); + + var allNames = new HashSet(StringComparer.OrdinalIgnoreCase); + + foreach (KnownPoint point in knownPoints) + { + allNames.Add(point.Name.Trim()); + } + + foreach (Observation obs in observations) + { + allNames.Add(obs.FromPoint.Trim()); + allNames.Add(obs.ToPoint.Trim()); + } + + int unknownCount = allNames.Count(name => !knownNames.Contains(name)); + int observationCount = observations.Count; + int redundancy = observationCount - unknownCount; + + result.Infos.Add($"多余观测数：r = {observationCount} - {unknownCount} = {redundancy}"); + + if (observationCount < unknownCount) + { + result.Errors.Add("观测数少于未知点数，无法平差。"); + } + else if (redundancy == 0) + { + result.Warnings.Add("多余观测数为 0，无法进行精度评定。"); + } + } + } +} \ No newline at end of file diff --git a/Services/DataStore.cs b/Services/DataStore.cs new file mode 100644 index 0000000..dfc416e --- /dev/null +++ b/Services/DataStore.cs @@ -0,0 +1,56 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using kcsj.Models; + +namespace kcsj.Services +{ + public static class DataStore + { + public static List KnownPoints { get; private set; } = new(); + public static List Observations { get; private set; } = new(); + + public static string DataSource { get; private set; } = ""; + + public static bool HasData + { + get + { + return KnownPoints.Count > 0 && Observations.Count > 0; + } + } + + public static void SetData( + List knownPoints, + List observations, + string dataSource) + { + if (knownPoints == null || knownPoints.Count == 0) + { + throw new ArgumentException("已知点数据为空，不能保存数据。"); + } + + if (observations == null || observations.Count == 0) + { + throw new ArgumentException("观测数据为空，不能保存数据。"); + } + + KnownPoints = knownPoints.ToList(); + Observations = observations.ToList(); + DataSource = dataSource; + + LogService.AddLog($"数据已更新，来源：{dataSource}"); + LogService.AddLog($"已知点数量：{KnownPoints.Count}"); + LogService.AddLog($"观测数据数量：{Observations.Count}"); + } + + public static void Clear() + { + KnownPoints.Clear(); + Observations.Clear(); + DataSource = ""; + + LogService.AddLog("数据已清空。"); + } + } +} \ No newline at end of file diff --git a/Services/LeastSquares.cs b/Services/LeastSquares.cs new file mode 100644 index 0000000..3e41ee7 --- /dev/null +++ b/Services/LeastSquares.cs @@ -0,0 +1,242 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using kcsj.Models; + +namespace kcsj.Services +{ + public static class LeastSquaresAdjustmentService + { + /// + /// 水准网间接平差 + /// + /// 观测关系： + /// h_AB = H_B - H_A + /// + /// 误差方程： + /// v = Bx - L + /// + /// 法方程： + /// N x = W + /// N = B^T P B + /// W = B^T P L + /// x = N^-1 W + /// + public static LeastSquaresResult Adjust( + List knownPoints, + List observations) + { + if (knownPoints == null || knownPoints.Count == 0) + { + throw new ArgumentException("已知点不能为空。"); + } + + if (observations == null || observations.Count == 0) + { + throw new ArgumentException("观测数据不能为空。"); + } + + Dictionary knownElevations = + knownPoints.ToDictionary( + p => p.Name.Trim(), + p => p.Elevation, + StringComparer.OrdinalIgnoreCase); + + HashSet allPointNames = new(StringComparer.OrdinalIgnoreCase); + + foreach (Observation obs in observations) + { + if (string.IsNullOrWhiteSpace(obs.FromPoint) || + string.IsNullOrWhiteSpace(obs.ToPoint)) + { + throw new ArgumentException("观测数据中存在空点名。"); + } + + allPointNames.Add(obs.FromPoint.Trim()); + allPointNames.Add(obs.ToPoint.Trim()); + } + + List unknownNames = allPointNames + .Where(name => !knownElevations.ContainsKey(name)) + .OrderBy(name => name) + .ToList(); + + if (unknownNames.Count == 0) + { + throw new InvalidOperationException("没有未知点，不需要进行间接平差。"); + } + + int observationCount = observations.Count; + int unknownCount = unknownNames.Count; + + if (observationCount < unknownCount) + { + throw new InvalidOperationException("观测数小于未知数，无法进行最小二乘平差。"); + } + + Dictionary unknownIndex = new(StringComparer.OrdinalIgnoreCase); + + for (int i = 0; i < unknownNames.Count; i++) + { + unknownIndex[unknownNames[i]] = i; + } + + MatrixOperations B = new MatrixOperations(observationCount, unknownCount); + MatrixOperations L = new MatrixOperations(observationCount, 1); + MatrixOperations P = new MatrixOperations(observationCount, observationCount); + + double[] weights = new double[observationCount]; + + for (int i = 0; i < observationCount; i++) + { + Observation obs = observations[i]; + + string from = obs.FromPoint.Trim(); + string to = obs.ToPoint.Trim(); + + // h_AB = H_B - H_A + // FromPoint 是未知点，系数为 -1 + if (unknownIndex.ContainsKey(from)) + { + B[i, unknownIndex[from]] = -1.0; + } + + // ToPoint 是未知点，系数为 +1 + if (unknownIndex.ContainsKey(to)) + { + B[i, unknownIndex[to]] = 1.0; + } + + // 已知点贡献：H_B(已知) - H_A(已知) + double knownContribution = 0.0; + + if (knownElevations.ContainsKey(to)) + { + knownContribution += knownElevations[to]; + } + + if (knownElevations.ContainsKey(from)) + { + knownContribution -= knownElevations[from]; + } + + // L = h观测 - 已知点贡献 + L[i, 0] = obs.HeightDiff - knownContribution; + + // 水准测量常用定权：p = 1 / S + // 距离越长，权越小 + double weight; + + if (obs.Distance <= 0) + { + weight = 1.0; + } + else + { + weight = 1.0 / obs.Distance; + } + + P[i, i] = weight; + weights[i] = weight; + } + + MatrixOperations Bt = B.Transpose(); + + MatrixOperations N = Bt.Multiply(P).Multiply(B); + MatrixOperations W = Bt.Multiply(P).Multiply(L); + + MatrixOperations Qxx = N.Inverse(); + MatrixOperations X = Qxx.Multiply(W); + + // v = Bx - L + MatrixOperations V = B.Multiply(X).Sub(L); + + // V^T P V + MatrixOperations VtPV = V.Transpose().Multiply(P).Multiply(V); + + int redundancy = observationCount - unknownCount; + + double sigma0; + + if (redundancy > 0) + { + sigma0 = Math.Sqrt(VtPV[0, 0] / redundancy); + } + else + { + sigma0 = double.NaN; + } + + LeastSquaresResult result = new LeastSquaresResult + { + ObservationCount = observationCount, + UnknownCount = unknownCount, + Redundancy = redundancy, + Sigma0 = sigma0 + }; + + foreach (KnownPoint point in knownPoints) + { + result.AdjustedElevations[point.Name.Trim()] = point.Elevation; + } + + for (int i = 0; i < unknownCount; i++) + { + string pointName = unknownNames[i]; + double elevation = X[i, 0]; + + result.UnknownElevations[pointName] = elevation; + result.AdjustedElevations[pointName] = elevation; + + if (double.IsNaN(sigma0)) + { + result.UnknownElevationErrors[pointName] = double.NaN; + } + else + { + result.UnknownElevationErrors[pointName] = + sigma0 * Math.Sqrt(Math.Abs(Qxx[i, i])); + } + } + + for (int i = 0; i < observationCount; i++) + { + Observation obs = observations[i]; + + double residual = V[i, 0]; + double adjustedHeightDiff = obs.HeightDiff + residual; + + result.Residuals.Add(residual); + result.AdjustedHeightDiffs.Add(adjustedHeightDiff); + + result.ObservationResults.Add(new ObservationAdjustmentResult + { + Index = i + 1, + FromPoint = obs.FromPoint, + ToPoint = obs.ToPoint, + ObservedHeightDiff = obs.HeightDiff, + Distance = obs.Distance, + Weight = weights[i], + Residual = residual, + AdjustedHeightDiff = adjustedHeightDiff + }); + } + + LogService.AddLog("间接平差计算完成。"); + LogService.AddLog($"观测数：{observationCount}"); + LogService.AddLog($"未知点数：{unknownCount}"); + LogService.AddLog($"多余观测数：{redundancy}"); + + if (double.IsNaN(sigma0)) + { + LogService.AddLog("单位权中误差：无法计算，多余观测数为 0。"); + } + else + { + LogService.AddLog($"单位权中误差：{sigma0:F6}"); + } + + return result; + } + } +} \ No newline at end of file diff --git a/Services/MatrixOperations.cs b/Services/MatrixOperations.cs new file mode 100644 index 0000000..bfa795f --- /dev/null +++ b/Services/MatrixOperations.cs @@ -0,0 +1,197 @@ +using System; +using System.Collections.Generic; +using System.Linq; +using System.Text; +using System.Threading.Tasks; + +namespace kcsj.Services +{ + public class MatrixOperations + { + private const double eps = 1e-2; + private readonly double[,] values; + + private int RowCount { get; } + private int ColumnCount { get; } + + public double this[int row, int column] + { + get => values[row, column]; + set => values[row, column] = value; + } + + // 初始化矩阵大小 + public MatrixOperations(int rowCount, int columnCount) + { + if (rowCount <= 0 || columnCount <= 0) + { + throw new ArgumentException("矩阵行列数必须大于 0。"); + } + + RowCount = rowCount; + ColumnCount = columnCount; + values = new double[rowCount, columnCount]; + } + + // 从二维数组初始化矩阵 + public MatrixOperations(double[,] source) + { + RowCount = source.GetLength(0); + ColumnCount = source.GetLength(1); + values = new double[RowCount, ColumnCount]; + + for (int i = 0; i < RowCount; i++) + { + for (int j = 0; j < ColumnCount; j++) + { + values[i, j] = source[i, j]; + } + } + } + + public MatrixOperations Add(MatrixOperations other) + { + if (RowCount != other.RowCount || ColumnCount != other.ColumnCount) + { + throw new ArgumentException("矩阵维度不匹配，无法相加。"); + } + MatrixOperations result = new MatrixOperations(RowCount, ColumnCount); + for (int i = 0; i < RowCount; i++) + { + for (int j = 0; j < ColumnCount; j++) + { + result[i, j] = this[i, j] + other[i, j]; + } + } + return result; + } + + public MatrixOperations Sub (MatrixOperations other) + { + if (RowCount != other.RowCount || ColumnCount != other.ColumnCount) + { + throw new ArgumentException("矩阵维度不匹配，无法相减。"); + } + MatrixOperations result = new MatrixOperations(RowCount, ColumnCount); + for (int i = 0; i < RowCount; i++) + { + for (int j = 0; j < ColumnCount; j++) + { + result[i, j] = this[i, j] - other[i, j]; + } + } + return result; + + } + + public MatrixOperations Multiply(MatrixOperations other) + { + if (ColumnCount != other.RowCount) + { + throw new ArgumentException("矩阵维度不匹配，无法相乘。"); + } + MatrixOperations result = new MatrixOperations(RowCount, other.ColumnCount); + for (int i = 0; i < RowCount; i++) + { + for (int j = 0; j < other.ColumnCount; j++) + { + double sum = 0; + for (int k = 0; k < ColumnCount; k++) + { + sum += this[i, k] * other[k, j]; + } + result[i, j] = sum; + } + } + return result; + } + + public MatrixOperations Transpose() + { + MatrixOperations result = new MatrixOperations(ColumnCount, RowCount); + for (int i = 0; i < RowCount; i++) + { + for (int j = 0; j < ColumnCount; j++) + { + result[j, i] = this[i, j]; + } + } + return result; + } + + public MatrixOperations Inverse() + { + if (RowCount != ColumnCount) + { + throw new ArgumentException("只有方阵才有逆矩阵。"); + } + int n = RowCount; + MatrixOperations augmented = new MatrixOperations(n, 2 * n); + // 构造增广矩阵 [A | I] + for (int i = 0; i < n; i++) + { + for (int j = 0; j < n; j++) + { + augmented[i, j] = this[i, j]; + } + augmented[i, n + i] = 1; // 添加单位矩阵部分 + } + // 使用高斯消元法将左侧变为单位矩阵 + for (int i = 0; i < n; i++) + { + // 寻找主元素 + int pivotRow = i; + for (int row = i + 1; row < n; row++) + { + if (Math.Abs(augmented[row, i]) > Math.Abs(augmented[pivotRow, i])) + { + pivotRow = row; + } + } + if (Math.Abs(augmented[pivotRow, i]) < eps) + { + throw new InvalidOperationException("矩阵不可逆。"); + } + // 交换行 + if (pivotRow != i) + { + for (int col = 0; col < 2 * n; col++) + { + double temp = augmented[i, col]; + augmented[i, col] = augmented[pivotRow, col]; + augmented[pivotRow, col] = temp; + } + } + // 将主元素归一化 + double pivotValue = augmented[i, i]; + for (int col = 0; col < 2 * n; col++) + { + augmented[i, col] /= pivotValue; + } + // 消去其他行的当前列 + for (int row = 0; row < n; row++) + { + if (row != i) + { + double factor = augmented[row, i]; + for (int col = 0; col < 2 * n; col++) + { + augmented[row, col] -= factor * augmented[i, col]; + } + } + } + } + // 提取右侧的逆矩阵 + MatrixOperations inverse = new MatrixOperations(n, n); + for (int i = 0; i < n; i++) + { + for (int j = 0; j < n; j++) + { + inverse[i, j] = augmented[i, n + j]; + } + } + return inverse; + } + + } +}