水准网程序设计

误差理论与测量平差基础
（MATLAB）
实习报告
学号：
姓名：
班级：
专业：测绘工程
课程名称：误差理论与测量平差基础
任课老师：
20XX 年5 月
一、任务概述
利用MATLAB或者C++编程间接平差程序，通过该程序读取观测数据文件，并计算出平差结果。

C++
// 水准网平差.cpp : 定义控制台应用程序的入口点。

//
#include "stdafx.h"
#include <cstdio>
#include "math.h"
#include <iostream>
#include "string.h"
using namespace std;
bool inverse(double a[],int n)
{
double *a0=new double[n];
for(int k=0;k<n;k++)
{
double a00=a[0];
if(a00+1.0==1.0)
{
delete []a0;
return false;
}
for(int i=1;i<n;i++)
{
double ai0 = a[i*(i+1)/2];
if(i<=n-k-1)a0[i]= -ai0/a00;
else a0[i]= ai0/a00;
for(int j=1;j<=i;j++)
{
a[(i-1)*i/2+j-1]=a[i*(i+1)/2+j]+ai0*a0[j];
}
}
for(int i=1;i<n;i++)
{
a[(n-1)*n/2+i-1]=a0[i];
}
a[n*(n+1)/2-1]=1.0/a00;
}
delete []a0;
return true;
}
int matrix_s(int i,int j)
{
return (i>=j)? i*(i+1)/2+j :j*(j+1)/2+i;
}
class CLevelingAdjust
{
public:
CLevelingAdjust();
virtual ~CLevelingAdjust();
int m_Lnumber; //高差总数
int m_Pnumber; //总点数
int m_knPnumber; //已知点数
int unPnumber; //未知点个数
double m_pvv; //[pvv]
FILE *resultfp; //文件指针,输出计算结果
double m_Sigma; //验前单位权中误差（粗差探测、闭合差计算用）int *StartP; //高差起点号
int *EndP; //高差终点号
char **Pname; //点名地址数组
double *L; //观测值数组
double *Height; //高程值数组
double *P; //观测值的权
double *ATPA,*ATPL; //法方程系数矩阵与自由项
double *dX; //参数平差值（高程改正数）
double *V; //残差数组
double m_mu; //验后单位权中误差
void Inputdata(char *datafile);//输入原始数据函数
void Printdata(); //输出原始数据函数
int GetStationNumber(char *name); //获取点号函数
void ca_H0(); //近似高程计算函数
void ca_ATPA(); //法方程组成函数
void ca_dX(); //平差值计算函数
void PrintResult(); //精度估计与平差值输出函数
double ca_V(); //残差计算函数
};
CLevelingAdjust::CLevelingAdjust()
{
m_Lnumber=0;
m_Pnumber=0;
}
CLevelingAdjust::~CLevelingAdjust()
{
if(m_Lnumber>0)
{
delete []StartP;
delete []EndP;
delete []L;
delete []P;
delete []V;
}
if(m_Pnumber>0)
{
delete []Height;
delete []ATPA;
delete []ATPL;
delete []dX;
for(int i=0; i<m_Pnumber;i++)
if(Pname[i]!=NULL)delete[](Pname[i]);
delete []Pname;
}
}
// 原始数据输入函数
void CLevelingAdjust::Inputdata(char *datafile)
{
FILE *fr;
if((fr=fopen(datafile,"r"))==NULL)
{
cout<<"\n 数据文件打不开!"<<endl;
exit(0);
}
fscanf(fr,"%d%d%d",&m_Lnumber,&m_Pnumber,&m_knPnumber); unPnumber=m_Pnumber-m_knPnumber;
Height=new double [m_Pnumber];
dX=new double [m_Pnumber];
ATPA=new double [m_Pnumber*(m_Pnumber+1)/2];
ATPL=new double [m_Pnumber];
StartP=new int [m_Lnumber];
EndP=new int [m_Lnumber];
L=new double [m_Lnumber];
V=new double [m_Lnumber];
P=new double [m_Lnumber];
//fscanf(fp,"%lf",&m_Sigma);
Pname=new char* [m_Pnumber];
for(int i=0;i<m_Pnumber;i++)
{
// GetStationNumber函数根据Pname[i]是否为NULL
// 确定Pname[i]是否为点名地址
Pname[i] = NULL;
}
char buffer[100];
// 读取已知高程数据
for(int i=0;i<=m_knPnumber-1;i++)
{
fscanf(fr,"%s",buffer);
int c=GetStationNumber(buffer);
fscanf(fr,"%lf",&Height[c]);
}
// 读取观测数据
for(int i=0;i<m_Lnumber;i++)
{
fscanf(fr,"%s",buffer); //读取高程起点名
StartP[i]=GetStationNumber(buffer);
if(StartP[i]<0)
{
fprintf(resultfp,"\n数据文件出错：");
fprintf(resultfp,"\n第%d个高差的起始点名为\"%s\"",i+1,buffer);
fclose(resultfp);
exit(0);
}
fscanf(fr,"%s",buffer);//读取高程终点
EndP[i]=GetStationNumber(buffer);
if(EndP[i]<0)
{
fprintf(resultfp,"\n数据文件出错：");
fprintf(resultfp,"\n第%d个高差终点的点名为\"%s\"",i+1,buffer);
fclose(resultfp);
exit(0);
}
fscanf(fr,"%lf%lf",&L[i],&P[i]); //读取高差值与路线长度
P[i]=1.0/P[i];
}
fclose(fr);
}
// 原始数据写到结果文件
void CLevelingAdjust::Printdata()
{
int i;
fprintf(resultfp,"\n 观测值总数: %d 总点数: %d 已知点数：%d \n",
m_Lnumber, m_Pnumber,m_knPnumber);
//fprintf(resultfp,"\n 验前单位权中误差：%lf" ,m_Sigma);
fprintf(resultfp,"\n\n 已知高程:\n");
for(i=0;i<=m_knPnumber-1;i++)
{
fprintf(resultfp,"\n%5d %8s ",i+1,Pname[i]);
fprintf(resultfp,"%10.4lf ",Height[i]);
}
fprintf(resultfp,"\n\n 高差观测值:\n");
for(i=0;i<=m_Lnumber-1;i++)
{
fprintf(resultfp,"\n%5d %8s %8s",i+1,Pname[StartP[i]],Pname[EndP[i]]);
fprintf(resultfp,"%12.4lf %10.3lf",L[i],1.0/P[i]);
}
}
// 点名存贮，返回点名对应的点号
int CLevelingAdjust::GetStationNumber(char *name)
{
for(int i=0; i<m_Pnumber; i++)
{
if(Pname[i]!=NULL)
{
//将待查点名与已经存入点名数组的点名比较
if(strcmp(name,Pname[i])==0)return i;
}
else
{
//待查点名是一个新的点名，将新点名的地址放到Pname数组中
int len = strlen(name);
Pname[i] = new char[len+1];
strcpy(Pname[i], name);
return i;
}
}
return -1; //Pname数组已经存满，且没有待查点名
}
//高程近似值计算
void CLevelingAdjust::ca_H0()
{
for(int i=m_knPnumber;i<m_Pnumber;i++)Height[i]=-9999.9;
int jj=0; //计算出近似高程的点数
for(int ii=1;;ii++)
{
for(int i=0;i<m_Lnumber;i++)
{
int k1=StartP[i]; //高差起点号
int k2=EndP[i]; //高差终点号
if(Height[k1]>-9999.0 && Height[k2]<-9999.0)
{
Height[k2]=Height[k1]+L[i];
jj++;
}
else if(Height[k1]<-9999.0 && Height[k2]>-9999.0)
{
Height[k1]=Height[k2]-L[i];
jj++;
}
}
if(jj==(m_Pnumber-m_knPnumber))break;
if(ii>(m_Pnumber-m_knPnumber))
{
fprintf(resultfp,"\n\n下列点无法计算出概略高程:\n");
for(int i=0;i<m_Pnumber;i++)
{
if(Height[i]<-9999.0) printf("\n%s",Pname[i]);
}
cout<<"近似高程计算失败！"<<endl;
fclose(resultfp);
exit(0);
}
}
}
void CLevelingAdjust::ca_ATPA()
{
int t=m_Pnumber;
for(int i=0; i<t*(t+1)/2; i++) ATPA[i]=0.0;
for(int i=0; i<t; i++) ATPL[i]=0.0;
for(int i=0;i<m_knPnumber;i++)
{
ATPA[matrix_s(i,i)]=1.0E08;
}
for(int k=0; k<m_Lnumber; k++)
{
int i=StartP[k];
int j=EndP[k];
double Pk=P[k];
double Lk=L[k]-(Height[j]-Height[i]);
ATPL[i]-=Pk*Lk;
ATPL[j]+=Pk*Lk;
ATPA[matrix_s(i,i)]+=Pk;
ATPA[matrix_s(j,j)]+=Pk;
ATPA[matrix_s(i,j)]-=Pk;
}
}
void CLevelingAdjust::ca_dX()
{
if(!inverse(ATPA,m_Pnumber))
{
cout<<"法方程系数矩阵降秩!"<<endl;
fclose(resultfp);
exit(0);
}
for(int i=0; i<m_Pnumber; i++)
{
double xi=0.0;
for(int j=0; j<m_Pnumber; j++)
{
xi+=ATPA[matrix_s(i,j)]*ATPL[j];
}
dX[i]=xi;
Height[i]+=xi;
}
}
double CLevelingAdjust::ca_V()
{
double pvv=0.0;
for(int i=0;i<=m_Lnumber-1;i++)
{
int k1=StartP[i];
int k2=EndP[i];
V[i]=Height[k2]-Height[k1]-L[i];
pvv+=V[i]*V[i]*P[i];
}
m_mu=sqrt(pvv/(unPnumber));
return(m_mu);
}
void CLevelingAdjust::PrintResult()
{
fprintf(resultfp,"\n\n ==== 高程平差值及其精度 ====\n");
fprintf(resultfp,"\n 点名近似高程改正数高程平差值中误差\n");
for(int i=0; i<m_Pnumber; i++)
{
fprintf(resultfp,"\n %5s ",Pname[i]);
double dx=dX[i];
double qii=ATPA[matrix_s(i,i)];
fprintf(resultfp,"%12.4lf%9.4lf%12.4lf%9.4lf",
Height[i]-dx,dx,Height[i],sqrt(qii)*m_mu);
}
}
int _tmain(int argc, _TCHAR* argv[])
{
char *datafile ="算例\\data.txt"; char *resultfile ="算例\\result.txt";
printf("\n\n水准网平差\n");
printf("数据文件： %s\n",datafile); printf("结果文件： %s\n",resultfile);
FILE *fp=fopen(resultfile,"w");
if (fp==NULL)
{
cout<<"创建结果文件失败！"<<endl;
return -1;
}
CLevelingAdjust level;
level.resultfp=fp;
level.Inputdata(datafile);
level.Printdata();
level.ca_H0();
level.ca_ATPA();
level.ca_dX();
level.ca_V();
level.PrintResult();
fclose(fp);
cout<<"解算成功！"<<endl;
char jj;
cin>>jj;
return 0;
MATLAB
图一图二
图三图四三、水准网图
四、输入的数据格式
数据格式为TXT文件，如图所示：
TXT文件格式说明：
（1）第一行格式
第一行分别表示观测个数5个，水准点数4个，未知点3个，
已知点1个，所有数据用英文逗号隔开
（2）已知点数据格式
第二行开始是已知点点号和高程，一行列一个已知点点号和高
程，由于该水准网只有一个已知点，所有只能列出一行。

图中
表示已知点点号为1，高程为237.483m
（3）测站起始点号格式
（4）测站终点点号格式（5）高差格式
（6）距离格式该部分表示测站的起始点点号
该部分表示测站的终点点号该部分表示各测站的高差
该部分表示各测站的距离S 五、流程图
六、附件代码
function SDJianJiePingCha()
[FileName,PathName] = uigetfile('*.txt','打开水准观测数据');%打开文件
f=csvread( strcat(PathName,FileName));%打开文件并存在矩阵f中
point=f(1,2);%获取所有水准点个数
n=f(1,1);%获得观测个数n
t=f(1,3);%获得必要观测个数t
y=f(1,4);%获得已知点个数y
XX=zeros(point,1);%初始化XX阵等于0，方便下面把已知点高程和未知点参数估值放到XX阵
B=zeros(n,t);%初始化B阵，方便下面求V=Bx-l中的系数阵B；
for j=1:y
XX(j,1)=f(j+1,2);%把已知点高程放到XX阵中
end
data=f((2+y):end,:);%从文件中获取观测数据，并放到data阵中
h=data(:,3);%从data中获取观测高差，并放到h阵中
P=zeros(n);%初始化权阵P
for j=1:n
P(j,j)=10/data(j,4);%以10km观测值为单位权误差计算权阵P
end
for i=1:n%通过循环求B阵
point1=data(i,1);%获取某个测站的起始点号
point2=data(i,2);%获取某个测站的终点点号
if point1>y&&point2>y%当某测站起始点和终点高程都未知时，求B阵第i行B(i,point1-y)=-1;
B(i,point2-y)=1;
elseif point1<=y&&point2>y%当起始点高程已知和终点高程未知时，求B阵第i 行
B(i,point2-y)=1;
XX(point2,1)=XX(point1,1)+h(i,1);%求第i个参数估值
elseif point1>y&&point2<=y%当起始点高程未知和终点高程已知时，求B阵第i 行
B(i,point1-y)=-1;
XX(point1,1)=XX(point2,1)-h(i,1);%求第i个参数估值
end
end
l=zeros(n,1);%初始化小l阵，方便下面求V=Bx-l中的系数阵l；
for i=1:n%通过循环求小l
point1=data(i,1);
point2=data(i,2);
l(i,1)=-(XX(point2,1)-XX(point1,1)-h(i,1));
end
%带入间接平差数学模型公式进行计算：
r=n-t;%求多余观测数N=B'*P*B;
W=B'*P*l;
x=N\W;
X=XX((y+1):end,1)+x; V=B*x-l;
L=h+V;
a0=sqrt(V'*P*V/r); Qxx=inv(N);
Dxx=a0*a0*inv(N);
%输出计算结果：disp('参数改正数：') x=x'
disp('参数平差值：') X=X'
disp('观测值改正数：') V=V'
disp('观测值平差值：') L=L'
disp('协方差阵：') Dxx
disp('单位权方差：') a0
disp('协因数阵：') Qxx
B
l
P
N
W
end。