水准网间接平差程序设计(C++)

////////////////////////////////////////////////////// visual C++6.0 编译通过 /////////////////////////////////////////////////////////////////////////////////////////////////////////// 参考资料 //// 部分网络资料 //// 宋力杰《测量平差程序设计》 ////连壁《基于matlab的控制网平差程序设计》 /////////////////////////////////////////////////////#include<iostream>#include<fstream>#include <stdlib.h>#include<math.h>#include <iomanip>using namespace std;//////////////////////////////////////////////////////////////////////////class class SZWPC{private:int gcz_zs; //高差总数int szd_zs; //总点数int yz_szd_zs; //已知点数double m_pvv; //[pvv]int *qsd_dh; //高差起点号int *zd_dh; //高差终点号char **dm; //点名地址数组double *gcz; //观测值数组double *szd_gc; //高程值数组double *P; //观测值的权double *ATPA,*ATPL; //法方程系数矩阵与自由项double *dX; //高程改正数、平差值double *V; //残差double m_mu; //单位权中误差public:SZWPC();~SZWPC();int ij(int i,int j);//对称矩阵下标计算函数bool inverse(double a[],int n);//对称正定矩阵求逆(仅存下三角元素)(参考他人）void inputdata(char *datafile);//输入原始数据函数int dm_dh(char *name); //点名转点号void ca_H0(); //近似高程计算函数void ca_ATPA(); //法方程组成函数void ca_dX(); //高程平差值计算函数void printresult(char *resultfile); //精度估计与平差值输出函数double ca_V(); //残差计算函数void zxecpc(char *resultfile);//最小二乘平差函数};////////////////////////////////////////////////////////////////////// // 构造函数SZWPC::SZWPC(){gcz_zs=0;szd_zs=0;yz_szd_zs=0;}////////////////////////////////////////////////////////////////////// // 析构函数SZWPC::~SZWPC(){if(gcz_zs>0){delete []qsd_dh;delete []zd_dh;delete []gcz;delete []P;delete []V;}if(szd_zs>0){delete []szd_gc;delete []ATPA;delete []ATPL;delete []dX;for(int i=0; i<szd_zs;i++)if(dm[i]!=NULL)delete[](dm[i]);delete []dm;}}////////////////////////////////////////////////////////////////////////// // 对称矩阵下标计算函数int SZWPC::ij(int i,int j){return (i>=j)? i*(i+1)/2+j :j*(j+1)/2+i;}////////////////////////////////////////////////////////////////////////// // 对称正定矩阵求逆(仅存下三角元素)(参考他人）bool SZWPC::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(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;}/////////////////////////////////////////////////////////////////////// 原始数据输入函数void SZWPC::inputdata(char *datafile){ifstream infile(datafile,ios::in);if(! infile){cerr<<" Open error!"<<endl;}infile>>gcz_zs>>szd_zs>>yz_szd_zs;int unPnumber=szd_zs-yz_szd_zs;szd_gc=new double [szd_zs];dX=new double [szd_zs];ATPA=new double [szd_zs*(szd_zs+1)/2];ATPL=new double [szd_zs];qsd_dh=new int [gcz_zs];zd_dh=new int [gcz_zs];gcz=new double [gcz_zs];V=new double [gcz_zs];P=new double [gcz_zs];dm=new char* [szd_zs];for(int i=0;i<szd_zs;i++){dm[i] = NULL;// dm_dh函数根据dm[i]是否为NULL确定dm[i]是否为点名地址}char buffer[128]; //临时数组，保存从文件中读到的点名for( i=0;i<=yz_szd_zs-1;i++)// 读取已知高程数据{infile>>buffer;int c=dm_dh(buffer);infile>>szd_gc[i];}for(i=0;i<gcz_zs;i++)// 读取观测数据{infile>>buffer; //读取高程起点名qsd_dh[i]=dm_dh(buffer);infile>>buffer;//读取高程终点zd_dh[i]=dm_dh(buffer);infile>>gcz[i]>>P[i]; //读取高差值与路线长度P[i]=1.0/P[i];//线路长转化为观测值的权}infile.close();}//////////////////////////////////////////////////////////////////// 点名转点号，返回点名对应的点号int SZWPC::dm_dh(char *name){for(int i=0; i<szd_zs; i++){if(dm[i]!=NULL){if(strcmp(name,dm[i])==0)return i;//将待查点名与已经存入点名数组的点名比较，若存在返回点号}else{int len = strlen(name);//判断点名长度dm[i] = new char[len+1];//为点名申请存储空间strcpy(dm[i], name);//待查点名是一个新的点名，将新点名的地址放到dm 数组中return i;//返回点号}}return -1; //dm数组已经存满，且没有待查点名}////////////////////////////////////////////////////////////////////////////高程近似值计算void SZWPC::ca_H0(){for(int i=yz_szd_zs;i<szd_zs;i++)szd_gc[i]=-10000.9;//为计算机设置辨别未知高程点的标志for(int j=1;;j++){int k=0; //计算出近似高程的点数for(i=0;i<gcz_zs;i++){int k1=qsd_dh[i]; //高差起点号int k2=zd_dh[i]; //高差终点号if(szd_gc[k1]>-10000.0 && szd_gc[k2]<-10000.0)//k1点高程或高程近似值已知，k2点高程或高程近似值未知{szd_gc[k2]=szd_gc[k1]+gcz[i];//计算近似高程k++;}elseif(szd_gc[k1]<-10000.0 && szd_gc[k2]>-10000.0)//k2点高程或高程近似值已知，k1点高程或高程近似值未知{szd_gc[k1]=szd_gc[k2]-gcz[i];//计算近似高程k++;}}if(k==(szd_zs-yz_szd_zs))break;//所有的近似高程计算完成，退出}}//////////////////////////////////////////////////////////////////////////// 组成法方程void SZWPC::ca_ATPA(){//int t=szd_zs;for(int i=0; i<szd_zs*(szd_zs+1)/2; i++) ATPA[i]=0.0;//赋初值for(i=0; i<szd_zs; i++) ATPL[i]=0.0;//赋初值for(int k=0; k<gcz_zs; k++){int i=qsd_dh[k];//获取点号int j=zd_dh[k];//获取点号double Pk=P[k];//获取权值double lk=gcz[k]-(szd_gc[j]-szd_gc[i]);//获得第k个自由项ATPL[i]-=Pk*lk;//获得法方程自由项ATPL[j]+=Pk*lk;ATPA[ij(i,i)]+=Pk;//获得法方程系数矩阵ATPA[ij(j,j)]+=Pk;ATPA[ij(i,j)]-=Pk;}}////////////////////////////////////////////////////////////////////////// // 高程平差值计算void SZWPC::ca_dX(){for(int i=0;i<yz_szd_zs;i++) ATPA[ij(i,i)]=1.0e30;//处理已知点if(!inverse(ATPA,szd_zs))//矩阵求逆{cerr<<"法方程系数矩阵降秩!"<<endl;//矩阵为奇异矩阵，无法求逆exit(0);//退出程序}for(i=0; i<szd_zs; i++)//计算高程改正数{double xi=0.0;for(int j=0; j<szd_zs; j++){xi+=ATPA[ij(i,j)]*ATPL[j];}dX[i]=xi;szd_gc[i]+=xi;//计算高程平差值}}////////////////////////////////////////////////////////////////////////// // 残差计算double SZWPC::ca_V(){double pvv=0.0;for(int i=0;i<=gcz_zs-1;i++){int k1=qsd_dh[i];int k2=zd_dh[i];V[i]=szd_gc[k2]-szd_gc[k1]-gcz[i];pvv+=V[i]*V[i]*P[i];}return(pvv);}//////////////////////////////////////////////////////////////////////////// 原始数据和平差值输出void SZWPC::printresult(char *resultfile){double pvv=ca_V(); // 残差计算ofstream outfile(resultfile,ios::out);//以输出方式打开文件，若文件不存在，创建文件//输出原始观测数据outfile<<endl<<"观测总数："<<gcz_zs<<" "<<"总点数:"<<szd_zs;outfile<<" "<<"已知点数："<<yz_szd_zs<<endl;outfile<<endl<<"===================== 已知高程====================="<<endl;//输出原始观测数据已知点点号、高程for(int i=0;i<=yz_szd_zs-1;i++){outfile<<" "<<dm[i];outfile<<setiosflags(ios::fixed);outfile<<setw(10)<<setprecision(4)<<szd_gc[i]<<endl;}outfile<<endl<<endl<<"===================== 高差观测值====================="<<endl<<endl;//输出原始观测数据高程观测值与路线长outfile<<"起始点名"<<" "<<"终点点名"<<" "<<"高差观测值(m)"<<" "<<"两点间距离(km)"<<endl;for(i=0;i<=gcz_zs-1;i++){outfile<<" "<<dm[qsd_dh[i]]<<setw(9)<<dm[zd_dh[i]];outfile<<setiosflags(ios::fixed);outfile<<setw(16)<<setprecision(4)<<gcz[i];outfile<<setiosflags(ios::fixed);outfile<<setw(16)<<setprecision(4)<<1.0/P[i]<<endl;}m_mu=sqrt(pvv/(gcz_zs-(szd_zs-yz_szd_zs)));//计算单位权中误差outfile<<endl<<"===================== 单位权中误差====================="<<endl;//输出单位权中误差outfile<<endl<<"σ0="<<m_mu<<endl;outfile<<endl<<"===================== 高程平差值及其精度====================="<<endl<<endl;//输出高程平差值及其精度outfile<<"点名近似高程改正数高程平差值中误差"<<endl;for( i=0; i<szd_zs; i++){outfile<<setw(2)<<dm[i];double dx=dX[i];double qii=ATPA[ij(i,i)];outfile<<setiosflags(ios::fixed);outfile<<setw(12)<<setprecision(4)<<szd_gc[i]-dx;outfile<<setiosflags(ios::fixed);outfile<<setw(10)<<setprecision(4)<<dx;outfile<<setiosflags(ios::fixed);outfile<<setw(11)<<setprecision(4)<<szd_gc[i];outfile<<setiosflags(ios::fixed);outfile<<setw(10)<<setprecision(4)<<sqrt(qii)*m_mu<<endl;}//输出观测值平差值及其精度outfile<<endl<<endl<<"===================== 观测值平差值及其精度====================="<<endl<<endl;outfile<<"起点终点观测高差ｖ"<<" 高差平差值观测权中误差"<<endl;for(i=0;i<=gcz_zs-1;i++){int k1=qsd_dh[i];int k2=zd_dh[i];double qii=ATPA[ij(k1,k1)];double qjj= ATPA[ij(k2,k2)] ;double qij=ATPA[ij(k1,k2)];double ml=sqrt(qii+qjj-2.0*qij)*m_mu;outfile.width(2);outfile<<dm[k1];outfile.width(7);outfile<<dm[k2];outfile<<setiosflags(ios::fixed);outfile<<setw(12)<<setprecision(4)<<gcz[i];outfile<<setiosflags(ios::fixed);outfile<<setw(10)<<setprecision(4)<<V[i];outfile<<setiosflags(ios::fixed);outfile<<setw(10)<<setprecision(4)<<gcz[i]+V[i];outfile<<setiosflags(ios::fixed);outfile<<setw(10)<<setprecision(4)<<P[i];outfile<<setiosflags(ios::fixed);outfile<<setw(10)<<setprecision(4)<<ml<<endl;}outfile.close();}////////////////////////////////////////////////////////////////////////// // 水准网最小二乘平差void SZWPC::zxecpc(char *resultfile){ca_H0(); //近似高程计算ca_ATPA(); // 组成法方程ca_dX(); // 高程平差值计算}///////////////////////////////////////////////////////////////////////int main(){char *datafile ="算例\\Data.txt";//原始数据文件存储地址指针char *resultfile ="算例\\Result.txt";//平差结果输出地址指针cout<<endl<<endl<<"水准网经典间接平差"<<endl<<endl;cout<<"原数据文件位置："<<datafile<<endl;cout<<"平差结果文件位置："<<resultfile<<endl<<endl;SZWPC new_net;//定义新的对象new_net.inputdata(datafile);//输入原始数据new_net.zxecpc(resultfile);//最小二乘平差计算new_net.printresult(resultfile);//输出平差结果return 0;}。