|
| 中国矿业大学 |
| China University of Mining and Technology |
数据结构与测绘软件开发
边角网测量平差程序设计
| 姓? ???????名： | FXXX |
| 班? ?????级： | 测绘工程XXXX班 |
| 学 ????????号： | XXXXXX |
| 指导老师： | XXXX |
中国矿业大学环境与测绘学院
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开发日志
本次整个程序设计从动员会12月07号开始，听完老师的讲解后，限于迷茫，没有清晰地思路，所以没有立刻开始。在思考几天之后，便开始了从学习C++文件流开始，以下是我的课设学习及开发日志。
12月7号～10号。这几天在思考思路，看书，从C++的基本语法复习，着重学习了C++文件流读取txt文件。理清楚了文件读取时“<<”、“>>”的用法。
12月11号。思考采用什么数据结构，复习“数据结构”知识，上网自学等。在百度“C++矩阵运算”时，Eigen映入我的眼帘。之后便学习了Eigen的用法。
12月12号。Eigen学习完了，思考怎么进行数据的平差，遍历，循环查找等问题。学习了基于Graph的算法等问题，但是觉得都是数学运算，于是想可不可以要都放在矩阵里存储，用我自己设计矩阵对应的数据的查询关系应该也就可以计算了。
12月13号。初步设计了类结构，对读取Data.txt数据进行设计。当中读取数据没遇到问题，但是控制台和数据输出至txt时中文出现了乱码，我又面向“百度”，花了好几个小时。但是按照网上教程并不能解决我的问题，于是就请教了一位同学。
12月14号。问题解决后，我开始了计算类的设计。复习间接平差知识，这天完成了未知点近似坐标的计算。
12月15号。这天完成了系数阵B，自由项l的计算。
12月16号。今天完成了精度评定以及误差椭圆的计算，设置输出PutData.txt文件，但是遇到了输出的数据对齐整不美观，当时采用的是数据后面跟着输出制表符‘\t’,空格，来控制对齐等格式，但是差强人意，因为制表符‘\t’是占8个空格，如果数据宽度大于8，就会出现对不齐的问题。
12月17号～25号。这几天感觉主要问题都完成了，就休息了。
12月25号。准备开始写课设报告。于是就试着回去解决输出对齐等问题，贯彻面向“百度”，查资料，学习等。找到了#include <iomanip>，然后学习了一下用法，特别是填充setfill()函数，以及设置域宽setw()，设置有效数字位数setprecision()，固定小数位数fixed()等函数，通过尝试练习，我完成了输出相关数据对齐，以及数据的小数位数等问题。
12月26～28号，撰写报告。
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
7.2主函数入口main.cpp
#include <Eigen/Dense>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <string>
#include <cmath>
#include "indata.h"
#include "anglechange.h"
#include "x0y0.h"
#include "ajust.h"
#include "outData.h"
using namespace Eigen;
using namespace std;
int main()
{
//设置输出数字有效位数,不影响计算精度，主要是方便查看输出数据以便检查
cout.precision(10);
OutData get;//定义OutData类的对象get,调用OutData类及父类成员
??? get.DoInData();//从Data.txt读取信息，存储到InData类
??? cout<<"··························从Data.txt读取数据已完成##########"<<endl;
??? get.DoAngleChange();//弧度转换cout<<"································角度格式转换已完成##########"<<endl;
??? get.DoX0Y0();//计算未知点近似坐标
??? cout<<"···························未知点近似坐标计算已完成##########"<<endl;
??? get.DoForBl();//计算误差方程系数阵B、自由项矩阵l
??? cout<<"····················误差方程系数阵B、自由项矩阵l计算已完成#####"<<endl;
??? get.DoForP();//计算权矩阵P
??? cout<<"····························权矩阵P计算已完成##############"<<endl;
??? get.DoForVx();//计算改正数V、待求参数x
??? cout<<"························改正数V、待求参数x计算已完成#########"<<endl;
get.Accuracy_assess();//精度评定计算
cout<<"·····························精度评定计算已完成#############"<<endl;
??? get.Error_ellipse();//计算误差椭圆参数
??? cout<<"··························误差椭圆参数计算已完成#############"<<endl;
??? get.DoPutData();//输出数据到文件PutData.txt
??? cout<<endl;
??? cout<<"本次平差已结束，结果请查看项目下PutData.txt文件"<<endl;
??? system("pause");
??? return 0;
}
7.3类的.h和.cpp文件
1.读取Data.txt类InData
（1）indata.h
#ifndef INDATA_H
#define INDATA_H
#include <Eigen/Dense> //Eigen是一个开源的C++模板库，引用用于数值分析和矩阵运算
#include <iostream>
#include <fstream>
#include <string>
#include <cmath>
using namespace Eigen;
using namespace std;
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类
};
#endif // INDATA_H
（2）indata.cpp
#include "indata.h"
void InData::DoInData()
{
??? ifstream FileRead;
??? FileRead.open("Data.txt",ios::in);//读取项目相对路径下的Data.txt文件
??? if(!FileRead)
??? {
??????? cerr<<"文件读取失败，请检查数据和路径!"<<endl;
??????? exit(1);
??? }
??? string help_str;//定义一个字符串，无实际作用，主要用来读到中文描述语时跳过。
??? FileRead>>help_str;//跳过第1行中文描述语
??? FileRead>>knownpoint;
??? FileRead>>unknownpoint;
??? FileRead>>Num_Distance;
??? FileRead>>Num_Angle;
??? cout<<"边角网数据概况:"<<endl;
??? cout<<"已知点数:"<<knownpoint<<"? 未知点数："<<unknownpoint<<
????????? "? 边长观测数："<<Num_Distance<<"? 角度观测数：" <<Num_Angle<<endl;
??? FileRead>>help_str;//跳过第2行描述语
??? FileRead>>Pre_Angle;
??? FileRead>>Pre_Distance_A;
??? FileRead>>Pre_Distance_B;
??? cout<<"边角网精度指标:"<<endl;
??? cout<<"角度观测精度（秒）:"<<Pre_Angle<<"? 固定误差A："<<Pre_Distance_A
?????? <<"? 比例误差B：" <<Pre_Distance_B<<endl;
??? //读入已知点的ID编号以及X坐标和Y坐标存入KnownPointXY
??? FileRead>>help_str;//跳过第3行描述语
??? KnownPointXY.resize(knownpoint,3);//设置矩阵KnownPointXY的大小
??? for(int i=0;i<knownpoint;i++)
??? {
??????? for(int j=0;j<3;j++)
??????? {
??????????? FileRead>>KnownPointXY(i,j);
??????? }
??? }
??? cout<<"已知点信息如下，请查阅："<<endl<<"已知点的ID编号 X Y"<<endl;
??? cout<<KnownPointXY<<endl;//输出矩阵KnownPointXY，供用户检查数据是否正确
??? //读入未知点的ID编号，存入KnownPointX0Y0
??? FileRead>>help_str;//跳过第4行描述语
??? UnKnownPointX0Y0.resize(unknownpoint,3);//设置矩阵UnKnownPointX0Y0的大小
??? for(int i=0;i<unknownpoint;i++)
??? {
??????? FileRead >> UnKnownPointX0Y0(i,0);
??????? for(int j=1;j<3;j++)//近似坐标X0和Y0初始为0
??????? {
??????????? UnKnownPointX0Y0(i,j)=0;
??????? }
??? }
??? cout<<"未知点信息如下，请查阅："<<endl<<"未知点的ID编号 X0 Y0(未计算近似坐标默认0)"<<endl;
??? cout<<UnKnownPointX0Y0<<endl;//检查输出矩阵UnKnownPointX0Y0，检查是否正确
??? //读入边长观测值的ID编号、边长/M、起点、终点到矩阵Distance_Obs
??? FileRead>>help_str;//跳过第5行描述语
??? Distance_Obs.resize(Num_Distance,4);//Num_Distance*4的矩阵
??? for(int i=0;i<Num_Distance;i++)
??? {
??????? for(int j=0;j<4;j++)
??????? {
??????????? FileRead >> Distance_Obs(i,j);
??????? }
??? }
??? cout<<"边长观测值信息如下，请查阅："<<endl;
??? cout<<Distance_Obs<<endl;//检查输出矩阵Distance_Obs，检查是否正确
??? //读入测站ID，后视照准ID_1，前视照准ID_2，夹角观测值,测站ID和ID_2对应的边长ID到Angle_Obs
??? FileRead>>help_str;//跳过第6行描述语
??? Angle_Obs.resize(Num_Angle,5);//Num_Angle*5的矩阵
??? for(int i=0;i<Num_Angle;i++)
??? {
??????? for(int j=0;j<5;j++)
??????? {
??????????? FileRead >> Angle_Obs(i,j);
??????? }
??? }
??? FileRead.close();//文件读取完毕，关闭Data.txt文件
??? cout<<"角度观测值信息如下，请查阅："<<endl<<"测站ID? 后视ID_1? 前视ID_2? 夹角值? 边长ID"<<endl;
??? cout<<Angle_Obs<<endl;//检查输出矩阵Angle_Obs，检查是否正确
}
2.角度格式转换类AngleChange
（1）anglechange.h
#ifndef ANGLECHANGE_H
#define ANGLECHANGE_H
#include "indata.h"
class AngleChange:public InData
{
public:
??? Eigen::MatrixXd Angle_Rad;//定义存放角度转换为弧度后的矩阵
??? void DoAngleChange();//弧度转换
};
#endif // ANGLECHANGE_H
（2）anglechange.cpp
#include "anglechange.h"
#define PI 3.14159265358979312//π值
void AngleChange::DoAngleChange()
{
??? Angle_Rad=Angle_Obs;//设置矩阵Angle_Rad替换Angle_Obs其他元素
??? for(int i=0;i<Num_Angle;i++)//将输入角度（度.分秒）转化为弧度
??? {
??????? double help_rad =Angle_Obs(i,3);
??????? double deg=floor(help_rad);
??????? double min1=((help_rad*100)-(deg*100));
??????? double min2=floor(min1);
??????? double min=min2/60;
??????? double second=(min1-min2)/36;
??????? help_rad=deg+min+second;
??????? help_rad=help_rad*PI/180;
??????? Angle_Rad(i,3)=help_rad;//将弧度值重新载入矩阵Angle_Rad
??? }
??? cout<<"角度观测值转化为弧度："<<endl<<Angle_Rad<<endl;
}
3.计算近似坐标类X0Y0
（1）x0y0.h
#ifndef X0Y0_H
#define X0Y0_H
#include "anglechange.h"
class X0Y0:public AngleChange
{
public:
??? Eigen::MatrixXd AllPointX0Y0;//存放已知点坐标和未知点计算出来后的近似坐标
??? void DoX0Y0();//计算未知点近似坐标
};
#endif // X0Y0_H
（2）x0y0.cpp
#include "x0y0.h"
#include <cmath>
#define PI 3.14159265358979312//π值
void X0Y0::DoX0Y0()
{
??? AllPointX0Y0.resize(knownpoint+unknownpoint,3);//设置矩阵AllPointX0Y0大小
??? for(int i=0;i<knownpoint;i++)//将矩阵KnownPointXY的数据放到矩阵AllPointX0Y0
??? {
??????? for(int j=0;j<3;j++)
??????? AllPointX0Y0(i,j)=KnownPointXY(i,j);
??? }
??? for(int i=0;i<unknownpoint;i++)//将矩阵UnKnownPointX0Y0的数据放到矩阵AllPointX0Y0
??? {
??????? for(int j=0;j<3;j++)
??????? AllPointX0Y0(i+knownpoint,j)=UnKnownPointX0Y0(i,j);
??? }
??? for(int i=0;i<unknownpoint;i++)//计算未知点坐标，并将坐标对应存放到矩阵AllPointX0Y0
??? {
??????? for(int j=0;j<Num_Angle;j++)//循环查找未知点计算所需元素
??????? {
??????????? if(AllPointX0Y0(i+knownpoint,1)!=0)//检索并判断当前的未知点的X坐标是否为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坐标
??????????? {
????????????????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);
??????????????? AllPointX0Y0(i+knownpoint,1)=X0;//之后的未知点计算需要用到之前的未知点坐标，
??????????????? AllPointX0Y0(i+knownpoint,2)=Y0;//所以将计算出来的近似坐标存放回矩阵AllPointX0Y0
??????????? }
??????? }
??? }
??? cout<<"未知点的近似坐标："<<endl<<AllPointX0Y0<<endl;
}
4.平差模块类Ajust
（1）ajust.h
#ifndef AJUST_H
#define AJUST_H
#include <cmath>
#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();//误差椭圆参数计算
};
#endif // AJUST_H
（2）ajust.cpp
#include "ajust.h"
#define p11 206264.80624709636//ρ″值
#define PI 3.14159265358979312//π值
void Ajust::DoForBl()//计算系数阵B和自由项矩阵l
{
??? B.setZero(Num_Angle+Num_Distance,unknownpoint*2);//设置系数阵B阵大小并初始化为0
??? l.setZero(Num_Angle+Num_Distance,1);//同上
??? for(int i=0;i<Num_Distance;i++)//求边长误差方程的系数B，自由项l
??? {
??????? int IDstart=Distance_Obs(i,2);//边长的起点ID
??????? int IDend=Distance_Obs(i,3);//边长的终点ID
??????? double dX0=AllPointX0Y0(IDend-1,1)-AllPointX0Y0(IDstart-1,1);
??????? double dY0=AllPointX0Y0(IDend-1,2)-AllPointX0Y0(IDstart-1,2);
??????? double S0=sqrt(pow(dX0,2)+pow(dY0,2));
??????? double li=Distance_Obs(i,1)-S0;//计算边长自由项li的值
??????? l(i,0)=li*1000;//将自由项li的值存入矩阵l
??????? if(IDstart>knownpoint)//计算未知点待求参数对应的系数
??????? {
??????????? double Bi_start_x=-dX0/S0;
??????????? double Bi_start_y=-dY0/S0;
??????????? B(i,(IDstart-knownpoint)*2-2)=Bi_start_x;//系数的值存入矩阵B
??????????? B(i,(IDstart-knownpoint)*2-1)=Bi_start_y;
??????? }
??????? if(IDend>knownpoint)
??????? {
??????????? double Bi_end_x=dX0/S0;
??????????? double Bi_end_y=dY0/S0;
??????????? B(i,(IDend-knownpoint)*2-2)=Bi_end_x;//系数的值存入矩阵B
??????????? B(i,(IDend-knownpoint)*2-1)=Bi_end_y;
??????? }
??? }
??? for(int j=0;j<Num_Angle;j++)//求角度误差方程的系数B，自由项l
??? {
??????? int IDstation=Angle_Rad(j,0);//角度观测的测站ID
??????? int IDlast=Angle_Rad(j,1);//角度观测的后视ID
??????? int IDnext=Angle_Rad(j,2);//角度观测的前视ID
??????? //测站点IDstation至前视照准点IDnext的dX,dY,距离平方dS^2
??????? double dXs_n=AllPointX0Y0(IDnext-1,1)-AllPointX0Y0(IDstation-1,1);
??????? double dYs_n=AllPointX0Y0(IDnext-1,2)-AllPointX0Y0(IDstation-1,2);
??????? double dS2s_n=pow(dXs_n,2)+pow(dYs_n,2);
??????? //测站点IDstation至后视照准点IDlast的dX,dY,距离平方dS^2
??????? double dXs_l=AllPointX0Y0(IDlast-1,1)-AllPointX0Y0(IDstation-1,1);
??????? double dYs_l=AllPointX0Y0(IDlast-1,2)-AllPointX0Y0(IDstation-1,2);
??????? double dS2s_l=pow(dXs_l,2)+pow(dYs_l,2);
??????? //计算角度误差方程的自由项l
??????? double a_sn=atan2(dYs_n,dXs_n);//测站至前视的近似方位角
??????? if(a_sn<0)//使全为0~2π范围
??????? {
??????????? a_sn=2*PI+a_sn;
??????? }
??????? double a_sl=atan2(dYs_l,dXs_l);//测站至后视的近似方位角
??????? if(a_sl<0)//使全为0~2π范围
??????? {
??????????? a_sl=2*PI+a_sl;
??????? }
??????? double dAngle=a_sn-a_sl;
??????? if(dAngle<0)//使全为0~2π范围
??????? {
??????????? dAngle=2*PI+dAngle;
??????? }
??????? double lj=p11*(Angle_Rad(j,3)-dAngle);//计算角度误差自由项lj
??????? l(j+Num_Distance,0)=lj;//将自由项lj的值存入矩阵l，注意单位化为秒
??????? if(IDstation>knownpoint)//计算未知点待求参数对应的系数
??????? {
??????????? double Bj_station_x=(dYs_n/dS2s_n - dYs_l/dS2s_l)*p11;
??????????? double Bj_station_y=(dXs_l/dS2s_l - dXs_n/dS2s_n)*p11;
??????????? B(j+Num_Distance,(IDstation-knownpoint)*2-2)=Bj_station_x/1000;//系数的值存入矩阵B
??????????? B(j+Num_Distance,(IDstation-knownpoint)*2-1)=Bj_station_y/1000;
??????? }
??????? if(IDlast>knownpoint)
??????? {
??????????? double Bj_last_x= (dYs_l/dS2s_l)*p11;
??????????? double Bj_last_y= (-dXs_l/dS2s_l)*p11;
??????????? B(j+Num_Distance,(IDlast-knownpoint)*2-2)=Bj_last_x/1000 ;//系数的值存入矩阵B
??????????? B(j+Num_Distance,(IDlast-knownpoint)*2-1)=Bj_last_y/1000;
??????? }
??????? if(IDnext>knownpoint)
??????? {
??????????? double Bj_next_x=(dYs_n/dS2s_n)*(-p11);
??????????? double Bj_next_y= (dXs_n/dS2s_n)*p11;
??????????? B(j+Num_Distance,(IDnext-knownpoint)*2-2)=Bj_next_x/1000;//系数的值存入矩阵B
??????????? B(j+Num_Distance,(IDnext-knownpoint)*2-1)=Bj_next_y/1000;
??????? }
??? }
}
void Ajust::DoForP()//定权矩阵P
{
??? double P_distance;//观测边的权
??? double P_angle;//观测角的权
??? P.setZero(Num_Distance+Num_Angle,Num_Distance+Num_Angle);//设置权阵P大小并初始化为0
??? 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;
??? }
}
void Ajust::DoForVx()//计算改正数矩阵V、待定参数矩阵x
{
??? V.setZero(Num_Angle+Num_Distance,1);//设置改正数矩阵大小并初始化为0
??? x.setZero(unknownpoint*2,1);//设置待求参数阵x大小并初始化为0
??? MatrixXd BT=B.transpose();//系数阵B的转置
??? Nbb=BT*P*B;//法方程系数阵
??? MatrixXd W=BT*P*l;//法方程常数项
??? x=Nbb.inverse()*W;//待求参数阵x=Nbb的逆阵*法方程常数项W
??? V=B*x-l;//计算改正数矩阵
}
void Ajust::Accuracy_assess()//进行精度评定
{
??? MatrixXd VtPV=V.transpose()*P*V;//VT*P*V
??? double r=Num_Distance+Num_Angle-unknownpoint*2;//自由度，即多余观测数
??? double sigma0_2=VtPV(0,0)/r;//单位权方差
??? sigma0=sqrt(sigma0_2);//单位权中误差
??? QXX=Nbb.inverse();//平差参数的协因数阵
??? DXX=sigma0_2*QXX;//平差参数的协方差阵
}
void Ajust::Error_ellipse()//计算误差椭圆参数
{
??? Fhi0_E.setZero(unknownpoint,1);//设置矩阵Fhi0_E大小并初始化为0
??? E.setZero(unknownpoint,1);
??? Fhi0_F.setZero(unknownpoint,1);
??? F.setZero(unknownpoint,1);//同上
??? double _2Fhi0,_2Fhi0_180;//极值方向的两个根的两倍
??? for(int i=0;i<QXX.rows();i=i+2)
??? {
??????? double Qxx=QXX(i,i);
??????? double Qxy=QXX(i,i+1);
??????? double Qyy=QXX(i+1,i+1);
??????? _2Fhi0=atan2((2*Qxy),(Qxx-Qyy));//计算极值根的两倍
??????? if(_2Fhi0>=0)//若极值根>=0,另一个根+π
??????? {
??????????? _2Fhi0_180=_2Fhi0+PI;
??????? }
??????? if(_2Fhi0<0)//若极值根<0,另一个根+π
??????? {
??????????? _2Fhi0_180=_2Fhi0+PI;
??????????? _2Fhi0=_2Fhi0+2*PI;//若_2Fhi0根<0，将化为0~2π范围表示
??????? }
??????? double Fhi0=_2Fhi0/2;//极大值或极小值的方向
??????? double Fhi0_90=_2Fhi0_180/2;//极大值或极小值的方向
??????? double E_2,F_2;//极大值或极小值的平方
??????? E_2=pow(sigma0,2)*(Qxx*pow(cos(Fhi0),2)+Qyy*pow(sin(Fhi0),2)+Qxy*sin(2*Fhi0));
??????? F_2=pow(sigma0,2)*(Qxx*pow(cos(Fhi0_90),2)+Qyy*pow(sin(Fhi0_90),2)+Qxy*sin(2*Fhi0_90));
??????? //判断极大值与极小值
??????? if(E_2>=F_2)//若E>F，Fhi0这个根就是极大值E的方向，若E=F，误差椭圆即为圆，这里归类到这一类
??????? {
??????????? Fhi0_E(i/2,0)=Fhi0;
??????????? E(i/2,0)=sqrt(E_2);
??????????? Fhi0_F(i/2,0)=Fhi0_90;
??????????? F(i/2,0)=sqrt(F_2);
??????? }
??????? if(E_2<F_2)//若E<F，Fhi0这个根就是极小值F的方向
??????? {
??????????? Fhi0_E(i/2,0)=Fhi0_90;
??????????? E(i/2,0)=sqrt(F_2);
??????????? Fhi0_F(i/2,0)=Fhi0;
??????????? F(i/2,0)=sqrt(E_2);
??????? }
??? }
}
5.数据输出PutData.txt类OutData
（1）outdata.h
#ifndef OUTDATA_H
#define OUTDATA_H
#include "ajust.h"
class OutData:public Ajust
{
public:
??? void DoPutData();//将数据输出至PutData.txt文本文档
};
#endif // OUTDATA_H
（2）outdata.cpp
#include "outdata.h"
#define PI 3.14159265358979312//π值
#include <iomanip>
//I/O流控制头文件,用于控制输出数据的有效数字位数、指定小数位数、域宽等格式
void OutData::DoPutData()
{
??? ofstream FileWrite;//输出文件流对象
??? FileWrite.open("PutData.txt",ios::out|ios::trunc);//打开项目下的PutData.txt文件
??? if(!FileWrite)
??? {
??????? cout<<"输出文件打开时出错，请检查路径"<<endl;
??????? exit(1);
??? }
??? FileWrite<<"---------------------平差成果综述如下--------------------"<<'\n';
??? FileWrite<<"1.平差值计算："<<'\n';
??? FileWrite<<"(1)坐标平差值："<<'\n';
??? //设置标头及合适的宽度
??? FileWrite<<setfill('\0')<<setw(10)<<"点ID";
??? FileWrite<<setfill('\0')<<setw(20)<<"近似X0/m";
??? FileWrite<<setfill('\0')<<setw(15)<<"Vx/mm";
??? FileWrite<<setfill('\0')<<setw(20)<<"平差X/m";
??? FileWrite<<setfill('\0')<<setw(20)<<"近似Y0/m";
??? FileWrite<<setfill('\0')<<setw(15)<<"Vy/mm";
??? FileWrite<<setfill('\0')<<setw(20)<<"平差Y/m"<<'\n';
??? for(int i=1;i<=unknownpoint;i++)
??? {
??????? FileWrite<<setfill('\0')<<setw(10)<<fixed<<setprecision(0)<<AllPointX0Y0(i+knownpoint-1,0);//未知点ID
??????? FileWrite<<setfill('\0')<<setw(20)<<fixed<<setprecision(4)<<AllPointX0Y0(i+knownpoint-1,1);//近似坐标X0
??????? FileWrite<<setfill('\0')<<setw(15)<<fixed<<setprecision(1)<<x(2*i-2,0);//Vx
??????? FileWrite<<setfill('\0')<<setw(20)<<fixed<<setprecision(4)<<AllPointX0Y0(i+knownpoint-1,1)+0.001*x(2*i-2,0);//平差值X
??????? FileWrite<<setfill('\0')<<setw(20)<<fixed<<setprecision(4)<<AllPointX0Y0(i+knownpoint-1,2);//近似坐标Y0
??????? FileWrite<<setfill('\0')<<setw(15)<<fixed<<setprecision(1)<<x(2*i-1,0);//Vy
??????? FileWrite<<setfill('\0')<<setw(20)<<fixed<<setprecision(4)<<AllPointX0Y0(i+knownpoint-1,2)+0.001*x(2*i-1,0)<<'\n';//平差值Y
??? }
??? FileWrite<<'\n';
??? FileWrite<<"(2)边长平差值："<<'\n';
??? //设置标头及合适的宽度
??? FileWrite<<setfill('\0')<<setw(10)<<"边ID";
??? FileWrite<<setfill('\0')<<setw(20)<<"观测值/m";
??? FileWrite<<setfill('\0')<<setw(15)<<"改正数/mm";
??? FileWrite<<setfill('\0')<<setw(20)<<"平差值/m"<<'\n';
??? for(int i=0;i<Num_Distance;i++)
??? {
??????? FileWrite<<setfill('\0')<<setw(10)<<fixed<<setprecision(0)<<Distance_Obs(i,0);//边长ID
??????? FileWrite<<setfill('\0')<<setw(20)<<fixed<<setprecision(4)<<Distance_Obs(i,1);//边长观测值
??????? FileWrite<<setfill('\0')<<setw(15)<<fixed<<setprecision(1)<<V(i,0);//改正数/mm
??????? FileWrite<<setfill('\0')<<setw(20)<<fixed<<setprecision(4)<<Distance_Obs(i,1)+0.001*V(i,0)<<'\n';//边长平差值
??? }
??? FileWrite<<'\n';
??? FileWrite<<"(3)角度平差值："<<'\n';
??? //设置标头及合适的宽度
??? FileWrite<<setfill('\0')<<setw(10)<<"角ID";
??? FileWrite<<setfill('\0')<<setw(25)<<"观测值(度.分秒)";
??? FileWrite<<setfill('\0')<<setw(18)<<"改正数(″)";
??? FileWrite<<setfill('\0')<<setw(25)<<"平差值(度.分秒)"<<'\n';
??? for(int i=0;i<Num_Angle;i++)
??? {
??????? FileWrite<<setfill('\0')<<setw(10)<<fixed<<setprecision(0)<<i+1;//角ID
??????? FileWrite<<setfill('\0')<<setw(20)<<fixed<<setprecision(5)<<Angle_Obs(i,3);//观测值(度.分秒)
??????? FileWrite<<setfill('\0')<<setw(15)<<fixed<<setprecision(2)<<V(i+Num_Distance,0);//改正数(″)
??????? FileWrite<<setfill('\0')<<setw(20)<<fixed<<setprecision(5）<<Angle_Obs(i,3)+0.0001*V(i+Num_Distance,0)<<'\n';//平差值
??? }
??? FileWrite<<'\n';
??? FileWrite<<"2.精度评定："<<'\n';
??? FileWrite<<"(1)单位权中误差："<<fixed<<setprecision(1)<<sigma0<<"″"<<'\n'<<'\n';
??? FileWrite<<"(2)坐标平差值的协因数阵QXX:"<<'\n';
??? FileWrite<<fixed<<setprecision(3)<<QXX<<'\n'<<'\n';
??? FileWrite<<"(3)坐标平差值的X、Y方向中误差，点位中误差："<<'\n';
??? //设置标头及合适的宽度
??? FileWrite<<setfill('\0')<<setw(10)<<"ID";
??? FileWrite<<setfill('\0')<<setw(10)<<"Mx/mm";
??? FileWrite<<setfill('\0')<<setw(10)<<"My/mm";
??? FileWrite<<setfill('\0')<<setw(10)<<"M/mm"<<'\n';
??? for(int i=0;i<unknownpoint;i++)
??? {
??????? double Mx2=DXX(2*i,2*i);//坐标平差值的X方向中误差
??????? double My2=DXX(2*i+1,2*i+1);//坐标平差值的Y方向中误差
??????? FileWrite<<setfill('\0')<<setw(10)<<fixed<<setprecision(0)<<UnKnownPointX0Y0(i,0);
??????? FileWrite<<setfill('\0')<<setw(15)<<fixed<<setprecision(1)<<sqrt(Mx2);
??????? FileWrite<<setfill('\0')<<setw(15)<<fixed<<setprecision(1)<<sqrt(My2);
??????? FileWrite<<setfill('\0')<<setw(15)<<fixed<<setprecision(1)<<sqrt(Mx2+My2)<<'\n';
??? }
??? FileWrite<<'\n';
??? FileWrite<<"3.误差椭圆参数："<<'\n';
??? FileWrite<<setfill('\0')<<setw(10)<<"ID";
??? FileWrite<<setfill('\0')<<setw(10)<<"E/mm";
??? FileWrite<<setfill('\0')<<setw(20)<<"φE(度.分秒)";
??? FileWrite<<setfill('\0')<<setw(10)<<"F/mm";
??? FileWrite<<setfill('\0')<<setw(20)<<"φF(度.分秒)"<<'\n';
??? for(int i=0;i<unknownpoint;i++)
??? {
??????? FileWrite<<setfill('\0')<<setw(10)<<fixed<<setprecision(0)<<UnKnownPointX0Y0(i,0);
??????? FileWrite<<setfill('\0')<<setw(10)<<fixed<<setprecision(1)<<E(i,0);
??????? double fhi_E=Fhi0_E(i,0);//将极大值E方向由弧度转为度.分秒
??????? double deg=fhi_E*180/PI;//度数
??????? double min=(deg-floor(deg))*60;//分数
??????? double second=(min-floor(min))*60;//秒数
??????? double deg_min_second=floor(deg)+0.01*floor(min)+0.0001*floor(second);//度.分秒
??????? FileWrite<<setfill('\0')<<setw(20)<<fixed<<setprecision(4)<<deg_min_second;
??????? FileWrite<<setfill('\0')<<setw(10)<<fixed<<setprecision(1)<<F(i,0);
??????? double fhi_F=Fhi0_F(i,0);//将极小值F方向由弧度转为度.分秒
??????? double _deg=fhi_F*180/PI;//度数
??????? double _min=(_deg-floor(_deg))*60;//分数
??????? double _second=(_min-floor(_min))*60;//秒数
??????? double _deg_min_second=floor(_deg)+0.01*floor(_min)+0.0001*floor(_second);//度.分秒
??????? FileWrite<<setfill('\0')<<setw(20)<<fixed<<setprecision(4)<<_deg_min_second<<'\n';
??? }
??? FileWrite.close();//关闭PutData.txt文件
}