1、 /观测值的权 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(); /近似高
2、程计算函数 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() if(gcz_zs0) delete qsd_dh; delete zd_dh; delete gcz; delete P; de
3、lete V; if(szd_zs delete szd_gc; delete ATPA; delete ATPL; delete dX; for(int i=0; i=j)? i*(i+1)/2+j :j*(j+1)/2+i;/ 对称正定矩阵求逆(仅存下三角元素)(参考他人)bool SZWPC:inverse(double a,int n) double *a0=new doublen; for(int k=0;kn;k+) double a00=a0; if(a00+1.0=1.0) delete a0; return false; for(int i=1;i double ai0 =
4、ai*(i+1)/2; if(i=n-k-1)a0i= -ai0/a00; else a0i= ai0/a00; for(int j=1;j=i;j+) a(i-1)*i/2+j-1=ai*(i+1)/2+j+ai0*a0j; for(i=1; a(n-1)*n/2+i-1=a0i; an*(n+1)/2-1=1.0/a00; delete a0; return true;/ 原始数据输入函数void SZWPC:inputdata(char *datafile) ifstream infile(datafile,ios:in); if(! infile) cerrgcz_zsszd_zsyz
5、_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; dmi =
6、NULL;/ dm_dh函数根据dmi是否为NULL确定dmi是否为点名地址 char buffer128; /临时数组,保存从文件中读到的点名 for( i=0;=yz_szd_zs-1;i+)/ 读取已知高程数据 infilebuffer; int c=dm_dh(buffer);szd_gci; for(i=0;gcz_zs;i+)/ 读取观测数据 /读取高程起点名 qsd_dhi=dm_dh(buffer);/读取高程终点 zd_dhi=dm_dh(buffer);gcziPi; /读取高差值与路线长度 Pi=1.0/Pi;/线路长转化为观测值的权 infile.close();/ 点
7、名转点号,返回点名对应的点号dm_dh(char *name) i+) if(dmi!=NULL) if(strcmp(name,dmi)=0)return i;/将待查点名与已经存入点名数组的点名比较,若存在返回点号 else int len = strlen(name);/判断点名长度 dmi = new charlen+1;/为点名申请存储空间 strcpy(dmi, name);/待查点名是一个新的点名,将新点名的地址放到dm数组中 return i;/返回点号 return -1; /dm数组已经存满,且没有待查点名/高程近似值计算ca_H0() for(int i=yz_szd_z
8、s;i+)szd_gci=-10000.9;/为计算机设置辨别未知高程点的标志 for(int j=1; int k=0; /计算出近似高程的点数 for(i=0; int k1=qsd_dhi; int k2=zd_dhi; if(szd_gck1-10000.0 & szd_gck2-10000.0)/k1点高程或高程近似值已知,k2点高程或高程近似值未知 szd_gck2=szd_gck1+gczi;/计算近似高程 k+; else if(szd_gck1-10000.0)/k2点高程或高程近似值已知,k1点高程或高程近似值未知 szd_gck1=szd_gck2-gczi; k+; i
9、f(k=(szd_zs-yz_szd_zs)break;/所有的近似高程计算完成,退出/ 组成法方程ca_ATPA() /int t=szd_zs;szd_zs*(szd_zs+1)/2; i+) ATPAi=0.0;/赋初值 i+) ATPLi=0.0; k k+) int i=qsd_dhk;/获取点号 int j=zd_dhk; double Pk=Pk;/获取权值 double lk=gczk-(szd_gcj-szd_gci);/获得第k个自由项 ATPLi-=Pk*lk;/获得法方程自由项 ATPLj+=Pk*lk; ATPAij(i,i)+=Pk;/获得法方程系数矩阵 ATPAi
10、j(j,j)+=Pk; ATPAij(i,j)-=Pk;/ 高程平差值计算ca_dX()i+) ATPAij(i,i)=1.0e30;/处理已知点inverse(ATPA,szd_zs)/矩阵求逆法方程系数矩阵降秩!/矩阵为奇异矩阵,无法求逆 exit(0);/退出程序 i+)/计算高程改正数 double xi=0.0; for(int j=0; j j+) xi+=ATPAij(i,j)*ATPLj; dXi=xi; szd_gci+=xi;/计算高程平差值/ 残差计算double SZWPC:ca_V() double pvv=0.0;=gcz_zs-1; int k1=qsd_dhi;
11、 int k2=zd_dhi; Vi=szd_gck2-szd_gck1-gczi; pvv+=Vi*Vi*Pi; return(pvv);/ 原始数据和平差值输出printresult(char *resultfile) double pvv=ca_V(); / 残差计算 ofstream outfile(resultfile,ios:out);/以输出方式打开文件,若文件不存在,创建文件 /输出原始观测数据 outfileendl观测总数:gcz_zs 总点数:已知点数:yz_szd_zs= 已知高程 =/输出原始观测数据已知点点号、高程 outfiledmi;setiosflags(io
12、s:fixed);setw(10)setprecision(4)szd_gci= 高差观测值=/输出原始观测数据高程观测值与路线长起始点名终点点名高差观测值(m)两点间距离(km)dmqsd_dhisetw(9)dmzd_dhi;setw(16)gczi;1.0/Pi m_mu=sqrt(pvv/(gcz_zs-(szd_zs-yz_szd_zs);/计算单位权中误差= 单位权中误差=/输出单位权中误差0=m_mu= 高程平差值及其精度 =/输出高程平差值及其精度点名 近似高程 改正数 高程平差值 中误差setw(2) double dx=dXi; double qii=ATPAij(i,i)
13、;setw(12)szd_gci-dx;dx;setw(11)sqrt(qii)*m_mu /输出观测值平差值及其精度= 观测值平差值及其精度 =起 点 终 点 观测高差 高差平差值 观测权 中误差 double qii=ATPAij(k1,k1); double qjj= ATPAij(k2,k2) ; double qij=ATPAij(k1,k2); double ml=sqrt(qii+qjj-2.0*qij)*m_mu; outfile.width(2);dmk1; outfile.width(7);dmk2;Vi;gczi+Vi;ml outfile.close();/ 水准网最小
14、二乘平差zxecpc(char *resultfile) ca_H0(); /近似高程计算 ca_ATPA(); / 组成法方程 ca_dX(); / 高程平差值计算 /int main() char *datafile =算例Data.txt/原始数据文件存储地址指针 char *resultfile =算例Result.txt/平差结果输出地址指针 cout水准网经典间接平差原数据文件位置:datafile平差结果文件位置:resultfile SZWPC new_net;/定义新的对象 new_net.inputdata(datafile);/输入原始数据 new_net.zxecpc(resultfile);/最小二乘平差计算 new_net.printresult(resultfile);/输出平差结果 return 0;
copyright@ 2008-2022 冰豆网网站版权所有
经营许可证编号:鄂ICP备2022015515号-1