数据结构与测绘软件开发
课程设计
边角网测量平差程序设计
| 姓 名： | FXXX |
| 班 级： | 测绘工程XXXX班 |
| 学 号： | XXXXXX |
| 指导老师： | XXXX |
XXXXXXXX)XXX
2021年12月
前言：主要分享了一个基于C++文本输入输出流读取边角网数据，进行测量平差的程序。限于本人水平有限，望读者恳请指正。
目录
一、课设概况
1.1课设目标
课程设计的任务是进行边角网测量平差程序设计，旨在通过C++文件流数据读写的方式，并进行边角网平差计算，以实现数据的批量输入，经过程序化处理后输出的目标。
1.2开发环境
本次课设中我的应用程序开发框架选择的是：开源版、Windows平台下的Qt5.9.0。
集成开发环境（IDE）则为Qt5.9.0相对应的Qt Creator4.3.0。开发中我选择创建简单的控制台应用程序。
除了开发环境的选择外，还引用了Eigen3.4.0。Eigen库是一个开源的C++模板库，在这次实习中主要是用来进行矩阵运算和数值分析。
Eigen目前最新的版本是3.4.0，除了C++标准库以外，不需要任何其他的依赖包。Eigen使用的CMake建立配置文件和单元测试，并自动安装。如果使用Eigen库，只需包特定模块的的头文件即可，简单了说它就是一个C++版本的matlab包。一般情况下，使用include<Eigen/Dense>即可。
表1 Eigen库模块表
| Module | Header file | Contents |
| Core | #include<Eigen/Core> | Matrix和Array类，基础的线性代数运算和数组操作 |
| Geometry | #include<Eigen/Geometry> | 旋转、平移、缩放、2维和3维的各种变换 |
| LU | #include<Eigen/LU> | 求逆，行列式，LU分解 |
| Cholesky | #include <Eigen/Cholesky> | LLT和LDLT Cholesky分解 |
| Householder | #include<Eigen/Householder> | 豪斯霍尔德变换，用于线性代数运算 |
| SVD | #include<Eigen/SVD> | SVD分解 |
| QR | #include<Eigen/QR> | QR分解 |
| Eigenvalues | #include<Eigen/Eigenvalues> | 特征值，特征向量分解 |
| Sparse | #include<Eigen/Sparse> | 稀疏矩阵的存储和一些基本的线性运算 |
| 稠密矩阵 | #include<Eigen/Dense> | Core/Geometry/LU/Cholesky/SVD/QR/Eigenvalues模块 |
| 矩阵 | #include<Eigen/Eigen> | 包括Dense和Sparse(整合库) |
不同的IDE下，引用Eigen库方式不同，Qt Creator在项目工程文件中添加Eigen的路径，添加一行代码INCLUDEPATH += Eigen库的路径即可，如下图所示：
图1 Qt环境下添加Eigen库
至于其他环境，如VS，VC++等，在此不再说明，面向百度即可，在整个开发过程中，我遇到很多不懂的方法，以及bug等等，也是通过查看书籍、面向“百度”、以及向的优秀同学们请教等最终完成了这次边角网数据平差的课设目标。
二、平差模型
本程序采用间接平差模型。令未知点的坐标平差值为参数。
平差有两种方式，一种是条件平差，一种是间接平差。条件方程式较少，占用内存小，但条件方程与网形有关，很难编制通用的程序。间接平差一个观测值就是一个误差方程，占用内存相对较大，但易于编程解算，容易编制通用的程序。从程序编制和程序的通用性考虑，所以选择间接平差的方法。
2.1误差方程
（1）边长误差方程：
设Ski为k号点至i号点的边长观测值，则有：
线性化：
（2）观测角度误差方程
2.2平差原理
（1）将每一个观测量的平差值分别表达成所选参数的函数：
（2）若函数非线性要将其线性化，列出误差方程
（3）由误差方程系数B和自由项组成法方程：
2.3精度评定
（1）单位权中误差
（2）坐标平差值的协因数阵：
QXX=(ATPA)-1
（3）K点坐标平差值的中误差、点位中误差：
2.4误差椭圆参数计算
计算出两个根后，直接带入公式计算E和F，再判断E和F谁是极大值，极小值：
选择φ方向位差公式，因为用它程序计算方便直观，别的方式需要进行其他的判断。
三、数据文件格式
将起算数据、观测精度指标、观测值及相关数据存入Data.txt文本文档中，进行数据读取操作，因此将涉及设计他们之间的数据关系以及结构设计。以如下图2边角网为例，说明数据文件的格式。
已知点为A、B、C、D、E、F共6个，未知点编号7，8,9,10·······17，18共12个，观测了14条边长S1～S14，观测了16个转角β1～β16，测角精度10″，测边精度10mm+15ppm。
图2 示例边角网
1.边角网概况信息
此部分，说明边角网的整体控制信息，包含，已知点数目、未知点数目、边长观测数、角度观测数，在本例中，它们分别为6、12、14、16。
以及精度指标：角度观测精度（秒），光电测距边长精度Amm+Bppm*L的固定误差A，比例误差B。本例中，分别为10、10、15。
在Data.txt文档中的中文描述语要为一整句，即中间须有标点符号连接，不能有空格，否则数据读取时会出错，对于数据的数字之间则用一个空格隔开，如下图3所示：
图3 边角网概况信息数据格式
2.已知点数据格式
已知点的ID编号以及X坐标和Y坐标。ID编号必须为数字，且从1开始编号，不能为字母。在本程序中，设计了一个大小为 ：已知点数*3的矩阵（数组），第一列是点的ID编号、第二列是X坐标，第三列是Y坐标。因为统一放在同一个矩阵，所以ID编号必须为数字，当然也可以用类型不同的数组存放，字符型的或是其他类型，但是在我程序编写的过程中，我自认为那样变得更加繁琐冗杂了，另外为了从存放已知点信息的这个矩阵调用已知点数据更加方便，所以从1开始编号。
本例中已知点为A、B、C、D、E、F分别令他们ID编号为1、2、3、4、5、6。一行输入已知点的ID编号以及X坐标和Y坐标，之间用空格隔开，如下图4所示：
图4 已知点数据格式
3.未知点数据格式
同已知点格式类似，但是只需要ID编号即可，未知点的ID编号和已知点连续，本例中，已知点有6个，那未知点的ID编号就从7开始。
本例中未知点,分别令他们ID编号为7～18，如下图5所示：
图5 未知点数据格式
4.边长观测数据格式
边长数据格式设计为：边长观测数 * 4的矩阵，每一行数据存储边长的ID编号，从1开始；对应的边长，单位：米/m；起点、终点的ID编号，这里ID编号即是已知点或者未知点的对应编号ID。
本例中，观测14条边，数据格式如下图6所示：
图6 边长观测数据格式
5.角度观测数据格式
角度数据格式设计为：角度观测数 * 5的矩阵，每一行数据存储角度观测值的测站ID，后视照准ID_1，前视照准ID_2，夹角观测值，测站ID和前视照准ID_2对应的边的ID编号。
其中夹角按照顺时针记录，数字用“度.分秒”表示，如23.20205格式表示23度20分20.5秒。测站ID，后视照准ID_1，前视照准ID_2指的是当前对应的已知点或未知点的ID编号。测站ID和前视照准ID_2对应的边的ID编号，这里边的ID编号指的是测站为起点（终点），前视照准为终点（起点）的那条边的ID编号。
参考图2示例边角网，以本例中的夹角β3具体说明：
测站ID=8
后视照准ID_1=7
前视照准ID_2=9
角值=250.1811
测站ID和前视照准ID_2对应的边的ID编号=3。其起点（终点）为8，终点（起点）为9，对应的边是S3，在边长数据格式中该边的ID编号是3。若测站ID和前视照准ID_2对应的边是两个已知点之间，那对应边ID此时写为0即可。
另外，该数据编写时，尽量从一已知方向出发，编写完该方向的所有观测角数据，再换另外的方向，否则，数据不清晰虽然并不影响程序执行，但是影响可读性。
本例角度观测数据格式如下图7所示：
图7 角度观测数据格式
通过以上的数据结构设计，将所有数据存放到Data.txt，并将Data.txt存储在编译环境下的工程项目中，以方便C++文件流以相对路径打开文件读写，当然也可以不考虑Data.txt的文件位置，在文件读取时以绝对路径即可。
但我认为，相对路径的方式更好，绝对路径用户的电脑命名不同，可能导致数据读取出错，另外，文件、程序的路径尽量不要包含有中文，横杠等特殊字符，可能导致不同的电脑不能运行。不同的IDE存放相对路径目录不同，在此仅以Qt为例，在编译器构建生成的同级目录下。如下图8所示：
图8 Date.txt相对路径位置
本例的数据格式全貌如下图9所示：
图9 Date.txt数据
四、类及功能设计
本程序我设计了读取Data.txt文本文件的类InData、角度格式转换的类AngleChange、计算近似坐标的类X0Y0、平差模块的Ajust、将平差处理后的相关数据输出的类OutData。它们之间的关系如下所示：
4.1 读取Data.txt数据类InData
类体定义在indta.h头文件，声明函数具体实现indata.cpp文件。源代码见附录。在此说明定义了什么变量，函数以及它们关系。
（1）indata.h
class InData//边角网数据读取类
{
public://定义一些外部接口成员
int knownpoint,unknownpoint,Num_Distance,Num_Angle;
//已知点数，未知点数，边长观测数，角度观测数
double Pre_Angle,Pre_Distance_A,Pre_Distance_B;
//角度观测精度（秒），光电测距边长精度Amm+Bppm*L的固定误差A，比例误差B
Eigen::MatrixXd KnownPointXY;//定义矩阵存放已知点的信息
Eigen::MatrixXd UnKnownPointX0Y0;//定义矩阵存放未知点的信息
Eigen::MatrixXd Distance_Obs;//定义矩阵存放边长观测值的信息
Eigen::MatrixXd Angle_Obs;//定义矩阵存放角度观测值的信息
public:
void DoInData();//从Data.txt读取信息，存储到InData类
};
（2）indata.cpp
在此处具体实现函数void DoInData()的功能，在mian函数时调用该函数，以实现数据读取。
定义文件流对象FileRead，用于读取项目下的Data.txt文件。
ifstream FileRead;
FileRead.open("Data.txt",ios::in);//读取项目相对路径下的Data.txt文件
具体代码间附录。可以将读到的数据输出到控制台应用，即显示屏上，结果如下图10。并且可供检查Data.txt数据是否正确。
图10 读取Data数据存放到成员
4.2角度格式转换类AngleChange
由于程序处理运算使用的三角函数等是弧度值，因此需要将Data.txt数据中的度.分秒转换为0～2π的弧度值。AngleChange需要使用InData中读取到的数据，因此InData类的成员属性定义为public，并且AngleChange类public继承了InData类。
class AngleChange:public InData
{
public:
Eigen::MatrixXd Angle_Rad;//定义存放角度转换为弧度后的矩阵
void DoAngleChange();//弧度转换
};具体代码见附录。
4.3计算未知点近似坐标类X0Y0
利用读取到的数据计算近似坐标计算，并存储到新定义的矩阵AllPointX0Y0，存放已知点坐标和未知点计算出来后的近似坐标。同样，X0Y0类public继承了AngleChange类。类体如下：
#include "anglechange.h"
class X0Y0:public AngleChange
{
public:
//存放已知点坐标和未知点计算出来后的近似坐标
Eigen::MatrixXd AllPointX0Y0;
void DoX0Y0();//计算未知点近似坐标
};
以下图11中未知点7计算近似坐标说明：
图11 未知点7计算示例图
（1）在角β1的观测数据中，即矩阵Angle_Obs中可知：
测站ID=2，后视ID_1=1 ，前视ID_2=7 ，夹角值（弧度）=163.4504（2.857996043），边长ID=1
（2）判断前视ID_2是否为未知点ID。
//检索并判断当前的未知点的X坐标是否为0；
if(AllPointX0Y0(i+knownpoint,1)!=0)
continue;
//若不等于0，说明上一次检索已经计算过当前未知点X和Y坐标，就不再执行下面的语句
int Station_ID=Angle_Rad(j,0);//测站ID
int Station_ID1=Angle_Rad(j,1);//后照ID_1
double Angle=Angle_Rad(j,3);//角度观测值（弧度）
int Dis_ID=Angle_Rad(j,4);//测站ID和前照ID_2对应的边长编号ID
if(Angle_Rad(j,2)==UnKnownPointX0Y0(i,0)&&Dis_ID!=0)//检索并判断矩阵Angle_Rad第3列数据，即前照ID_2是否等于当前要计算的未知点ID编号，且对应的边长编号ID不为0，若条件成立，通过矩阵Angle_Rad中该前照ID_2对应的测站ID，后照ID_1计算未知点近似X和Y坐标
（3）该示例，前视ID_2=7是未知点ID编号，就执行计算。
double X1,Y1,X2,Y2;
X1=AllPointX0Y0(Station_ID1-1,1);//测站后照ID_1对应X坐标
Y1=AllPointX0Y0(Station_ID1-1,2);//测站后照ID_1对应Y坐标
X2=AllPointX0Y0(Station_ID-1,1);//测站ID对应X坐标
Y2=AllPointX0Y0(Station_ID-1,2);//测站ID对应Y坐标
double a0=atan2((Y2-Y1),(X2-X1));//计算上一条边方位角
if(a0<0)//使方位角全为0~2PI范围
{
a0=a0+2*PI;}
double a=a0+Angle-PI;//计算测站ID和前照ID_2对应的边的方位角
double X0,Y0;//计算当前未知点的对应近似X0和Y0坐标
X0=X2+Distance_Obs(Dis_ID-1,1)*cos(a);
Y0=Y2+Distance_Obs(Dis_ID-1,1)*sin(a);
（4）计算完未知点7坐标，更新矩阵AllPointX0Y0
之后的未知点计算需要用到之前的未知点坐标，如接下来未知点8的近似坐标计算，就需要未知点7的近似坐标，所以将计算出来的近似坐标存放回矩阵AllPointX0Y0。只需要如下简单的两句语句即可：
AllPointX0Y0(i+knownpoint,1)=X0;
AllPointX0Y0(i+knownpoint,2)=Y0;
该示例计算成功后的结果如下图12：其中1～6是已知点数据。
图12 近似坐标计算成果
4.4平差模块类Ajust
计算出未知点近似坐标后，就可以误差方程计算系数阵B，自由项等，平差模块类Ajust主要定义了平差所需的相关矩阵和变量，以及计算的函数，同样，需要用到之前类的数据，所以Ajust类public继承了X0Y0类。
类体如下：
#include "x0y0.h"
class Ajust:public X0Y0
{public:
Eigen::MatrixXd V;//改正数矩阵
Eigen::MatrixXd B;//误差方程系数矩阵
Eigen::MatrixXd x;//待定参数矩阵
Eigen::MatrixXd l;//自由项矩阵
Eigen::MatrixXd P;//权矩阵
Eigen::MatrixXd Nbb;//法方程系数阵
double sigma0;//单位权中误差
Eigen::MatrixXd QXX;//平差参数X的协因数阵
Eigen::MatrixXd DXX;//平差参数X的协方差阵
Eigen::MatrixXd Fhi0_E;//极大值方向
Eigen::MatrixXd Fhi0_F;//极大值方向
Eigen::MatrixXd E;//误差椭圆极大值
Eigen::MatrixXd F;//误差椭圆极小值
public:
void DoForBl();//计算误差方程系数阵B、自由项矩阵l
void DoForP();//计算权矩阵P
void DoForVx();//计算改正数矩阵V、待定参数矩阵x
void Accuracy_assess();//精度评定函数
void Error_ellipse();//误差椭圆参数计算
};具体实现见附录源代码。
4.5数据输出至PutData.txt类OutData
此类主要定义了函数void DoPutData()进行相关数据及其成果的输出，通过文件流输出至项目相对路径下的文本文档PutData.txt中，因为要使用到之前数据成果，所以OutData类public继承了Ajust类。类体如下
#include "ajust.h"
class OutData:public Ajust
{
public:
void DoPutData();//将数据输出至PutData.txt文本文档
};具体实现见附录源代码。
五、运算实例
此部分选择《测量平差》-第二版-中国矿业大学出版社，教材中几个例题，按照本程序的思路，编写Data.txt文本文档，就程序运行进行实例展示。
例题一：课本103页，例【4-2】；
例题二：课本139页，例【5-3】；
例题三：课本245页，第四题。
5.1例题一
（1）此题为，课本103页，例【4-2】，补充说明：
测边精度：0mm+15ppm;详情如下图13
图13 例【4-2】图
（2）编写Data.txt数据。
图14 例【4-2】Data.txt数据
（3）将数据文档命名为“Data.txt”放置于平差程序项目路径下。然后开始运行。运行结果如下图15：
图15 控制台运行窗口结果
（4）提示运行成功后，查看项目路径下PutData.txt文件，查看平差结果。本例题结果如下图16。
图16 例【4-2】平差数据PutData.txt
参考课本上结果，平差结果成果，且数据比对相同，正确。
5.2例题二
（1）此题为，课本139页，例【5-3】。
网形详情如下图17：
图17 例题二网形图
（2）编写Data.txt数据
图18 例题二Data.txt数据
（3）将数据文档命名为“Data.txt”放置于平差程序项目路径下。然后开始运行。运行结果如下图19：
图19 运行结果
（4）提示运行成功后，查看项目路径下PutData.txt文件，查看平差结果。本例题二结果如下图20。
图20 例题二平差数据PutData.txt
参考课本上结果，平差结果成果，且数据比对相同，正确。
5.3例题三
（1）此题为课本245页，第四题。
图21 例题三网型
已知点为A、B、C、D、E、F共6个，未知点编号7，8,9,10·······17，18共12个，观测了14条边长S1～S14，观测了16个转角β1～β16，测角精度10″，测边精度1.0*
/mm。
可见，此题的测边精度指标不是光电测距仪精度指标Amm+Bppm形式;仅需要在程序中简单修改以下权矩阵P中的边长的权即可。
以下是测边精度指标为Amm+Bppm时的计算P矩阵代码：
void Ajust::DoForP()//定权矩阵P
{
double P_distance;//观测边的权
double P_angle;//观测角的权
//设置权阵P大小并初始化为0
P.setZero(Num_Distance+Num_Angle,Num_Distance+Num_Angle);
for(int i=0;i<Num_Distance;i++)
{
//对应边中误差
double delta=Pre_Distance_A+Pre_Distance_B*0.001*Distance_Obs(i,1);
P_distance=pow(Pre_Angle,2)/pow(delta,2);
P(i,i)=P_distance;
}
for(int j=0;j<Num_Angle;j++)
{
P_angle=1;//取角度观测精度Pre_Angle为验前单位权中误差
P(j+Num_Distance,j+Num_Distance)=P_angle;
}
}
测边精度 1.0* /mm，只需要更改定权时对应边中误差的计算方法即可。
将原来计算语句：
double delta=Pre_Distance_A+Pre_Distance_B*0.001*Distance_Obs(i,1)；
改为新的计算语句：
double delta=sqrt(Distance_Obs(i,1))；
当测边精度指标为 1.0* /mm时，权矩阵P的计算代码更改后如下：
void Ajust::DoForP()//定权矩阵P
{
double P_distance;//观测边的权
double P_angle;//观测角的权
//设置权阵P大小并初始化为0
P.setZero(Num_Distance+Num_Angle,Num_Distance+Num_Angle);
for(int i=0;i<Num_Distance;i++)
{
//对应边中误差
double delta=sqrt（Distance_Obs(i,1)）;
P_distance=pow(Pre_Angle,2)/pow(delta,2);
P(i,i)=P_distance;
}
for(int j=0;j<Num_Angle;j++)
{
P_angle=1;//取角度观测精度Pre_Angle为验前单位权中误差
P(j+Num_Distance,j+Num_Distance)=P_angle;
}
}
（2）编写Data.txt数据
图22 例题三Data.txt数据
（3）将数据文档命名为“Data.txt”放置于平差程序项目路径下。然后开始运行。运行结果如下图23：
图23 例题三运行结果
（4）提示运行成功后，查看项目路径下PutData.txt文件，查看平差结果。本例题二结果如下图24。
图24 例题三平差PutData.txt数据
六、总结
6.1开发日志
6.2感想
理论出真理，实践长才干。虽然这次课设目标勉强算是完成了，但这并不是学习的止步。通过这次课设，重拾了编程，体会到了程序化应用的优势，之后也会继续更深入的学习。
七、附录（源代码）
7.1项目工程管理文件day1.pro
QT += core
QT -= gui
INCLUDEPATH +=F:Eigen3.4.0
CONFIG += c++11
TARGET = day1
CONFIG += console
CONFIG -= app_bundle
TEMPLATE = app
SOURCES += main.cpp
anglechange.cpp
x0y0.cpp
ajust.cpp
indata.cpp
outdata.cpp
# The following define makes your compiler emit warnings if you use
# any feature of Qt which as been marked deprecated (the exact warnings
# depend on your compiler). Please consult the documentation of the
# deprecated API in order to know how to port your code away from it.
DEFINES += QT_DEPRECATED_WARNINGS
# You can also make your code fail to compile if you use deprecated APIs.
# In order to do so, uncomment the following line.
# You can also select to disable deprecated APIs only up to a certain version of Qt.
#DEFINES += QT_DISABLE_DEPRECATED_BEFORE=0x060000 # disables all the APIs deprecated before Qt 6.0.0
HEADERS +=
anglechange.h
x0y0.h
ajust.h
indata.h
outdata.h
。具体代码及配置调试可留言交流。
原来的QT C++版本的源码不在了，我后面用C#新写了一版本，逻辑差不多，源码资源在这里(https://download.csdn.net/download/From_fmyjj/89055841)
552
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
