水准网间接平差程序设计(C++)
- 格式:doc
- 大小:54.00 KB
- 文档页数:10
水准网平差c++代码水准平差网结果#include<iostream>#include<stdlib.h>#include<iomanip>#include<math.h>#define max 50Using namespace std;Class CMatrix{public:CMatrix(){row=0; column=0;}; // 默构造函数认认认认认CMatrix(int i, int j){row=i;column=j;} // 构造函数一CMatrix(constCMatrix& m); // 认制构造函数~CMatrix(void){/*cout<<"认认认认认认认认认认认使用,矩所占空以放:"<<endl;*/} // 默析构函数认认认认认CMatrix& operator=(constCMatrix& m); // 认认运算符bool operator==(constCMatrix& m); // 比括运算符认认认认bool operator!=(constCMatrix& m); // 比括运算符认认认认CMatrix operator+(constCMatrix& m); // 加运算符CMatrix operator-(constCMatrix& m); // 减运算符CMatrix& operator+=(constCMatrix& m); // 自加运算符CMatrix& operator-=(constCMatrix& m); // 自减运算符CMatrix operator-();// 取数认认CMatrix& operator*(constCMatrix& m); // 乘法运算符void input(); //认认入矩void outputMatrix(); // 认认认认出矩void setValue(int row, int column, double value) { A[row-1][column-1] = value; }// 认置(i,j)的认double getValue(int row, int column) const { return A[row-1][column-1]; }// 认认置行、列的voidsetRow(constint i) { row = i; }intgetRow() const { return row; }voidsetColunm(constint i) { column = i; }intgetColumn() const { return column; }CMatrix& change(int i, int j);//交矩的行认认认认认CMatrix& transpose(); // 矩置认认认CMatrix& inverse(); // 矩求逆认认认void find(int& f)const;// 判断矩是否可用于迭代求解认认认认认认认认认认认认friend void jocabi(constCMatrix& a) ; //迭代求解void lzys(); //列主元素法求解void solve(); //可逆性矩求解认认认认认认void qxnh(); //曲合认认认private:// 成量认认认double A[max][max];int row;// 行int column;// 列};void CMatrix::input() //认入{ cout<<"认认认认认认认始入矩:"<<endl;int i, j;double z;for(i=0;i<row;i++){ cout<<"认认入第"<<i+1<<"行的:认认"<<endl;for(j=0;j<column;j++){cin>>z;A[i][j]=z;}}cout<<endl;};CMatrix::CMatrix(constCMatrix& m) // 认制构造函数{ int i, j; for(i=0;i<m.row;i++)for(j=0;j<m.column;j++)this->A[i][j]=m.A[i][j];};CMatrix&CMatrix::operator=(constCMatrix& m) // 认认运算符{ inti,j;for(i=0;i<row;i++){for(j=0;j<column;j++)A[i][j]=m.A[i][j];}return *this;};boolCMatrix::operator ==(constCMatrix& m) // 比括运算符认认认认{ inti,j,k;for(i=0;i<m.row;i++){ for(j=0;j<m.column;j++)if(this->A[i][j]=m.A[i][j]) k=1;else k=0;}if(k=1) return true;else return false;};boolCMatrix::operator !=(constCMatrix& m) // 比括运算符认认认认{ inti,j,k;for(i=0;i<m.row;i++){ for(j=0;j<m.column;j++)if(this->A[i][j]=m.A[i][j]) k=1;else k=0;}if(k=0) return true;else return false;};CMatrixCMatrix::operator+(constCMatrix& m)// 加运算符{ inti,j; if((this->row==m.row)&&(this->column==m.column)){ for(i=0;i<m.row;i++)for(j=0;j<m.column;j++)this->A[i][j]+=m.A[i][j];}else { cout<<"此两矩不能相加,:认认认认认认认认认认"<<endl; } return *this;};CMatrixCMatrix::operator-(constCMatrix& m)// 减运算符{ inti,j;if((this->row==m.row)&&(this->column==m.column)){ for(i=0;i<m.row;i++)for(j=0;j<m.column;j++)this->A[i][j]-=m.A[i][j];}else { cout<<"此两矩不能相加,:认认认认认认认认认认"<<endl; }return *this;};CMatrix&CMatrix::operator+=(constCMatrix& m) //自加运算符{ inti,j;for(i=0;i<m.row;i++)for(j=0;j<m.column;j++)this->A[i][j]=2*m.A[i][j];return *this;};CMatrix&CMatrix::operator-=(constCMatrix& m) //自减运算符{ inti,j;for(i=0;i<m.row;i++)for(j=0;j<m.column;j++)this->A[i][j]=m.A[i][j]-m.A[i][j];return *this;};voidCMatrix::find(int& f)const{ int i;for(i=0;i<this->row;i++)if(this->A[i][i]!=0) f=1;else f=0;};CMatrixCMatrix::operator-() // 取数认认{ inti,j;for(i=0;i<this->row;i++)for(j=0;j<this->column;j++)this->A[i][j]=-this->A[i][j];return *this;};CMatrix&CMatrix::operator*(constCMatrix& m) // 乘法运算符{ inti,j,t;CMatrix n;if(this->column==m.row){for(i=0;i<this->row;i++)for(j=0;j<m.column;j++){double sum=0.0;for(t=0;t<m.row;t++)sum+=this->A[i][t]*m.A[t][j];n.A[i][j]=sum;}}else {cerr<<"此两矩不能相乘,:认认认认认认认认认认"<<endl; exit(1);} return n;};void CMatrix::outputMatrix()// 认认认认出矩{ inti,j;for(i=1;i<=row;i++){ for(j=1;j<=column;j++)cout<<A[i-1][j-1]<<" ";cout<<endl;}};CMatrix&CMatrix::transpose()// 矩置认认认{ inti,j;CMatrix m(this->column,this->row);for(i=0;i<this->row;i++)for(j=0;j<this->column;j++){ m.A[j][i]=this->A[i][j];}return m;};void jocabi(constCMatrix& a) //高斯迭代求解{ int f=1;a.find(f);if(f==0) {cerr<<"认认认认认认认认认认认认矩不足迭代求解条件:"<<endl; exit(1); }else{CMatrixx,w;x.setColunm(1);x.setRow(a.getColumn());w.setColunm(1);w.setRow(a.getColumn());int i;double z;for(i=1;i<=x.row;i++){ cout<<"认认认认认入等式右的第"<<i<<"个:认认"; cin>>z;w.A[i-1][0]=z;}for(i=1;i<=x.row;i++){ cout<<"认认入X"<<i<<"的近似:认认";cin>>z;x.setValue(i, 1, z);}i=1;while(i<=20){ int j, k;for(j=1;j<=x.row;j++){ double sum=0.0;for(k=1;k<j;k++){sum=sum-(a.A[j-1][k-1])*(x.A[k-1][0]);}for(k=j+1;k<=x.row;k++){sum=sum-(a.A[j-1][k-1])*(x.A[k-1][0]);}sum+=w.A[j-1][0];x.A[j-1][0]=sum/a.A[j-1][j-1];}i++;}for(i=1;i<=x.row;i++)cout<<"X"<<i<<" = "<<x.A[i-1][0]<<endl;}};double jdz(double a) //认认认{ if(a>0.0) return a;else return -a;};CMatrix&CMatrix::change(int i, int j)//交矩的行认认认认认{ int k;double z;for(k=1;k<=this->column;k++){z=A[i-1][k-1];A[i-1][k-1]=A[j-1][k-1];A[j-1][k-1]=z;}return *this;};void CMatrix::lzys() //列主元素法求解{ CMatrixx,w;x.setColunm(1);x.setRow(getColumn());w.setColunm(1);w.setRow(getColumn());int i;double z;for(i=0;i<row;i++){for(int j=0;j<row; j++)w.A[i][j]=0.0;}for(i=1;i<=x.row;i++)x.setValue(i, 1, 0.0);for(i=1;i<=x.row;i++){ cout<<"认认认认认入等式右的第"<<i<<"个:认认"; cin>>z;w.A[i-1][0]=z;}i=0;while(i<x.row-1){ intj,t,h=i;for(j=i;j<x.row-1;j++){if(jdz(A[j][i])<jdz(A[j+1][i])) h=j+1; }this->change(i+1,h+1);w.change(i+1,1);for(j=i+1;j<row;j++){ double k;k=A[j][i]/A[i][i];for(t=i;t<column;t++){A[j][t]=A[j][t]-k*A[i][t];}w.A[j][0]=w.A[j][0]-k*w.A[i][0];}i++;}if(this->A[row-1][column-1]==0){cerr<<"此矩的方程有无解:认认认认认认认认认认认认"<<endl; exit(1);} i=column-1;while(i>0){ int j, t;for(j=i-1;j>=0;j--){ double k;k=A[j][i]/A[i][i];for(t=i;t>j;t--){ A[j][t]=A[j][t]-k*A[i][t]; }w.A[j][0]=w.A[j][0]-k*w.A[i][0];}i--;}int j;for(j=0;j<x.row;j++){x.A[j][0]=w.A[j][0]/this->A[j][j];}for(i=1;i<=x.row;i++)cout<<"X"<<i<<" = "<<x.A[i-1][0]<<endl;};CMatrix&CMatrix::inverse() //矩求逆认认认{ if(this->row!=this->column) {cerr<<"认认认认认认认认认认矩不符合求逆条件:"<<endl; exit(1);}else{CMatrix w;w.setColunm(this->getRow());w.setRow(getColumn());int i;for(i=0;i<row;i++){for(int j=0;j<row; j++){w.A[i][j]=0.0;}w.A[i][i]=1;}i=0;while(i<column-1){ intj,t,h=i;for(j=i;j<row-1;j++){if(jdz(A[j][i])<jdz(A[j+1][i]))h=j+1;}this->change(i+1,h+1);w.change(i+1,h+1);for(j=i+1;j<row;j++){ double k;k=A[j][i]/A[i][i];for(t=i;t<column;t++){A[j][t]=A[j][t]-k*A[i][t];w.A[j][t]=w.A[j][t]-k*w.A[i][t];}}i++;}if(this->A[row-1][column-1]==0){cerr<<"此矩求逆不成功,其所的方程有无解:认认认认认认认认认认认认认认认认认认认认"<<endl; exit(1);}i=column-1;while(i>0){ int j, t;for(j=i-1;j>=0;j--){ double k;k=A[j][i]/A[i][i];for(t=i;t>j;t--){ A[j][t]=A[j][t]-k*A[i][t]; }for(t=column-1;t>=0;t--){w.A[j][t]=w.A[j][t]-k*w.A[i][t];}}i--;}intj,k;for(j=0;j<row;j++){for(k=0;k<w.column;k++){w.A[j][k]=w.A[j][k]/this->A[j][j];}this->A[j][j]=this->A[j][j]/this->A[j][j];}return w;}};void CMatrix::solve() //可逆性矩求解认认认认认认{CMatrixx,w,c;x.setColunm(1);x.setRow(getColumn());w.setColunm(1);w.setRow(getColumn());c.setColunm(getColumn());c.setRow(getRow());int i;double z;for(i=1;i<=x.row;i++){ cout<<"认认认认认入等式右的第"<<i<<"个:认认";cin>>z;w.A[i-1][0]=z;}c.operator=(this->inverse());x.operator=(c.operator *(w)); cout<<"求解果:认认认"<<endl;for(i=1;i<=x.row;i++)cout<<"X"<<i<<" = "<<x.A[i-1][0]<<endl;};void CMatrix::qxnh() //曲合认认认{cout<<"用矩曲合:认认认认认认"<<endl;int i, j,k,t;double x, y;cout<<"认认认认认认认认入已知坐点的个数:";cin>>i;CMatrix a(i,1), g(i,i) , w(i,1),c(i,i),q(i,i),p(i,i);for(j=1;j<=i;j++){ cout<<"认认入第"<<j<<"个已知点坐认(先横坐后坐认认认认认):"; cin>>x>>y;w.A[j-1][0]=y;g.A[j-1][0]=1.0;for(k=1;k<i;k++){ g.setValue(j,k+1,1.0);for(t=1;t<=k;t++)g.A[j-1][k]*=x;}}c=g.transpose();q=c*g;p=q.inverse();w=c*w;a=p*w;cout<<"Y="<<a.A[0][0];for(j=1;j<i;j++){ cout<<'+'<<'('<<a.A[j][0]<<')';for(k=0;k<j;k++)cout<<"*x";}cout<<endl;cout<<" 1 .求取其它点坐认 0. 退出 "<<endl; while(1){ cout<<"认认认:";cin>>k;switch(k){ case 1:{ double x, y,sum;cout<<"认认认认入待求点的横坐,x,:";cin>>x;y=a.A[0][0];for(j=1;j<i;j++){sum=a.A[j][0];for(k=0;k<j;k++){sum*=x;}y+=sum;}cout<<"认认认坐:y="<<y<<endl;break;}case 0:{cerr<<"认认使用:"<<endl; exit(1);break;}default: cout<<"认认认认认认认认有,:"<<endl; break; }}};structCElvDif{doublevalue;// 认认认double weight;// 认重longstartPoint;// 起始点号认认longendPoint;// 认认认点号};//水准点的认认认认structCLevelPoint{// 高程平差=认高程认+高程改正数认认认认long index;// 水准点号认认double eleValue;// 高程认double dv;// 高程改正数,初始化认认认认认认认认认0, boolisKnown;// 是否已知点认认认认};classCElevationNet{public:// 成函数认认认CElevationNet(){numElvDif=0;numPoints=0;numKnPoint=0;};// 构造函数~CElevationNet(){};// 析构函数void input();void output();voidjsgc();intgetgz(){return numElvDif;} //返回数目认认认认认intgetzs(){return numPoints;} //返回点数认认认intgetys(){return numKnPoint;} //返回已知点数intgetws(){return numPoints-numKnPoint;} //返回未知点数double geteleValue(int i){return this->lpVec[i-1].eleValue;} //认认得高程long getindex(int i){return lpVec[i-1].index;} //返回高程点号认认void seteleValue(int i, double value){lpVec[i-1].eleValue+=value;} //修改高程认void setdv(inti,double value){this->lpVec[i-1].dv=value;} //修改改正数double getdv(int i){return lpVec[i-1].dv;} //返回改正数friend void xishu(CMatrix& B, CMatrix&X,CElevationNet A); //求取系数矩和未知点高程矩认认认认认认认认认friend void quanzhen(CMatrix& Q, CElevationNet A); //求取认认friend void l_zhen( CMatrix& l ,CMatrix& L, CElevationNet A); //求取认认认L矩和认认l认private:// 成量认认认intnumElvDif;// 高差数认认intnumPoints;// 控制网中点的数目认认认intnumKnPoint;//控制网中已知点的数目CElvDifedVec[max];// 认认认认数CLevelPointlpVec[max]; // 高程数认认认};voidCElevationNet::input(){ inta,b,c,i;doublez,l,g;cout<<"认认认认认认认认认认认认认认认认认认认认认认认认认已知点和未知点行号,注意已知点的号小于未知点的号:"<<endl<<endl;cout<<"认认认认认认认认认认入控制网中水准点数和已知点数:";cin>>a>>b;numPoints=a;numKnPoint=b;for(i=0;i<numPoints;i++){ lpVec[i].index=i+1;lpVec[i].eleValue=0.0;lpVec[i].dv=0.0;lpVec[i].isKnown=0;}for(i=0;i<numKnPoint;i++){ cout<<"认认入第"<<i+1<<"个已知点高程:认认";cin>>z;lpVec[i].eleValue=z;lpVec[i].isKnown=1;}cout<<"认认认认认认认认入个数:";cin>>c;numElvDif=c;for(i=0;i<numElvDif;i++){ cout<<"认认入第"<<i+1<<"段段的高差认认认认认认认(认位米)、度认认(认位千米)、起始点号和点号:认认认认认认认认"<<endl;cin>>g>>l>>a>>b;edVec[i].value=g;edVec[i].weight=l;edVec[i].startPoint=a;edVec[i].endPoint=b;}};voidCElevationNet::output(){int i;for(i=0;i<this->numPoints;i++)cout<<"认认号"<<this->lpVec[i].index<<"的点的高程:认认认"<<this->lpVec[i].eleValue<<endl;for(i=0;i<this->numElvDif;i++)cout<<"认认段"<<i+1<<"的平差后的:认认认"<<this->edVec[i].value<<endl;};void CElevationNet::jsgc(){ inti,j;for(i=0;i<this->numElvDif;i++){ if((lpVec[edVec[i].startPoint-1].eleValue!=0.0)&&(lpVec[edVec[i].endPoint-1].eleValue==0.0)) {lpVec[edVec[i].endPoint-1].eleValue=(lpVec[edVec[i].startPoint-1].eleValue)+(edVec[i].value);}else if(lpVec[edVec[i].startPoint-1].eleValue==0.0&&lpVec[edVec[i].endPoint-1].eleValue!=0.0) {lpVec[edVec[i].startPoint-1].eleValue=(lpVec[edVec[i].endPoint-1].eleValue)-(edVec[i].value);}}for(j=numKnPoint;j<this->numPoints;j++){if(this->lpVec[j].eleValue==0.0){for(i=0;i<this->numElvDif;i++){ if(this->edVec[i].startPoint==lpVec[j].index&&this->lpVec[this->edVec[i].endPoint-1].eleValue!=0.0)this->lpVec[j].eleValue=this->lpVec[this->edVec[i].endPoint-1].eleValue-this->edVec[i].value;else if(this->edVec[i].endPoint==lpVec[j].index&&this->lpVec[this->edVec[i].startPoint-1].eleValue!=0.0){this->lpVec[j].eleValue=this->lpVec[this->edVec[i].startPoint-1].eleValue+this->edVec[i].value;}}}//cout<<"第"<<this->lpVec[j].index<<"个点的高程近似:认认"<<this->lpVec[j].eleValue<<endl;}};voidxishu(CMatrix& B, CMatrix&X,CElevationNet A){inti,j;for(i=1;i<=B.getRow();i++)for(j=1;j<=B.getColumn();j++){ B.setValue(i,j,0);}for(i=0;i<A.numElvDif;i++){ if(A.lpVec[A.edVec[i].startPoint-1].isKnown==1&&A.lpVec[A.edVec[i].endPoint-1].isKnown==0) {B.setValue(i+1,A.edVec[i].endPoint-A.numKnPoint,1);}else if(A.lpVec[A.edVec[i].startPoint-1].isKnown==0&&A.lpVec[A.edVec[i].endPoint-1].isKnown==1) { B.setValue(i+1,A.edVec[i].startPoint-A.numKnPoint,-1);} else if(A.lpVec[A.edVec[i].startPoint-1].isKnown==0&&A.lpVec[A.edVec[i].endPoint-1].isKnown==0) { B.setValue(i+1,A.edVec[i].endPoint-A.numKnPoint,1);B.setValue(i+1,A.edVec[i].startPoint-A.numKnPoint,-1);}}for(i=0;i<X.getRow();i++)X.setValue(i+1,1,A.lpVec[A.numKnPoint+i].eleValue);};void quanzhen(CMatrix& Q, CElevationNet A){inti,j;double sum=0.0;for(i=0;i<A.numElvDif;i++)sum+=A.edVec[i].weight;sum=sum/A.numElvDif;for(i=1;i<=Q.getRow();i++)for(j=1;j<=Q.getColumn();j++){ Q.setValue(i,j,0);}for(i=1;i<=Q.getRow();i++)Q.setValue(i,i,sum/A.edVec[i-1].weight);};void l_zhen( CMatrix& l ,CMatrix& L, CElevationNet A){int i;double m;for(i=0;i<A.numElvDif;i++)L.setValue(i+1,1,A.edVec[i].value);for(i=0;i<l.getRow();i++){m=(A.lpVec[A.edVec[i].startPoint-1].eleValue+A.edVec[i].value-A.lpVec[A.edVec[i].endPoint-1].eleValue)*1000;l.setValue(i+1,1,m);}};void main(){CElevationNet A;A.input();A.jsgc();CMatrixB(A.getgz(),A.getws()),Q(A.getgz(),A.getgz()),X(A.getws(),1),l(A.get gz(),1),b(A.getws(),A.getgz()),V(A.getgz(),1);CMatrixNBB(A.getws(),A.getws()),x(A.getws(),1),N(A.getws(),A.getws()),L(A.g etgz(),1),M(A.getws(),A.getgz()),R(A.getws(),A.getgz()), W(A.getws(),A.getgz());CMatrixvz(1,A.getgz()),vp(1,A.getgz()),p(1,1);xishu(B,X,A);quanzhen(Q,A);l_zhen(l,L,A);b.operator=(B.transpose());M.operator=(b.operator *(Q));NBB.operator=(M.operator *(B));N.operator=(NBB.inverse());R.operator=(N.operator *(b));W.operator=(R.operator *(Q));x.operator=(W.operator *(l));inti,j;for(i=A.getys()+1;i<=A.getzs();i++){A.setdv(i,(x.getValue(i-A.getys(),1)/1000));A.seteleValue(i,A.getdv(i));cout<<"平差后号认认认"<<A.getindex(i)<<"的点高程:认认"<<A.geteleValue(i)<<"米"<<endl;}V.operator=(B.operator *(x));V.operator-(l);for(j=0;j<V.getRow();j++)cout<<"第"<<j+1<<"段高差平差:认认认认认认认认认"<<(L.getValue(j+1,1)+V.getValue(j+1,1)/1000)<<"米"<<endl;vz.operator=(V.transpose());vp=(vz.operator *(Q));p=(vp.operator *(V));cout<<"平差后位中差:认认认认认认认认"<<"+_"<<sqrt(p.getValue(1,1)/A.getws())<<"毫米"<<endl;}。
水准网间接平差程序设计水准网间接平差是测量水准网中各测站的高程值,通过观测值的处理,进行计算来消除观测误差,得到准确的高程数据。
在进行水准网间接平差程序设计时,需要考虑观测值的处理方法、具体的计算步骤、误差的传递和消除等因素。
下面将详细介绍水准网间接平差程序设计的内容。
首先,在水准网间接平差的程序设计中,需要对观测值进行处理。
观测值的处理包括检查观测数据的精度、合理性及完整性,并进行数据的筛选和滤波处理。
在这一步骤中,需要使用适当的统计方法对观测数据进行筛选,剔除异常值和明显错误的数据,保留符合要求的观测值。
接下来,在进行水准网间接平差计算之前,需要对网络进行拟合,拟合过程即将观测值与已知高程值进行比较,并进行拟合计算得到误差。
网络拟合可以使用最小二乘法进行计算,即通过最小化观测值与已知高程值的差的平方和,来求得最优拟合结果。
然后,进行水准网的平差计算。
平差计算是根据测站之间的观测关系,通过一系列的计算公式,将所有观测值联立起来,并通过方程组进行求解,得到最终的平差结果。
在这个过程中,需要进行传递误差的计算,即通过误差传递公式计算各点高程值的精度,以评估平差结果的可靠性。
最后,在完成水准网间接平差计算之后,需要对平差结果进行检查和评估。
检查结果是否符合工程要求和精度要求,评估平差的可靠性。
如果结果不符合要求,需要重新进行观测值的处理和计算。
在进行水准网间接平差程序设计时,还需要注意以下几点:1.数据的输入与输出:程序需要提供方便的数据输入和输出方式,以便用户输入观测数据,并输出平差结果。
同时,需要考虑数据的存储和传输方式,确保数据的安全和完整性。
2.程序的可扩展性:设计程序时应考虑未来可能的数据规模扩大和功能的增加。
通过模块化设计和灵活的架构,使程序能够方便地扩展和添加新的功能。
3.用户友好性:程序应提供简单易用的操作界面,提供友好的用户交互方式。
用户应能够方便地输入观测数据和设置计算参数,并能够直观地查看和分析计算结果。
地球科学与环境工程学院水准间接平差实验报告书课程名:《误差理论与测量平差基础》学号:姓名:黄黎东指导老师:日期: 2015年12月7日一、任务概述利用MATLAB或者C++编程间接平差程序,通过该程序读取观测数据文件,并计算出平差结果。
二、计算结果截图:图一图二图三图四三、水准网图四、输入的数据格式数据格式为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);%获得观测个数nt=f(1,3);%获得必要观测个数ty=f(1,4);%获得已知点个数yXX=zeros(point,1);%初始化XX阵等于0,方便下面把已知点高程和未知点参数估值放到XX阵B=zeros(n,t);%初始化B阵,方便下面求V=Bx-l中的系数阵B;for j=1:yXX(j,1)=f(j+1,2);%把已知点高程放到XX阵中enddata=f((2+y):end,:);%从文件中获取观测数据,并放到data阵中h=data(:,3);%从data中获取观测高差,并放到h阵中P=zeros(n);%初始化权阵Pfor j=1:nP(j,j)=10/data(j,4);%以10km观测值为单位权误差计算权阵Pendfor 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个参数估值endendl=zeros(n,1);%初始化小l阵,方便下面求V=Bx-l中的系数阵l;for i=1:n%通过循环求小lpoint1=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('协方差阵:') Dxxdisp('单位权方差:') a0disp('协因数阵:') QxxBlPNWend。
《测绘程序设计》上机实验报告(Visual C++.Net)班级:测绘0901班学号: 04姓名:代娅琴2012年4月29日实验八平差程序设计基础一、实验目的巩固过程的定义与调用巩固类的创建与使用巩固间接平差模型及平差计算掌握平差程序设计的基本技巧与步骤二、实验内容水准网平差程序设计。
设计一个水准网平差的程序,要求数据从文件中读取,计算部分与界面无关。
1.水准网间接平差模型:2.计算示例:近似高程计算:3.水准网平差计算一般步骤(1)读取观测数据和已知数据;(2)计算未知点高程近似值;(3)列高差观测值误差方程;(4)根据水准路线长度计算高差观测值的权;(5)组成法方程;(6)解法方程,求得未知点高程改正数及平差后高程值;(7)求高差观测值残差及平差后高差观测值;(8)精度评定;(9)输出平差结果。
4.水准网高程近似值计算算法5.输入数据格式示例实验代码:#pragma onceclass LevelControlPoint{public:LevelControlPoint(void);~LevelControlPoint(void);public:CString strName;trName=pstrData[0];m_pKnownPoint[i].strID=pstrData[0];m_pKnownPoint[i].H=_tstof(pstrData[1]);m_pKnownPoint[i].flag=1;trName=pstrData[i];m_pUnknownPoint[i].strID=pstrData[i];m_pUnknownPoint[i].H=0;lag=0;pBackObj=SearchPointUsingID(pstrData[0]);pFrontObj=SearchPointUsingI D(pstrData[1]);ObsValue=_tstof(pstrData[2]);ist=_tstof(pstrData[3]);trID==ID){return &m_pKnownPoint[i];}}return NULL;}trID==ID){return &m_pUnknownPoint[i];}}return NULL;}LevelControlPoint* AdjustLevel::SearchPointUsingID(CString ID){LevelControlPoint* cp;cp=SearchKnownPointUsingID(ID);if(cp==NULL)cp=SearchUnknownPointUsingID(ID);return cp;}void AdjustLevel::ApproHeignt(void)lag!=1){pFrontObj->strID==m_pUnknownPoint[i].strID)&& m_pDhObs[j].cpBackObj->flag==1 ){ =m_pDhObs[i].cpBackObj->H - m_pDhObs[i].ObsValue;*/m_pUnknownPoint[i].H=m_pDhObs[j].cpBackObj->H + m_pDhObs[j].HObsValue;m_pUnknownPoint[i].flag=1;break;}}if(m_pUnknownPoint[i].flag!=1)pBackObj->strID==m_pUnknownPoint[i].strID)&& m_pDhObs[j].cpFrontObj->flag==1 ){ =m_pDhObs[j].cpFrontObj->H-m_pDhObs[j].HObsValue;/* m_pUnknownPoint[i].H=m_pDhObs[i].cpFrontObj->H+m_pDhObs[i].ObsValue;*/ m_pUnknownPoint[i].flag=1;break;}}}}if(i==m_iUnknownPointCount-1)lag!=1)ist);p(i,i)=value;}return p;}void AdjustLevel::FormErrorEquation(CMatrix &B, CMatrix &L){(m_iDhObsCount,m_iUnknownPointCount);(m_iDhObsCount,1);for(int i=0;i<m_iDhObsCount;i++)pBackObj->strID);tmpFront=SearchPointUsingID(m_pDhObs[i].cpFrontObj->strID);trID==tmpBack->strID)trID==tmpFront->strID)bsValue-(m_pDhObs[i].cpBackObj->H-m_pDhObs[i].cpFrontO bj->H);*/L(i,0)=m_pDhObs[i].HObsValue-(m_pDhObs[i].cpFrontObj->H - m_pDhObs[i].cpBackObj->H);(_T("%.3f"),L(i,0));L(i,0)=_tstof(tmp);L(i,0)=L(i,0)*1000;+=x(i,0);xt"));xt"));if()==IDCANCEL) return;CString strFileName=();setlocale(LC_ALL,"");CStdioFile sf;if(!(strFileName, CFile::modeCreate|CFile::modeWrite)) return;(LevleContent);();UpdateData(FALSE);}void CIndircLelveDlg::OnBnClickedComputelevel(){f\r\n"), [i].strID,[i].H);LevleContent+=Temp;}(_T("单位权中误差:%.1f mm\r\n"),r0*1000);LevleContent+=Temp;LevleContent+=_T("未知点高程中误差(mm):\r\n");for(int i=0;i< ;i++){();(_T("%s,%.1f\r\n"),[i].strName,Qx[i]*1000);LevleContent+=Temp;}UpdateData(false);}void CIndircLelveDlg::OnBnClickedSavelevleresult(){xt"));if()==IDCANCEL) return;CString strFileName=();setlocale(LC_ALL,"");CStdioFile sf;if(!(strFileName, CFile::modeCreate|CFile::modeWrite)) return;(LevleContent);();UpdateData(FALSE);}三、实验结果打开文件数据:平差结果:四、实验心得这从实验是我们测绘程序设计的最后一次实验,虽然这个学期我们做了好几次相关的实验,但是我却发现自己学的东西也越来越模糊,感觉很多内容都不理解。
水准网平差程序设计水准测量是现今城市测量工作中测定高程和对建筑物进行变形观测常用的测量方法。
水准网平差的目的在于依据最小二乘原理,消除观测之间的矛盾和不符值,进而求出点的最后高程及评定精度,编制相应的计算程序,实现水准网的自动平差,才能提高计算效率,然而近似高程的自动推算、闭合差的计算、以及误差方程、法方程的列立解算是程序设计的重点和难点。
标签:水准网平差;高程;程序设计1 数据库的建立1.1 建立建网信息表格建网信息表用于保存建网信息。
1.2 建立观测高差录入表格记录高差和高差起点终点点号和两点之间的距离。
1.3建立已知高程录入表格记录高程用于以后计算。
1.4建立平差成果表用于记录和输出结果。
2 应用软件用户界面的设计本程序共设计了五个窗体,分别为主窗体,建网信息输入窗体,测量高差输入窗体,已知高程输入窗体,平差计算结果显示窗体。
其中主窗体是程序的启动窗体,所有窗体都从主窗体界面弹出。
建网信息输入窗体用于输入项目名称,测量单位,负责人,测量范围,备注信息,并于数据库连接。
测量高差输入窗体用于输入高差观测值和路线长度,并和数据库连接。
已知高程输入窗体用于输入已知高程,并连接数据库。
平差计算结果显示窗体聚集了大部分的算法,包括待定点高程计算,误差方程和法方程的列立,法方程的运算等等,并通过DataGrid控件输出。
3 程序算法的基本思想及部分代码3.1 高程控制网间接平差的步骤(1)计算待定点的近似高程X=(x1,x2 …)T;(2)列出误差方程;(3)组成法方程;(4)解算法方程,求得dX=(dx1,dx2 …)T;(5)求得平差后的高程X=X+dX;(6)精度评定。
3.2 网形的编号及部分变量的定义3.2.1 网形编号为了使编程更加方便,首先约定各高程点编号由小到大按自然数的顺序编码,不可重复也不可缺少。
已知点优先于待定点,靠近已知点的待定点编号要尽量小。
3.2.2 部分变量定义Option Base 1Public IKP As Integer’已知点点个数;Public IUP As Integer’未知点点个数;Public IQ As Integer’总点个数;Public IOH As Integer’高差观测数;Public IZ As Integer’必要观测数。
水准网平差程序设计水准网平差是测绘工程中的一项重要工作,它涉及到对水准测量数据进行处理,以达到测量精度的要求。
水准网平差程序设计通常需要考虑以下几个方面:1. 程序设计的目的和意义水准网平差程序设计的主要目的是通过计算机软件对水准测量数据进行自动化处理,提高数据处理的效率和准确性。
这对于大型工程测量、城市基础设施建设、土地管理等领域具有重要意义。
2. 程序设计的基本要求- 准确性:程序需要能够准确地处理水准测量数据,减少人为误差。
- 稳定性:程序在运行过程中应具有较高的稳定性,避免因系统崩溃等原因导致数据丢失。
- 用户友好性:程序应具备良好的用户界面,使得非专业用户也能方便地使用。
- 扩展性:程序设计应考虑未来可能的功能扩展,以适应不断变化的测量需求。
3. 程序设计的理论基础水准网平差程序设计的理论基础主要包括:- 水准测量原理:了解水准测量的基本原理,包括视线高、转点高、已知点高程等概念。
- 误差理论:掌握测量误差的来源、分类及其对测量结果的影响。
- 最小二乘法:水准网平差通常采用最小二乘法进行数据处理,需要理解其数学原理和应用方法。
4. 程序设计的关键技术- 数据输入:设计高效的数据输入界面,支持多种数据格式的导入。
- 数据处理:实现数据的自动校验、筛选和处理功能。
- 平差计算:编写平差计算算法,包括闭合差计算、误差分配等。
- 结果输出:设计结果输出模块,支持多种输出格式,如文本、图表等。
5. 程序设计的实现步骤1. 需求分析:明确程序设计的目标和用户需求。
2. 系统设计:设计程序的整体架构,包括模块划分、数据流等。
3. 编码实现:根据设计文档进行编码,实现各个功能模块。
4. 测试验证:对程序进行测试,确保其准确性和稳定性。
5. 用户手册编写:编写用户手册,指导用户如何使用程序。
6. 程序设计的注意事项- 数据安全:确保程序在处理数据时的安全性,防止数据泄露。
- 异常处理:程序应能妥善处理各种异常情况,如数据格式错误、计算溢出等。
////////////////////////////////////////////////////// visual C++6.0 编译通过/////////////////////////////////////////////////////////////////////////////////////////////////////////// 参考资料//// 部分网络资料//// 宋力杰《测量平差程序设计》////姚连壁《基于matlab的控制网平差程序设计》/////////////////////////////////////////////////////#include<iostream>#include<fstream>#include <stdlib.h>#include<math.h>#include <iomanip>using namespace std;//////////////////////////////////////////////////////////////////////////c lassclass 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=ATP A[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<<"起点终点观测高差v"<<" 高差平差值观测权中误差"<<endl;for(i=0;i<=gcz_zs-1;i++){int k1=qsd_dh[i];int k2=zd_dh[i];double qii=ATP A[ij(k1,k1)];double qjj= ATPA[ij(k2,k2)] ;double qij=ATP A[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;}。