gps.rar_GPS网平差_gps平差

//#define MAX 100 #define PI 3.14159265358979312 #define rou (180.0*60*60/PI) #include"MyClassLibrary0504.h" #include<iostream.h> //************************************************************************************************ class adj : public comadj{ // adjustment model : // V=AX-L P // A[m][n], P[m][m], L[m][1] public: adj(){} adj(int k,int w):comadj(k,w) {} void ksetadj(); int fsetadj(char *name) ; //file_data input int doadj(int known,int r); // 极大权法最小二乘平差,known--控制点已知数据个数; // r－固定数据个数＋测站数（平面网） int doadj(); int rubust(); void adjdis(); // 屏幕输出 int foutadj(char *name); // 文件输出 private: }; //******************************************************************************************* int adj::doadj() // 普通最小二乘平差 { MAT APA,AT; AT=A.T(); APA=AT*P*A; N_1=APA.inverse1(); int flag; if(APA.R()==APA.GetRow()) flag=1; else flag=0; if(flag!=1) { this->flag=0; return 0; } MAT AX=A.T()*P*l; X=N_1*AX; AX=A*X; for(int i=0;i<m;i++) V=AX-l; MAT VPV=V.T()*P*V; double cc=VPV.GetElem(0,0); this->m0=sqrt(cc/(m-n)); this->flag=1; return 1; } //************************************************************************************************ void adj::ksetadj() //keyboard and screen input { cout<<" input the object name: "; cin>>this->name;cout<<endl; cout<<" input XYobservation number: "; cin>>this->m;cout<<endl; cout<<" input unknown data number: "; cin>>n;cout<<endl; t=m-n; A.SetRow(m); A.SetRank(n); P.SetRow(m); P.SetRank(m); N_1.SetRow(n); N_1.SetRank(n); X.SetRow(n); X.SetRank(1); l.SetRow(m); l.SetRank(1); V.SetRow(m); V.SetRank(1); mX.SetRow(n); mX.SetRank(1); N.SetRow(n); N.SetRank(n); cout<<" input data of mat A["<<m<<"]["<<n<<"]"<<endl; for(int i=0;i<this->m;i++) for(int j=0;j<this->n;j++) {double cc; cin>>cc;A.SetElem(i,j,cc);} cout<<" input data of mat P["<<this->m<<"]["<<this->m<<"]"<<endl; for(i=0;i<this->m;i++) for(int j=0;j<this->m;j++) {double cc; cin>>cc;P.SetElem(i,j,cc);} cout<<" input data of mat L["<<this->m<<"][1]"<<endl; for(i=0;i<this->m;i++) {double cc; cin>>cc;l.SetElem(i,0,cc);} } //************************************************************************************************ int adj::fsetadj(char *name) //keyboard and screen input { ifstream in(name,ios::nocreate); if(!in) return 0; // 输入平差项目名 in>>this->name; // input XYobservation number: in>>this->m; // input unknown data number: in>>this->n; t=m-n; A.SetRow(m); A.SetRank(n); P.SetRow(m); P.SetRank(m); N_1.SetRow(n); N_1.SetRank(n); X.SetRow(n); X.SetRank(1); l.SetRow(m); l.SetRank(1); V.SetRow(m); V.SetRank(1); mX.SetRow(n); mX.SetRank(1); N.SetRow(n); N.SetRank(n); // input data of mat A[m][n]: for(int i=0;i<this->m;i++) for(int j=0;j<this->n;j++) {double cc; in>>cc;A.SetElem(i,j,cc);} // input data of mat P[m][m]: for(i=0;i<this->m;i++) for(int j=0;j<this->m;j++) {double cc; in>>cc;P.SetElem(i,j,cc);} // input data of mat L[m][1]: for(i=0;i<this->m;i++) {double cc; in>>cc;l.SetElem(i,0,cc);} in.close(); return 1; } //******************************************************************************************* int adj::doadj(int known,int r) // 极大权法最小二乘平差,known--控制点已知数据个数; { // r－固定数据个数＋测站数（平面网） MAT APA,AT; AT=A.T(); APA=AT*P*A; double add(0); for(int i=0;i<this->n;i++) if(add<=APA.GetElem(i,i)) add=APA.GetElem(i,i); add*=1000; for(i=0;i<known;i++) // 对已知点方法程系数阵的处理： APA.SetElem(i,i,APA.GetElem(i,i)+add); // 控制点点号对应矩阵主元加平均权的10000倍 N_1=APA.inverse1(); if(N_1.R()==N_1.GetRow()) flag=1; else flag=0; if(flag!=1) // 不可求逆的判断及处理 { this->flag=0; return 0;} MAT AX; // 极大权平差过程 AX=AT*P*l; X=N_1*AX; AX=A*X; V=AX-l; MAT VPV=V.T()*P*V; double cc=VPV.GetElem(0,0); // VTPV计算 this->m0=sqrt(cc/(this->m-this->n+known-r)); // 单位权中误差计算，必要观测数为（this->n-knownp） this->flag=1; // 平差完毕 return 1; } //******************************************************************************************* void adj::adjdis() { cout<<" 平差项目名 ："<<this->name<<endl<<endl; cout<<" 误差方程阵："<<endl; this->A.print(); cout<<endl<<endl<<" 观测值权矩阵："<<endl; this->P.print(); cout<<endl<<endl<<" 常数项："<<endl; l.print(); cout<<endl<<endl; cout<<endl<<endl<<" 未知数解："<<endl; X.print(); cout<<endl<<endl<<" 未知数协因数阵："<<endl; N_1.print(); cout<<endl<<endl<<" 改正数："<<endl; V.print(); cout<<endl<<endl<<" 单位权中误差：+-"<<this->m0<<endl; } //************************************************************************************************ int adj::foutadj(char *name) { ofstream on(name); if(!on) return 0; // 输出平差项目名 on<<this->name<<endl<<endl; // output XYobservation number: on<<this->m<<endl; // output unknown data number: on<<this->n<<endl<<endl; //out unknown data results for(int i=0;i<n;i++) on<<X.GetElem(i,0)<<endl; on<<endl; // 未知数协因数阵："<<endl; for(i=0;i<n;i++) { for(int j=0;j<n;j++) on<<N_1.GetElem(i,j)<<" "; on<<endl; } on<<endl; //output 改正数：; for(i=0;i<m;i++) on<<V.GetElem(i,0)<<endl; on<<endl; on<<m0<<endl; on.close(); return 1; } //************************************************************************************************ int adj :: rubust() // 定义抗差估计 { // 1.最小二乘平差 if(!doadj()) return 0; double *xold,*xnew,small(0),n0(0); MAT P(m,1); xold=new double[n]; int num(0); xnew=new double[n]; for(int i=0;i<m;i++) P.SetElem(i,0,this->P.GetElem(i,i)); do{ small=0; for(i=0;i<n;i++) *(xold+i)=X.GetElem(i,0); // 2.等权处理 for(i=0;i<m;i++) { if(fabs(V.GetElem(i,0))>2.5*m0) { this->P.SetElem(i,i,0);n0++;} else if(fabs(this->V.GetElem(i,0))<1.8*this->m0*sqrt(1/this->P.GetElem(i,i))) ; else this->P.SetElem(i,i,this->P.GetElem(i,i)/1.8/this->m0); } cout<<"this->m0="<<this->m0<<endl; P.print(); doadj(); for(i=0;i<this->n;i++) { *(xnew+i)=this->X.GetElem(i,0)-*(xold+i); if(small<fabs(*(xnew+i))) small=fabs(*(xnew+i)); } num++; }while(small>0.0000000001); if(n0) this->m0*=sqrt((this->m-this->n)/(this->m-this->n-n0)); delete [] xold; delete [] xnew; return 1; } //********************************************************************************** //**************** 高程网平差程序 **************************************** class Hpnt{ // 高程点结构 public: char name[20]; double H; double H0; double mH; int fixed; // known=1;else =0; int i; // point number }; class line{ // 水准线路结构 public: Hpnt *startp; Hpnt *endp; double length; double h; int style; // 水准测量=1;三角高程=2 }; class Hnet{ // 高程网结构 public: Hnet(char *fname); // 文件输入高程网函数 Hnet(){obnum=allpnum=fixpnum=0;} Hnet(int allp,int fixp,int obn):aa(obn,allp-fixp) { allpnum=allp; fixpnum=fixp; obnum=obn; Pt=new Hpnt[allpnum]; L=new line[obnum]; }