#include <iostream>
#include <fstream>
#include <vector>
#include <string>
#include <cmath>
#include <map>
#include <set>
#include <queue>
#include <algorithm>
#include <
stdexcept>
#include <iomanip>
#include <sstream>
#include <unordered_set>
#include <stack>
#include <functional>
#ifdef _WIN32
#include <windows.h> // 用于GetCurrentDirectoryA
#else
#include <uni
std.h> // 用于getcwd
#endif
constexpr double EPSILON = 1e-10;
class Point {
public
:
std::string id;
double elevation;
bool isFixed;
Point(const
std::string& _id, double _elevation = 0.0, bool _isFixed = false)
: id(_id), elevation(_elevation), isFixed(_isFixed) {
}
};
class Observation {
public
:
std::string from;
std::string to;
double dh;
double weight;
double correction;
double adjustedValue;
Observation(const
std::string& _from, const
std::string& _to, double _dh, double _weight = 1.0)
: from(_from), to(_to), dh(_dh), weight(_weight), correction(0.0), adjustedValue(0.0) {
}
};
class Condition {
public
:
std::vector<size_t> obsIndices;
std::vector<double> coefficients;
double w;
Condition()
: w(0.0) {}
void addObservation(size_t obsIdx, double coeff) {
obsIndices.push_back(obsIdx);
coefficients.push_back(coeff);
}
void calculateW(const
std::vector<Observation>& observations) {
w = 0.0;
for (size_t i = 0; i < obsIndices.size(); ++i) {
if (obsIndices[i] < observations.size()) {
w += coefficients[i] * observations[obsIndices[i]].dh;
}
}
w = -w;
}
};
class Matrix {
private
:
size_t rows;
size_t cols;
std::vector<double> data;
public
:
Matrix(size_t r, size_t c)
: rows(r), cols(c), data(r* c, 0.0) {}
// 访问
元素
double& operator()(size_t i, size_t j) {
if (i >= rows || j >= cols) {
throw
std::out_of_range("矩阵索引越界");
}
return data[i * cols + j];
}
const double& operator()(size_t i, size_t j) const {
if (i >= rows || j >= cols) {
throw
std::out_of_range("矩阵索引越界");
}
return data[i * cols + j];
}
size_t getRows() const { return rows; }
size_t getCols() const { return cols; }
// 矩阵乘法
Matrix operator*(const Matrix& other) const {
if (cols != other.rows) {
throw
std::invalid_argument("矩阵维度不匹配");
}
Matrix result(rows, other.cols);
for (size_t i = 0; i < rows; ++i) {
for (size_t j = 0; j < other.cols; ++j) {
for (size_t k = 0; k < cols; ++k) {
result(i, j) += (*this)(i, k) * other(k, j);//C(i，j）=A（i*k)*B(k*j）
}
}
}
return result;
}
// 矩阵转置
Matrix transpose() const {
Matrix result(cols, rows);
for (size_t i = 0; i < rows; ++i) {
for (size_t j = 0; j < cols; ++j) {
result(j, i) = (*this)(i, j);//算A^(T）
}
}
return result;
}
// 更稳定的高斯消元法
static
std::vector<double> solve(const Matrix& A, const
std::vector<double>& b) {
if (A.rows != A.cols || A.rows != b.size()) {
throw
std::invalid_argument("矩阵维度不匹配");
}
size_t n = A.rows;
Matrix augmented(n, n + 1);
for (size_t i = 0; i < n; ++i) {
for (size_t j = 0; j < n; ++j) {
augmented(i, j) = A(i, j);
}
augmented(i, n) = b[i];
}
// 部分选主元
for (size_t i = 0; i < n; ++i) {
// 寻找主元
size_t maxRow = i;
for (size_t k = i + 1; k < n; ++k) {
if (
std::abs(augmented(k, i)) >
std::abs(augmented(maxRow, i))) {
maxRow = k;
}
}
// 交换行
if (maxRow != i) {
for (size_t j = i; j <= n; ++j) {
std::swap(augmented(i, j), augmented(maxRow, j));
}
}
// 奇异矩阵检查
if (
std::abs(augmented(i, i)) < EPSILON) {
throw
std::runtime_error("矩阵奇异，无法
求解");
}
// 消元
for (size_t k = i + 1; k < n; ++k) {
double factor = augmented(k, i) / augmented(i, i);
for (size_t j = i; j <= n; ++j) {
augmented(k, j) -= factor * augmented(i, j);
}
}
}
// 回代
std::vector<double> x(n);
for (int i = static_cast<int>(n) - 1; i >= 0; --i) {
x[i] = augmented(i, n);
for (size_t j = i + 1; j < n; ++j) {
x[i] -= augmented(i, j) * x[j];
}
x[i] /= augmented(i, i);
}
return x;
}
};
class LevelingNetwork {
private
:
std::vector<Point> points;
std::vector<Observation> observations;
std::vector<Condition> conditions;
std::map<
std::string, size_t> pointIndex;
double unitWeightVariance;
std::vector<double> pointPrecisions;
std::vector<double> initialElevations; // 保存初始高程
// 构建图结构
struct GraphNode {
std::string id;
std::vector<
std::pair<
std::string, size_t>> neighbors; // 邻居点和观测索引
};
// 深度优先搜索识别路径
void dfs(const
std::string& current, const
std::string& end,
std::vector<size_t>& path,
std::unordered_set<
std::string>& visited,
std::vector<
std::vector<size_t>>& allPaths,
const
std::map<
std::string, GraphNode>& graph) {
visited.insert(current);
if (current == end) {
allPaths.push_back(path);
visited.erase(current);
return;
}
auto it = graph.find(current);
if (it != graph.end()) {
for (const auto& neighbor
: it->second.neighbors) {
if (visited.find(neighbor.first) == visited.end()) {
path.push_back(neighbor.second);
dfs(neighbor.first, end, path, visited, allPaths, graph);
path.pop_back();
}
}
}
visited.erase(current);
}
// 识别所有可能的路径（用于附合路线）
std::vector<
std::vector<size_t>> findAllPaths(const
std::string& start,
const
std::string& end,
const
std::map<
std::string, GraphNode>& graph) {
std::vector<
std::vector<size_t>> allPaths;
std::vector<size_t> path;
std::unordered_set<
std::string> visited;
if (graph.find(start) != graph.end() && graph.find(end) != graph.end()) {
dfs(start, end, path, visited, allPaths, graph);
}
return allPaths;
}
public
:
void addPoint(const Point& p) {
points.push_back(p);
pointIndex[p.id] = points.size() - 1;
}
void addObservation(const Observation& obs) {
observations.push_back(obs);
}
// 改进的条件方程识别算法
void identifyConditions() {
conditions.clear();
std::map<
std::string, GraphNode> graph;
// 构建图
for (size_t i = 0; i < observations.size(); ++i) {
const auto& obs = observations[i];
graph[obs.from].id = obs.from;
graph[obs.from].neighbors.push_back({ obs.to, i });
graph[obs.to].id = obs.to;
graph[obs.to].neighbors.push_back({ obs.from, i });
}
// 识别闭合环
std::vector<
std::string> knownPoints;
for (const auto& p
: points) {
if (p.isFixed) {
knownPoints.push_back(p.id);
}
}
// 1. 处理附合路线（已知点之间的所有可能路径）
if (knownPoints.size() >= 2) {
for (size_t i = 0; i < knownPoints.size(); ++i) {
for (size_t j = i + 1; j < knownPoints.size(); ++j) {
auto paths = findAllPaths(knownPoints[i], knownPoints[j], graph);
for (const auto& path
: paths) {
if (!path.empty()) {
Condition cond;
for (size_t obsIdx
: path) {
cond.addObservation(obsIdx, 1.0);
}
cond.calculateW(observations);
// 计算理论闭合差
double startElev = points[pointIndex[knownPoints[i]]].elevation;
double endElev = points[pointIndex[knownPoints[j]]].elevation;
cond.w = (endElev - startElev) - (-cond.w);
conditions.push_back(cond);
}
}
}
}
}
// 2. 处理闭合环（改进算法）
std::unordered_set<
std::string> visitedNodes;
for (const auto& node
: graph) {
if (visitedNodes.find(node.first) != visitedNodes.end())
continue;
std::unordered_set<
std::string> currentComponent;
std::queue<
std::string> bfsQueue;
bfsQueue.push(node.first);
currentComponent.insert(node.first);
// BFS遍历连通分量
while (!bfsQueue.empty()) {
std::string current = bfsQueue.front();
bfsQueue.pop();
for (const auto& neighbor
: graph[current].neighbors) {
if (currentComponent.find(neighbor.first) == currentComponent.end()) {
currentComponent.insert(neighbor.first);
bfsQueue.push(neighbor.first);
}
}
}
visitedNodes.insert(currentComponent.begin(), currentComponent.end());
// 为每个连通分量构建闭合环
for (const auto& startNode
: currentComponent) {
for (const auto& neighbor
: graph[startNode].neighbors) {
std::vector<
std::vector<size_t>> paths = findAllPaths(neighbor.first, startNode, graph);
for (const auto& path
: paths) {
if (!path.empty()) {
Condition cond;
std::unordered_set<size_t> obsInPath;
for (size_t obsIdx
: path) {
if (obsInPath.find(obsIdx) == obsInPath.end()) {
cond.addObservation(obsIdx, 1.0);
obsInPath.insert(obsIdx);
}
}
if (!cond.obsIndices.empty()) {
cond.calculateW(observations);
conditions.push_back(cond);
}
}
}
}
}
}
// 确保条件方程个数正确
size_t necessaryObs = 0;
for (const auto& p
: points) {
if (!p.isFixed) necessaryObs++;
}
if (necessaryObs == 0) necessaryObs = points.size() - 1;
size_t r = observations.size() - necessaryObs;
while (conditions.size() > r) {
conditions.pop_back();
}
std::cout << "生成 " << conditions.size() << " 个独立条件方程" <<
std::endl;
}
void calculateAdjustedElevations() {
// 保存初始高程
std::vector<double> initialElevations;
for (const auto& p
: points) {
initialElevations.push_back(p.elevation);
}
// 使用改正后的观测值计算高程
for (const auto& obs
: observations) {
size_t fromIdx = pointIndex[obs.from];
size_t toIdx = pointIndex[obs.to];
if (points[fromIdx].isFixed && !points[toIdx].isFixed) {
// 已知点到未知点：H_to = H_from + dh_adjusted
points[toIdx].elevation = points[fromIdx].elevation + obs.adjustedValue;
}
else if (!points[fromIdx].isFixed && points[toIdx].isFixed) {
// 未知点到已知点：H_from = H_to - dh_adjusted
points[fromIdx].elevation = points[toIdx].elevation - obs.adjustedValue;
}
}
// 处理多个已知点的情况
bool changed;
do {
changed = false;
for (const auto& obs
: observations) {
size_t fromIdx = pointIndex[obs.from];
size_t toIdx = pointIndex[obs.to];
if (points[fromIdx].isFixed && !points[toIdx].isFixed) {
double newElev = points[fromIdx].elevation + obs.adjustedValue;
if (fabs(newElev - points[toIdx].elevation) > EPSILON) {
points[toIdx].elevation = newElev;
changed = true;
}
}
else if (!points[fromIdx].isFixed && points[toIdx].isFixed) {
double newElev = points[toIdx].elevation - obs.adjustedValue;
if (fabs(newElev - points[fromIdx].elevation) > EPSILON) {
points[fromIdx].elevation = newElev;
changed = true;
}
}
else if (!points[fromIdx].isFixed && !points[toIdx].isFixed) {
// 如果有一个点已计算，计算另一个点
if (points[fromIdx].elevation > EPSILON && points[toIdx].elevation < EPSILON) {
points[toIdx].elevation = points[fromIdx].elevation + obs.adjustedValue;
changed = true;
}
else if (points[toIdx].elevation > EPSILON && points[fromIdx].elevation < EPSILON) {
points[fromIdx].elevation = points[toIdx].elevation - obs.adjustedValue;
changed = true;
}
}
}
} while (changed);
}
void evaluatePrecision() {
pointPrecisions.resize(points.size(), 0.0);
double mu0 =
std::sqrt(unitWeightVariance) * 1000; // 单位毫米
// 正确的高程中误差计算
for (size_t i = 0; i < points.size(); ++i) {
if (!points[i].isFixed) {
double weightSum = 0.0;
int connectionCount = 0;
// 计算与点相连的所有观测值的权
之和
for (const auto& obs
: observations) {
if (obs.from == points[i].id || obs.to == points[i].id) {
weightSum += obs.weight;
connectionCount++;
}
}
if (weightSum > EPSILON) {
// 正确公式：μ_i = μ₀ / √(ΣP)
double mu = mu0 /
std::sqrt(weightSum);
pointPrecisions[i] = mu;
}
}
}
}
bool readDataFromFile(const
std::string& filename) {
// 显示尝试读取的文件路径
std::cout << "尝试读取文件
: " << filename <<
std::endl;
std::ifstream file(filename);
if (!file.is_open()) {
// 显示详细错误信息
std::cerr << "错误：无法打开文件 '" << filename << "'" <<
std::endl;
// 显示当前工作目录
#ifdef _WIN32
char buffer[MAX_PATH];
if (GetCurrentDirectoryA(MAX_PATH, buffer)) {
std::cerr << "当前工作目录
: " << buffer <<
std::endl;
}
#else
char buffer[1024];
if (getcwd(buffer, sizeof(buffer)) != nullptr) {
std::cerr << "当前工作目录
: " << buffer <<
std::endl;
}
#endif
return false;
}
std::string line;
points.clear();
observations.clear();
conditions.clear();
pointIndex.clear();
if (!
std::getline(file, line)) {
std::cerr << "错误：文件格式不正确，无法读取点数" <<
std::endl;
file.close();
return false;
}
int numPoints;
try {
numPoints =
std::stoi(line);
}
catch (const
std::exception& e) {
std::cerr << "错误：无效的点数格式
: " << e.what() <<
std::endl;
file.close();
return false;
}
std::cout << "读取 " << numPoints << " 个点数据
:" <<
std::endl;
for (int i = 0; i < numPoints; ++i) {
if (!
std::getline(file, line)) {
std::cerr << "错误：文件格式不正确，点数据读取失败" <<
std::endl;
file.close();
return false;
}
std::istringstream iss(line);
std::string id;
double elevation;
int isFixedInt;
if (!(iss >> id >> elevation >> isFixedInt)) {
std::cerr << "错误：点数据格式不正确
: " << line <<
std::endl;
file.close();
return false;
}
addPoint(Point(id, elevation, isFixedInt != 0));
std::cout << " " << id << "
: " << elevation << "m ("
<< (isFixedInt != 0 ? "已知点"
: "待定点") << ")" <<
std::endl;
}
// 读取观测值数据
if (!
std::getline(file, line)) {
std::cerr << "错误：文件格式不正确，无法读取观测值数" <<
std::endl;
file.close();
return false;
}
int numObs;
try {
numObs =
std::stoi(line);
}
catch (const
std::exception& e) {
std::cerr << "错误：无效的观测值数格式
: " << e.what() <<
std::endl;
file.close();
return false;
}
std::cout << "读取 " << numObs << " 个观测值数据
:" <<
std::endl;
for (int i = 0; i < numObs; ++i) {
if (!
std::getline(file, line)) {
std::cerr << "错误：文件格式不正确，观测值数据读取失败" <<
std::endl;
file.close();
return false;
}
std::istringstream iss(line);
std::string from, to;
double dh, distance;
if (!(iss >> from >> to >> dh >> distance)) {
std::cerr << "错误：观测值数据格式不正确
: " << line <<
std::endl;
file.close();
return false;
}
double weight = 1.0 / distance;
addObservation(Observation(from, to, dh, weight));
std::cout << " " << from << " -> " << to
<< "
: " << dh << "m, 距离
: " << distance
<< " 权
: " << weight <<
std::endl;
}
// 验证观测值引用的点是否存在
for (const auto& obs
: observations) {
if (pointIndex.find(obs.from) == pointIndex.end()) {
std::cerr << "错误：观测值起点 '" << obs.from << "' 不存在" <<
std::endl;
return false;
}
if (pointIndex.find(obs.to) == pointIndex.end()) {
std::cerr << "错误：观测值终点 '" << obs.to << "' 不存在" <<
std::endl;
return false;
}
}
file.close();
return true;
}
void adjust() {
// 保存初始高程
initialElevations.clear();
for (const auto& p
: points) {
initialElevations.push_back(p.elevation);
}
std::cout << "\n开始水准网平差计算..." <<
std::endl;
if (points.empty() || observations.empty()) {
std::cerr << "错误：数据不完整" <<
std::endl;
return;
}
// 识别条件方程
identifyConditions();
size_t n = observations.size();
size_t r = conditions.size();
if (r == 0) {
std::cerr << "错误：未生成条件方程" <<
std::endl;
return;
}
// 构建条件方程矩阵
Matrix A(r, n);
for (size_t i = 0; i < r; ++i) {
const Condition& cond = conditions[i];
for (size_t j = 0; j < cond.obsIndices.size(); ++j) {
size_t obsIdx = cond.obsIndices[j];
if (obsIdx < n) {
A(i, obsIdx) = cond.coefficients[j];
}
}
}
// 输出条件方程矩阵
std::cout << "\n条件方程矩阵 A
:" <<
std::endl;
for (size_t i = 0; i < r; ++i) {
for (size_t j = 0; j < n; ++j) {
std::cout <<
std::fixed <<
std::setprecision(6) << A(i, j) << "\t";
}
std::cout << "| " << conditions[i].w <<
std::endl;
}
// 构建权对角向量
std::vector<double> P_diag(n);
for (size_t i = 0; i < n; ++i) {
P_diag[i] = observations[i].weight;
}
// 构建法方程 N = A * diag(P) * A^T
Matrix N(r, r);
for (size_t i = 0; i < r; ++i) {
for (size_t j = 0; j < r; ++j) {
double sum = 0.0;
for (size_t k = 0; k < n; ++k) {
sum += A(i, k) * P_diag[k] * A(j, k);
}
N(i, j) = sum;
}
}
// 构建闭合差向量
std::vector<double> w(r);
for (size_t i = 0; i < r; ++i) {
w[i] = conditions[i].w;
}
// 解法方程 N * k = -w
std::vector<double> minusW(r);
for (size_t i = 0; i < r; ++i) {
minusW[i] = -w[i];
}
std::vector<double> k;
try {
k = Matrix
::solve(N, minusW);
}
catch (const
std::exception& e) {
std::cerr << "法方程
求解失败
: " << e.what() <<
std::endl;
return;
}
// 输出联系数向量
std::cout << "\n联系数向量 k
:" <<
std::endl;
for (size_t i = 0; i < r; ++i) {
std::cout << "k[" << i << "] = " <<
std::fixed <<
std::setprecision(6) << k[i] <<
std::endl;
}
// 计算改正数 v = P^{-1} * A^T * k
std::vector<double> v(n, 0.0);
for (size_t i = 0; i < n; ++i) {
for (size_t j = 0; j < r; ++j) {
v[i] += A(j, i) * k[j] / P_diag[i];
}
observations[i].correction = v[i];
observations[i].adjustedValue = observations[i].dh + v[i];
// 输出改正数信息
std::cout << "观测值 " << i + 1 << " (" << observations[i].from << "->" << observations[i].to
<< ")
: 改正数 = " <<
std::fixed <<
std::setprecision(2) << v[i] /1000.0 << " mm"
<< ", 平差值 = " << observations[i].adjustedValue << " m" <<
std::endl;
}
// 计算单位权方差
double sumPVV = 0.0;
for (size_t i = 0; i < n; ++i) {
sumPVV += v[i] * v[i] * P_diag[i];
}
unitWeightVariance = sumPVV / r;
double mu0 =
std::sqrt(unitWeightVariance);
std::cout << "单位权中误差
: ±" <<
std::fixed <<
std::setprecision(2) << mu0 * 1000.0 << " mm" <<
std::endl;
// 计算高程平差值
calculateAdjustedElevations();
// 精度评定
evaluatePrecision();
// 输出结果
printResults();
}
// 输出计算结果
void printResults() const {
std::cout << "\n===== 水准网平差计算结果 =====" <<
std::endl;
// 1. 待定点高程平差值
std::cout << "\n1. 待定点高程平差值
:" <<
std::endl;
std::cout <<
std::left <<
std::setw(10) << "点号"
<<
std::setw(15) << "初始高程(m)"
<<
std::setw(15) << "平差后高程(m)"
<<
std::setw(15) << "中误差(mm)" <<
std::endl;
std::cout << "--------------------------------------------------------" <<
std::endl;
for (size_t i = 0; i < points.size(); ++i) {
const auto& p = points[i];
std::cout <<
std::left <<
std::setw(10) << p.id
<<
std::setw(15) <<
std::fixed <<
std::setprecision(6) << (p.isFixed ? p.elevation
: 0.0)
<<
std::setw(15) <<
std::fixed <<
std::setprecision(6) << p.elevation;
if (!p.isFixed) {
std::cout <<
std::setw(15) <<
std::fixed <<
std::setprecision(2) << pointPrecisions[i];
}
std::cout <<
std::endl;
}
// 观测值改正数（以毫米为单位）
for (size_t i = 0; i < observations.size(); ++i) {
const auto& obs = observations[i];
std::cout <<
std::left <<
std::setw(10) << (i + 1)
<<
std::setw(10) << obs.from
<<
std::setw(10) << obs.to
<<
std::setw(15) <<
std::fixed <<
std::setprecision(6) << obs.dh
<<
std::setw(15) <<
std::fixed <<
std::setprecision(2) << obs.correction /1000.0 // 毫米
<<
std::setw(15) <<
std::fixed <<
std::setprecision(6) << obs.adjustedValue <<
std::endl;
}
// 高程变化
std::cout << "\n3. 高程变化
:" <<
std::endl;
std::cout <<
std::left <<
std::setw(10) << "点号"
<<
std::setw(15) << "初始高程(m)"
<<
std::setw(15) << "平差后高程(m)"
<<
std::setw(15) << "变化量(mm)" <<
std::endl;
std::cout << "--------------------------------------------------------" <<
std::endl;
for (size_t i = 0; i < points.size(); ++i) {
const auto& p = points[i];
double change = (p.elevation - initialElevations[i]) / 1000.0; // 毫米
std::cout <<
std::left <<
std::setw(10) << p.id
<<
std::setw(15) <<
std::fixed <<
std::setprecision(6) << initialElevations[i]
<<
std::setw(15) <<
std::fixed <<
std::setprecision(6) << p.elevation
<<
std::setw(15) <<
std::fixed <<
std::setprecision(2) << change <<
std::endl;
}
// 3. 精度评定
std::cout << "\n3. 精度评定
:" <<
std::endl;
std::cout << "单位权中误差
: ±" <<
std::fixed <<
std::setprecision(2)
<<
std::sqrt(unitWeightVariance) / 1000 << " mm" <<
std::endl;
}
};
int main() {
try {
LevelingNetwork network;
if (!network.readDataFromFile("leveling.txt")) {
std::cerr << "错误：无法读取数据文件。程序将退出。" <<
std::endl;
return 1;
}
network.adjust();
return 0;
}
catch (const
std::exception& e) {
std::cerr << "错误
: " << e.what() <<
std::endl;
return 1;
}
}这是我写的程序，假如你是一个资深的测量平差以及c++的老师，面对我的水准网条件平差的程序，请帮我修改一下，使得该程序不仅能实现水准网平差过程通用程序，包括误差方程、法方程的组成与解算。得出平差后的高差平差值及各待定点的高程平差值；评定各待定点高程平差值的精度、观测差平差值的精度；；根据水准点间的相对关系，假定坐标，利用程序画出水准网示意图；还能够识别A-B，B-C，A-C，还能使得列出的所有条件方程足数（观测总数减去必要观测个数等于条件方程个数，如果有足够的已知观测数据，必要观测个数等于未知数个数；如果没有足够已知数据，必要观测个数等于总点数-1）且并不相关，不是奇异矩阵。（要
求不借助外部数据库，自己编写矩阵以及其他的函数），是在visual studio2022上实现，且只是c++11或14中。且如果是不需要的请帮我简化一下。重点检查一下改正数，单位权中误差，改正后高差，高程中误差的计算以及寻找闭合和符合导线是否有问题，请帮我修正
本文通过解决C++primer习题时遇到的问题，强调了在主函数调用函数前进行函数定义或声明的重要性。文章分享了一个实例，展示了未遵循此规则导致的调试困难及解决方案。
扫一扫
562
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
