C语言版水准控制测量平差程序
测绘程序设计—实验八-水准网平差程序设计报告
《测绘程序设计()》上机实验报告(Visual C++.Net)班级: 测绘0901班学号: **********姓名: 代娅琴4月29日实验八平差程序设计基础一、实验目旳•巩固过程旳定义与调用•巩固类旳创立与使用•巩固间接平差模型及平差计算•掌握平差程序设计旳基本技巧与环节二、实验内容水准网平差程序设计。
设计一种水准网平差旳程序, 规定数据从文献中读取, 计算部分与界面无关。
水准网间接平差模型:计算示例:近似高程计算:1.水准网平差计算一般环节(1)读取观测数据和已知数据;(2)计算未知点高程近似值;(3)列高差观测值误差方程;(4)根据水准路线长度计算高差观测值旳权;(5)构成法方程;(6)解法方程, 求得未知点高程改正数及平差后高程值;(7)求高差观测值残差及平差后高差观测值;(8)精度评估;(9)输出平差成果。
2.水准网高程近似值计算算法3.输入数据格式示例实验代码:#pragma onceclass LevelControlPoint{public:LevelControlPoint(void);~LevelControlPoint(void);public:CString strName;//点名CString strID;//点号float H;bool flag;//标记与否已经计算出近似高程值, 若计算出则为, 否则为};class CDhObs{public:CDhObs(void);~CDhObs(void);public:LevelControlPoint* cpBackObj;//后视点LevelControlPoint* cpFrontObj;//前视点double ObsValue;//高差值double Dist;//测站旳距离};#include"StdAfx.h"#include"LevelControlPoint.h"LevelControlPoint::LevelControlPoint(void){strName=_T("");strID=_T("");H=0;flag=0;}LevelControlPoint::~LevelControlPoint(void){}CDhObs::CDhObs(void){}CDhObs::~CDhObs(void){}#pragma once#include"LevelControlPoint.h"#include"Matrix.h"class AdjustLevel{public:AdjustLevel(void);~AdjustLevel(void);public:LevelControlPoint* m_pKnownPoint;//已知点数组int m_iKnownPointCount;//已知点个数LevelControlPoint* m_pUnknownPoint;//未知点数组int m_iUnknownPointCount;//未知点个数CDhObs* m_pDhObs;//高差观测值数组int m_iDhObsCount;//高差观测值个数public:void SetKnownPointSize(int size);//创立大小为size旳已知点数组void SetUnkonwnPointSize(int size);//创立大小为size旳未知点数组void SetDhObsSize(int size);//创立大小为size旳观测值数组bool LoadObsData(const CString& strFile);//读入观测文献CString* SplitString(CString str, char split, int& iSubStrs);void ApproHeignt(void);//计算近似值private:LevelControlPoint* SearchKnownPointUsingID(CString ID);LevelControlPoint* SearchUnknownPointUsingID(CString ID);LevelControlPoint* SearchPointUsingID(CString ID);CMatrix LevleWeight(void);//计算权矩阵public:void FormErrorEquation(CMatrix &B, CMatrix &L);//构成误差方程void EquationCompute(CMatrix &x);//计算法方程void Accuracy_Assessment(double &r0,CMatrix &Qxx);//精度评估void CompAdjust(double &r0,CMatrix Qx[]);};#include"StdAfx.h"#include"AdjustLevel.h"#include<locale.h>#include"LevelControlPoint.h"#include"math.h"AdjustLevel::AdjustLevel(void){m_pKnownPoint=NULL;//已知点数组m_iKnownPointCount=0;//已知点个数m_pUnknownPoint=NULL;//未知点数组m_iUnknownPointCount=0;//未知点个数m_pDhObs=NULL;//高差观测值数组m_iDhObsCount=0;//高差观测值个数}AdjustLevel::~AdjustLevel(void){if(m_pKnownPoint!=NULL){delete[] m_pKnownPoint;m_pKnownPoint=NULL;}if(m_pUnknownPoint!=NULL){delete[] m_pUnknownPoint;m_pUnknownPoint=NULL;}if(m_pDhObs!=NULL){delete[] m_pDhObs;m_pDhObs=NULL;}}void AdjustLevel::SetKnownPointSize(int size){m_pKnownPoint=new LevelControlPoint[size];//创立动态指针m_iKnownPointCount=size;}void AdjustLevel::SetUnkonwnPointSize(int size){m_pUnknownPoint=new LevelControlPoint[size];m_iUnknownPointCount=size;}void AdjustLevel::SetDhObsSize(int size){m_pDhObs=new CDhObs[size];m_iDhObsCount=size;//高差观测值个数}bool AdjustLevel::LoadObsData(const CString& strFile){CStdioFile sf;if(!sf.Open(strFile,CFile::modeRead)) return false;//创立并打开文献对象CString strLine;bool bEOF=sf.ReadString(strLine);//读取第一行, 即已知点旳数目SetKnownPointSize(_ttoi(strLine));//根据已知点旳数目, 创立已知点数组;int n=0;for(int i=0;i<m_iKnownPointCount;i++)//读取已知点旳点名和高程值{sf.ReadString(strLine);CString *pstrData=SplitString(strLine,',',n);m_pKnownPoint[i].strName=pstrData[0];m_pKnownPoint[i].strID=pstrData[0];m_pKnownPoint[i].H=_tstof(pstrData[1]);m_pKnownPoint[i].flag=1;//已知点不用平差, 故将其旳flag设立为delete[] pstrData;pstrData=NULL;}sf.ReadString(strLine);//读取未知点旳个数SetUnkonwnPointSize(_ttoi(strLine));//根据未知点旳个数创立未知点数组sf.ReadString(strLine);//读取未知点旳点名CString *pstrData=SplitString(strLine,',',n);for(int i=0;i<m_iUnknownPointCount;i++)//将未知点旳点名放入未知点数组{m_pUnknownPoint[i].strName=pstrData[i];m_pUnknownPoint[i].strID=pstrData[i];m_pUnknownPoint[i].H=0;//未知点旳高程值设立为m_pUnknownPoint[i].flag=0;//还没有求得近似高程, 故其flag设立为}if(pstrData!=NULL){delete[] pstrData;pstrData=NULL;}sf.ReadString(strLine);//读取观测值旳个数SetDhObsSize(_ttoi(strLine));//按照观测值旳大小, 创立观测值数组for(int i=0;i<m_iDhObsCount;i++)//分行读取观测值旳数据, 将其存入观测值数组{sf.ReadString(strLine);CString *pstrData=SplitString(strLine,',',n);m_pDhObs[i].cpBackObj=SearchPointUsingID(pstrData[0]);//后视点m_pDhObs[i].cpFrontObj=SearchPointUsingID(pstrData[1]);//前视点m_pDhObs[i].HObsValue=_tstof(pstrData[2]);//高差观测值m_pDhObs[i].Dist=_tstof(pstrData[3]);//距离观测值delete[] pstrData;pstrData=NULL;}sf.Close();return 1;}CString* AdjustLevel::SplitString(CString str, char split, int& iSubStrs) {int iPos = 0; //分割符位置int iNums = 0; //分割符旳总数CString strTemp = str;CString strRight;//先计算子字符串旳数量while (iPos != -1){iPos = strTemp.Find(split);if (iPos == -1){break;}strRight = strTemp.Mid(iPos + 1, str.GetLength());strTemp = strRight;iNums++;}if (iNums == 0) //没有找到分割符{//子字符串数就是字符串自身iSubStrs = 1;return NULL;}//子字符串数组iSubStrs = iNums + 1; //子串旳数量= 分割符数量+ 1CString* pStrSplit;pStrSplit = new CString[iSubStrs];strTemp = str;CString strLeft;for (int i = 0; i < iNums; i++){iPos = strTemp.Find(split);//左子串strLeft = strTemp.Left(iPos);//右子串strRight = strTemp.Mid(iPos + 1, strTemp.GetLength());strTemp = strRight;pStrSplit[i] = strLeft;}pStrSplit[iNums] = strTemp;return pStrSplit;}//LevelControlPoint* AdjustLevel::SearchKnownPointUsingID(CString ID) {for(int i=0;i<m_iKnownPointCount;i++){if(m_pKnownPoint[i].strID==ID){return &m_pKnownPoint[i];}}return NULL;}//LevelControlPoint* AdjustLevel::SearchUnknownPointUsingID(CString ID) {for(int i=0;i<m_iUnknownPointCount;i++){if(m_pUnknownPoint[i].strID==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)//用于计算高程近似值旳函数{for(int i=0;i<m_iUnknownPointCount;i++)//计算未知点高程值{if(m_pUnknownPoint[i].flag!=1){//先在未知点作为观测值旳前视点旳状况for(int j=0;j<m_iDhObsCount;j++)//从观测数组里找与未知点有关联旳点{//如果观测值旳前视点是未知点且其后视点已有高程值if((m_pDhObs[j].cpFrontObj->strID==m_pUnknownPoint[i].strID)&& m_pDhObs[j].cpBackObj->flag==1 ){ //前视点=后视点-高差/*m_pUnknownPoint[i].H=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)//如果通过上一环节未知点仍没有计算出近似值{for(int j=0;j<m_iDhObsCount;j++)//从观测数组里找与未知点有关联旳点 {//如果观测值旳后视点是未知点且其前视点已有高程值if((m_pDhObs[j].cpBackObj->strID==m_pUnknownPoint[i].strID)&& m_pDhObs[j].cpFrontObj->flag==1 ){ //后视点=前视点+高差m_pUnknownPoint[i].H=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)//如果已经计算到最后一种未知点{for(int a=0;a<m_iUnknownPointCount;a++){if(m_pUnknownPoint[i].flag!=1)//只要有一种未知点旳近似高程直没有计算{ //则要重新进行上面旳环节直到所有旳未知点旳近似高程值都计算出i=-1;break;}}}}}CMatrix AdjustLevel::LevleWeight(void){CMatrix p(m_iDhObsCount,m_iDhObsCount);p.Unit();double value;for(int i=0;i<m_iDhObsCount;i++){value=(1.0/m_pDhObs[i].Dist);p(i,i)=value;}return p;}void AdjustLevel::FormErrorEquation(CMatrix &B, CMatrix &L){B.SetSize(m_iDhObsCount,m_iUnknownPointCount);L.SetSize(m_iDhObsCount,1);for(int i=0;i<m_iDhObsCount;i++)//建立B系数阵{LevelControlPoint *tmpBack=NULL,*tmpFront=NULL;tmpBack=SearchPointUsingID(m_pDhObs[i].cpBackObj->strID);tmpFront=SearchPointUsingID(m_pDhObs[i].cpFrontObj->strID);//找到与第i个观测值有关旳未知点tmpBack->strID;for(int j=0;j<m_iUnknownPointCount;j++){if(m_pUnknownPoint[j].strID==tmpBack->strID)//如果是后视点则前面旳系数为-1{ B(i,j)=-1;continue;}if(m_pUnknownPoint[j].strID==tmpFront->strID)//如果是前视点则前面旳系数为{B(i,j)=1;}}}//建立L矩阵CString tmp;for(int i=0;i<m_iDhObsCount;i++){//l=高差观测值-(后视近似值-前视近似值)/*L(i,0)=m_pDhObs[i].ObsValue-(m_pDhObs[i].cpBackObj->H-m_pDhObs[i].cpFrontObj->H);*/ L(i,0)=m_pDhObs[i].HObsValue-(m_pDhObs[i].cpFrontObj->H -m_pDhObs[i].cpBackObj->H);tmp.Format(_T("%.3f"),L(i,0));L(i,0)=_tstof(tmp);L(i,0)=L(i,0)*1000;//将单位化为mm}}void AdjustLevel::EquationCompute(CMatrix &x)//计算法方程{CMatrix P,B,l;P=LevleWeight(); //P为权矩阵FormErrorEquation(B,l);ApproHeignt();CMatrix BT(m_iUnknownPointCount,m_iDhObsCount);BT=~B; //B旳转置矩阵CMatrix NBB(m_iUnknownPointCount,m_iUnknownPointCount);NBB=BT*P*B;CMatrix NBBl=NBB.Inv();x=NBBl*BT*P*l;for(int i=0;i<m_iUnknownPointCount;i++){m_pUnknownPoint[i].H+=x(i,0);//未知点高程值=近似值+改正数}}void AdjustLevel::Accuracy_Assessment(double &r0,CMatrix &Qxx)//精度评估{CMatrix B,l,P,x;P=LevleWeight(); //P为权矩阵FormErrorEquation(B,l);EquationCompute(x);CMatrix v(m_iDhObsCount,1);v=B*x-l;CMatrix vT(1,m_iDhObsCount);vT=~v;CMatrix r/*(1,l)*/;r=vT*P*v;r0=sqrt(r(0,0)/(m_iDhObsCount-m_iUnknownPointCount));//单位权中误差Qxx.SetSize(m_iUnknownPointCount,m_iUnknownPointCount);CMatrix BT(m_iUnknownPointCount,m_iDhObsCount);BT=~B;CMatrix NBB(m_iUnknownPointCount,m_iUnknownPointCount);NBB=BT*P*B;Qxx=NBB.Inv();}void AdjustLevel::CompAdjust(double &r0,CMatrix Qx[]){ApproHeignt();//计算未知点旳近似高程值并且存入数组CMatrix P(m_iDhObsCount,m_iDhObsCount);P=LevleWeight();//p为权矩阵CMatrix B,L;CMatrix x,Qxx;FormErrorEquation(B,L);//构成误差方程, B为系数矩阵, l为常数项EquationCompute(x);//计算法方程Accuracy_Assessment(r0,Qxx);//精度评估for(int i=0;i<m_iUnknownPointCount;i++)//未知点高程中误差{Qx[i]=sqrt(Qxx(i,i))*r0;}}#include"Matrix.h"#include"locale.h"#include"LevelControlPoint.h"#include"AdjustLevel.h"AdjustLevel LevelComput;CString* SplitString(CString str, char split, int& iSubStrs){int iPos = 0; //分割符位置int iNums = 0; //分割符旳总数CString strTemp = str;CString strRight;//先计算子字符串旳数量while (iPos != -1){iPos = strTemp.Find(split);if (iPos == -1){break;}strRight = strTemp.Mid(iPos + 1, str.GetLength());strTemp = strRight;iNums++;}if (iNums == 0) //没有找到分割符{//子字符串数就是字符串自身iSubStrs = 1;return NULL;}//子字符串数组iSubStrs = iNums + 1; //子串旳数量= 分割符数量+ 1CString* pStrSplit;pStrSplit = new CString[iSubStrs];strTemp = str;CString strLeft;for (int i = 0; i < iNums; i++){iPos = strTemp.Find(split);//左子串strLeft = strTemp.Left(iPos);//右子串strRight = strTemp.Mid(iPos + 1, strTemp.GetLength());strTemp = strRight;pStrSplit[i] = strLeft;}pStrSplit[iNums] = strTemp;return pStrSplit;}void CIndircLelveDlg::OnBnClickedOpendatafile(){// TODO: 在此添加控件告知解决程序代码UpdateData(TRUE);CFileDialog dlgFile(TRUE,_T("txt"),NULL,OFN_ALLOWMULTISELECT|OFN_EXPLORER,_T("(文本文献)|*.txt"));//创立文献对话框if(dlgFile.DoModal()==IDCANCEL) return;//如果选择取消按钮则返回CString strFileName=dlgFile.GetPathName();//打开获取文献文献名setlocale(LC_ALL,""); //设立语言环境CStdioFile sf;if(!sf.Open(strFileName, CFile::modeRead)) return;InputContent.Empty();//清空字符串str_openContent中旳内容CString strLine;BOOL bEOF=sf.ReadString(strLine);//读取第一行数据while(bEOF)//开始读取顶点数据{bEOF=sf.ReadString(strLine);if(bEOF)InputContent+=strLine+_T("\r\n");}sf.Close();UpdateData(FALSE);}void CIndircLelveDlg::OnBnClickedSavedata(){// TODO: 在此添加控件告知解决程序代码U pdateData(TRUE);CFileDialog dlgFile(FALSE,_T("txt"),NULL,OFN_EXPLORER,_T("(Level格式)|*.txt"));if(dlgFile.DoModal()==IDCANCEL) return;CString strFileName=dlgFile.GetPathName();setlocale(LC_ALL,"");CStdioFile sf;if(!sf.Open(strFileName, CFile::modeCreate|CFile::modeWrite)) return;sf.WriteString(LevleContent);sf.Close();UpdateData(FALSE);}void CIndircLelveDlg::OnBnClickedComputelevel(){// TODO: 在此添加控件告知解决程序代码UpdateData(TRUE);setlocale(LC_ALL,"");double *Qx=new double[LevelComput.m_iUnknownPointCount];double r0;pAdjust(r0,Qx);LevleContent.Format(_T("平差后高程值:\r\n"));CString Temp;for(int i=0;i<LevelComput.m_iUnknownPointCount;i++){Temp.Empty();Temp.Format(_T("%s,%.4f\r\n"),LevelComput.m_pUnknownPoint[i].strID,LevelComput.m_pUnknownPoint[i].H);LevleContent+=Temp;}Temp.Format(_T("单位权中误差: %.1f mm\r\n"),r0*1000);LevleContent+=Temp;LevleContent+=_T("未知点高程中误差(mm):\r\n");for(int i=0;i< LevelComput.m_iUnknownPointCount;i++){Temp.Empty();Temp.Format(_T("%s,%.1f\r\n"),LevelComput.m_pUnknownPoint[i].strName,Qx[i]*1000);LevleContent+=Temp;}UpdateData(false);}void CIndircLelveDlg::OnBnClickedSavelevleresult(){// TODO: 在此添加控件告知解决程序代码UpdateData(TRUE);CFileDialog dlgFile(FALSE,_T("txt"),NULL,OFN_EXPLORER,_T("(Level格式)|*.txt"));if(dlgFile.DoModal()==IDCANCEL) return;CString strFileName=dlgFile.GetPathName();setlocale(LC_ALL,"");CStdioFile sf;if(!sf.Open(strFileName, CFile::modeCreate|CFile::modeWrite)) return;sf.WriteString(LevleContent);sf.Close();UpdateData(FALSE);}三、实验成果打开文献数据:平差成果:四、实验心得这从实验是我们测绘程序设计旳最后一次实验, 虽然这个学期我们做了好几次有关旳实验, 但是我却发现自己学旳东西也越来越模糊, 感觉诸多内容都不理解。
水准网间接平差程序设计(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<<"起点终点观测高差v"<<" 高差平差值观测权中误差"<<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;}。
C语言实现测量数据误差处理
误差处理程序(C语言)电子信息工程学院通信100910211159高子豪实验目的实现对输入数据的误差处理:剔除粗大误差。
判断累进性系统误差和周期性系统误差。
计算平均值,方差,不确定度。
程序代码#include<stdio.h>#include<math.h>double SUM(double x[],int n);double AVRG(double x[],int n);double SD(double x[],int n);int PauTa(double x[],int n);int Chauvenet(double x[],int n);int Grubbs_1(double x[],int n);int Grubbs_2(double x[],int n);static int n;static double a[500];int main(){int i,choose,leap=1;double avg,sd,v[500],M=0,AH=0,vmax=0;doubleP,PX[]={12.706,4.303,3.182,2.776,2.571,2.447,2.365,2.306,2.262,2.228,2.131,2.086,2.060,2.042,2.021,2 .000,1.980,1.960};printf("请输入数据总个数:\n");scanf("%d",&n);printf("请输入数据:\n");for(i=0;i<n;i++)scanf("%lf",&a[i]);avg=AVRG(a,n);sd=SD(a,n);printf("\n输入数据的平均值为%lf,标准差为%lf\n",avg,sd);while(leap){printf("请选择粗大误差的检验法:\n1.莱特检验法\n2.肖维纳检验法\n3.格拉布斯检验法(置信概率99%%)\n4.格拉布斯检验法(置信概率95%%)\n5.停止检验\n");scanf("%d",&choose);if(choose==1&&n<10)printf("数据总量小于10,不能使用莱特检验法。
测量平差程序设计6
测量平差程序设计
n1 = 2 * (zds - yds) '未知数数目 n2 = n1 * (n1 + 1) / 2 '一维存储法方程系数 数组上限 ReDim NX(n2), UX(n1) ' 重新定义法方程 系数、常数数组 Call order(m(), yds) '对保存已知点序号的 m()数组排序
测量平差程序设计
If cha = "n" Then '照准方向点不是已知点 d = 2 * (h - k - 1) + 1 '计算照准方向点 x坐标未知数在未知数点集中的序号 nb(d) = -A: nb(d + 1) = -B '对误差方程 系数数组赋值 nc(d) = nc(d) + nb(d): nc(d + 1) = nc(d + 1) + nb(d + 1) '累积到和方程数组中 去 End If
测量平差程序设计
Private Function azimuth(x1 As Double, y1 As Double, x2 As Double, y2 As Double) '反算坐标方位角 dx = x2 - x1: dy = y2 - y1 azimuth = Atn(dy / dx) If dx < 0 Then azimuth = azimuth + pi Else If dy < 0 Then azimuth = azimuth + 2 * pi End If End Function
测量平差程序设计
ReDim nc(n1): ln = 0 ' 重新定义和方程系数数组,并且 每循环到一新测站前清零,ln是和方程常数项 For j = k1 To k2 '在i测站上按方向循环,求误差方程系数、 常数 ReDim nb(n1) ' 重新定义误差方程系数数组,并且每循 环到一新方向前清零 h = seqn(lb(j)) dx = x(h) - x(i): dy = y(h) - y(i) ss = Sqr(dx ^ 2 + dy ^ 2) '反算边长 t = azimuth(x(i), y(i), x(h), y(h)) '反算坐标方位角
水准网平差程序
AfxGetApp()->m_pMainWnd->MessageBox( "数据文件不存在或数据文件错!", "进程. . . . . .!!!",MB_OK|MB_ICONSTOP);
//重要说明:原始数据文件中,未知点的高程可以随意输入,也可以不输入空缺, 程序自动把待定点高程赋值为 0
} for(i=ne;i<nz;i++) {
fp.ReadString(buff,MAXLINE); sscanf(buff,"%d%s%lf",&dh,ch1,&gc); dm[i]=ch1;H[i]=0; } for(i=0;i<nn;i++) { fp.ReadString(buff,MAXLINE); sscanf(buff,"%d%d%lf%lf",
教师 评语
// Gckzwpc.cpp: implementation of the CGckzwpc class. // ////////////////////////////////////////////////////////////////////// #include "stdafx.h" #include "Survey.h" #include "Gckzwpc.h" #ifdef _DEBUG #undef THIS_FILE static char THIS_FILE[]=__FILE__; #define new DEBUG_NEW #endif ////////////////////////////////////////////////////////////////////// // Construction/Destruction ////////////////////////////////////////////////////////////////////// CGckzwpc::CGckzwpc() { } CGckzwpc::~CGckzwpc() { } bool CGckzwpc::ReadData(CString filename) {
C语言测量平差课程设计
C语言测量平差课程设计一、课程目标知识目标:1. 让学生掌握C语言中常用的数据类型、运算符和语法结构,并能将其应用于测量平差计算中。
2. 使学生了解测量平差的基本原理和常用算法,如最小二乘法等。
3. 帮助学生理解C语言在解决测量平差问题中的优势,如计算速度和精度。
技能目标:1. 培养学生运用C语言进行测量平差计算的能力,能独立编写、调试和优化相关程序。
2. 提高学生分析测量数据、选择合适算法解决问题的能力。
3. 培养学生团队合作精神和沟通能力,能共同完成测量平差项目。
情感态度价值观目标:1. 激发学生对测量平差和C语言的兴趣,培养主动学习和探究的精神。
2. 引导学生关注测量平差在工程、科研等领域的应用,认识到所学知识的社会价值。
3. 培养学生严谨、求实的科学态度,遵循学术道德,尊重他人成果。
本课程针对高年级学生,课程性质为理论与实践相结合。
在分析课程性质、学生特点和教学要求的基础上,将课程目标分解为具体的学习成果,为后续的教学设计和评估提供依据。
通过本课程的学习,学生将能够运用所学知识解决实际问题,提高C语言编程能力和测量平差理论水平。
二、教学内容1. C语言基础回顾:数据类型、变量、运算符、控制结构、函数等基本概念,重点复习数组、指针和结构体。
教材章节:第一章至第三章2. 测量平差基本原理:介绍测量平差的定义、目的、数学模型和常用算法(如最小二乘法)。
教材章节:第四章3. C语言实现测量平差算法:结合实际案例,讲解如何使用C语言编写测量平差程序。
教材章节:第五章4. 程序设计与调试:介绍编程规范,演示调试技巧,分析常见错误。
教材章节:第六章5. 测量平差实际应用:分析实际测量数据,运用所学知识解决具体问题。
教材章节:第七章6. 项目实践:分组进行项目设计,完成测量平差程序编写、调试和优化,撰写项目报告。
教材章节:第八章教学内容安排和进度:1. 第1周:C语言基础回顾。
2. 第2周:测量平差基本原理。
4.5水准网平差程序
本程序用于二、三、四等水准网平差计算和精度评定。采用间接观 测平差法,按结点平差法解算,以归算后的观测高差为平差元素, 计算待定点的高程平差值,计算单位权中误差(1公里)和每个待 定点的高程平差值的中误差。
3.编号说明 点的编号,先编已知点,接着编待定点,从1开始顺序编号。 观测高差(或路线)的编号,从1开始顺序编号。为推求高程的方 便,需从与某已知点相连的观测高差边开始编号,并照顾高程的传 递。观测高差方向可任意。 4.关于输入数据数据库的建立和计算成果输出的说明 将全部输入数据建立成Access数据库的三张表,并约定表名分 别为:table1、table2和table3,这里表名是主要的。
测量平差程序设计5 (2)
导线网近似坐标解算 p=0 For k = 1 To zds If x(k) < 1 Then p = 1 '查看是否还有没解算出坐标的点,有则对p赋值,使进入 查看是否还有没解算出坐标的点,有则对p
循环再次搜索计算
Next k Loop Until p = 0
导线网近似坐标解算 For i = 1 To yds
'查点名数组中那些点是已知点,用m(i)数组存其序号 查点名数组中那些点是已知点, m(i)数组存其序号
For j = 1 To zds If ym(i) = dm(j) Then m(i) = j Next j Next I
‘循环结束后,m(i)中按已知点在ym()中出现的的顺序,保 循环结束后,m(i)中按已知点在ym()中出现的的顺序,保 存着已知点在dm()数组中的序号,即ym(1)、ym(2)… 存着已知点在dm()数组中的序号,即ym(1)、ym(2)…, 在dm()中的序号。 dm()中的序号。
参考题 1、简述导线网近似坐标推算的算法? 2、什么情况下,两点之间查不到观测边? 3、坐标转换计算时,为什么执行语句: dx = x(i) - x1: dy = y(i) - y1,用 y1, dx、 dy代替旋转变换公式中的x(i) 、 dx、 dy代替旋转变换公式中的x(i) y(i) ? 4、坐标变换公式中,平移参数是什么?
环进行
For i = 1 To cds '按测站循环 If x(i) > 1 Then '在该测站假设坐标已计算出的似坐标解算 For j = nl(i - 1) + 1 To nl(i) ‘遍访i测站上所有方向 遍访i
值,查找一个已解算出的照准点,反算方位角。
测绘程序设计—实验八 水准网平差程序设计报告
《测绘程序设计》上机实验报告(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);}三、实验结果打开文件数据:平差结果:四、实验心得这从实验是我们测绘程序设计的最后一次实验,虽然这个学期我们做了好几次相关的实验,但是我却发现自己学的东西也越来越模糊,感觉很多内容都不理解。
测量平差程序设计4 (2)
读入观测数据
ts.Close cds = j: yds = k '用模块 '用模块 级变量cds 级变量cds 、 yds保存测站点总数、 yds保存测站点总数、 已知点总数 MsgBox "数据已成功读入", "数据已成功读入", 0 + 64 + 0, "数据输入" "数据输入" End If
读入观测数据
tr(1) = Mid(B, 1, m(1) - 1) '提取指定位置开始的指定数目字 符。 m(2) = InStr(m(1) + 1, B, ",") '从上一个找到的逗号位置起, '从上一个找到的逗号位置起, 查找下一个逗号的位置 tr(2) = Mid(B, m(1)+1,m(2)m(1)+1,m(2)M(1)M(1)-1)
Else j = j + 1:ReDim Preserve dm(j): ReDim Preserve nl(j): ReDim Preserve ns(j): dm(j) = B: nl(j) = nl(j - 1): ns(j) = ns(j 1):p=2 End if `m(1)=0,表明读到测站了,以后遇到 `m(1)=0,表明读到测站了,以后遇到 有逗号的行,只能是观测值了。所以 将识别符p赋值2 将识别符p赋值2,目的是使得提取已 知点信息的语句不再执行。
读入观测数据
B = ts.ReadLine ‘读数据文件 中一行,置入变量B 中一行,置入变量B B = Trim(B)'删除B可能有的前 Trim(B)'删除B 导和尾随空格
读入观测数据
实验一 利用平差易软件和程序设计语言完成水准网平差
实验一利用平差易软件和程序设计语言完成水准网平差一、实验目的1.掌握水准平差的基本原理。
2.掌握程序语言设计的基础。
3. 掌握平差易软件的基本功能。
4. 能够独立完成水准网观测数据的平差处理过程二、实验数据如图1所示水准网,有1个已知点,4个未知点,8个测段。
各已知数据及观测值见下表1,已知点高程H1=98m 。
图1表1.高差观测值(m)(3)求各待定点的高程;2号点、3号点、4号点的高程中误差。
三、控制网平差报告[控制网概况]1.本成果为按[平面]网处理的平差成果计算软件:南方平差易2002网名:xj 计算日期:2016/6/29 13:36观测人:xujiang 记录人:xujaing计算者:xujaing 测量单位:m备注:20144159高程控制网等级:等外水准2.控制网中:方向方位平距XY点待定8 5固定每公里高差中误差= 20.11 (mm)[距离观测成果表]测站照准距离(m) 改正数(m) 平差后值(m) 方位角(dms) 1 5 1380.0000 0.0000 1380.0000 360.0000003 1420.0000 0.0000 1420.0000 360.0000002 1810.0000 0.0000 1810.0000 360.0000002 3 940.0000 0.0000 940.0000 360.0000004 1760.0000 0.0000 1760.0000 360.0000003 4 1400.0000 0.0000 1400.0000 360.0000005 990.0000 0.0000 990.0000 360.000000 4 5 1350.0000 0.0000 1350.0000 360.000000[高差观测成果表]测站照准高差(m) 改正数(m) 平差后值(m) 备注1 5 3.1020 -0.0174 3.08463 3.5200 0.0335 3.55352 2.5420 -0.0199 2.52212 3 1.0340 -0.0025 1.03154 -1.5540 -0.0147 -1.56873 4 -2.6110 0.0108 -2.60025 -0.4820 0.0131 -0.46894 5 2.1320 -0.0008 2.1312[平面点位误差表]点名长轴(m) 短轴(m) 长轴方位dms 点位中误差m高程中误差m2 0.01813 0.01614 0.02015 0.0171[控制点成果表]点名X(m) Y(m) H(m) 备注1 0.0000 0.0000 98.0000 固定点2 0.0000 0.0000 100.54203 0.0000 0.0000 101.52004 0.0000 0.0000 98.98805 0.0000 0.0000 101.1020。
C语言实现测量数据误差处理
C语言实现测量数据误差处理数据误差处理是在实际测量中必须考虑的一个重要问题。
在科学研究和工程控制系统中,测量数据的准确性对结果分析和决策具有重要意义。
因此,在数据处理过程中,需要考虑数据误差的影响,并进行合理的处理和分析。
在C语言中,可以通过一些常用的方法和技巧来实现数据误差的处理。
下面将介绍几种常见的数据误差处理方法。
1.数据滤波数据滤波是一种常用的数据误差处理方法。
它通过对一系列相邻数据进行加权平均或滑动平均等操作,来降低误差的影响。
在C语言中,可以使用数组来存储测量数据,并通过循环的方式来实现滤波计算。
例如,可以定义一个数组来存储测量数据:```c#define DATA_SIZE 100float data[DATA_SIZE]; // 存储测量数据的数组```然后,可以通过循环为数组赋值,并进行滤波计算:```c#define FILTER_WINDOW_SIZE 5void smoothing_filter(float *input, int input_size, float *output)int i, j;for (i = FILTER_WINDOW_SIZE / 2; i < input_size -FILTER_WINDOW_SIZE / 2; i++)output[i] = 0;for (j = -FILTER_WINDOW_SIZE / 2; j <= FILTER_WINDOW_SIZE / 2; j++)output[i] += input[i + j];}output[i] /= FILTER_WINDOW_SIZE;}```2.数据整定数据整定是一种通过对测量数据进行修正来降低误差的方法。
它主要是通过校正系数或拟合公式等方式,将测量数据修正到一个更准确的值。
在C语言中,可以使用数学库函数来实现数据整定。
例如,可以使用线性拟合的方法来修正数据:```c#include <math.h>//线性拟合void linear_calibration(float *input, int input_size, float *output, float a, float b)int i;for (i = 0; i < input_size; i++)output[i] = a * input[i] + b;}```3.数据统计和分析数据统计和分析是一种通过对测量数据进行统计和分析,来确定误差范围和误差概率的方法。
水准平差程序说明
附、闭合水准平差程序说明一、概述本程序可完成各等级附、闭合、水准路线、水准支线的平差计算,加入尺长改正、正高改正对单一附闭合水准路线进行严密平差。
打印输出观测数据、线路闭合差、各点高程、每公里高差中误差。
二、程序概貌(1)、数据输入菜单实现观测数据文件的建立和数据文件路径的设定。
(2)、平差计算菜单完成平差计算。
(3)、成果输出菜单输出平差计算成果。
(4)、帮助菜单数据文件格式说明。
(5)、结束菜单结束程序。
三、操作方法(1)、建立观测数据文件本程序通过调用WINDOWS下记事本程序建立观测数据文件,也可直接用其它文本编辑软件建立观测数据文件。
选择数据文件菜单下的新建菜单。
数据文件格式为:测段数,尺长改正数,等级,是否支线(支线输1,否则输0)起点点名,纬度,起点高程,终点高程下一点点名,纬度,测段距离,往测测站数,返测测站数,往测高差,返测高差………………直到结束。
如本线路无返测高差则返测高差输0。
例:4,0.0003,4,0135j,24.28, 1295.2741,1295.27412yl1,24.25,1.2,12,14,20.34442,-20.346282yl2,24.22,1.0,10,10,77.30418,-77.302852yl3,24.19,0.8,8,8,55.57608,-55.57765135j,24.28,3.0,30,30,-153.2186,153.2196第一行起点名为135j,纬度为24°28′,起点高程1295.2741,终点高程为1295.2741。
第二行点名为2y11,纬度为24°25′,1.2为测段距离,12为测段的往测站数,14为测段的返测站数,20.34442为往测高差,-20.34628为返测高差。
………………第五行终点名为135j,纬度为24°28′,3.0为测段距离,30为测段的往测站数,30为测段的返测站数,-153.2186为往测高差,153.2196为返测高差。
windows绘图和水准网平差c++代码
Windows绘图#include "stdafx.h"#include "cbs.h"#include "cbsDoc.h"#include "cbsView.h"#ifdef _DEBUG#define new DEBUG_NEW#undef THIS_FILEstatic char THIS_FILE[] = __FILE__;#endif/////////////////////////////////////////// //////////////////////////////////// CCbsViewIMPLEMENT_DYNCREATE(CCbsView, CView)BEGIN_MESSAGE_MAP(CCbsView, CView)//{{AFX_MSG_MAP(CCbsView)ON_WM_MOUSEMOVE()ON_WM_LBUTTONDOWN()ON_WM_LBUTTONUP()//}}AFX_MSG_MAP// Standard printing commandsON_COMMAND(ID_FILE_PRINT,CView::OnFilePrint)ON_COMMAND(ID_FILE_PRINT_DIRECT, CView::OnFilePrint)ON_COMMAND(ID_FILE_PRINT_PREVIEW, CView::OnFilePrintPreview)END_MESSAGE_MAP()/////////////////////////////////////////// //////////////////////////////////// CCbsView construction/destruction CCbsView::CCbsView(){// TODO: add construction code herethis->m_bDragging=false;}CCbsView::~CCbsView(){}BOOL CCbsView::PreCreateWindow(CREATESTRUCT& cs){// TODO: Modify the Window class or styles here by modifying// the CREATESTRUCT csreturn CView::PreCreateWindow(cs);}/////////////////////////////////////////// //////////////////////////////////// CCbsView drawingvoid CCbsView::OnDraw(CDC* pDC){CCbsDoc* pDoc = GetDocument();ASSERT_VALID(pDoc);// TODO: add draw code for native data here /* pDC->SetMapMode(MM_TEXT);pDC->Rectangle(CRect(50, 50, 100, 100));// 直接基于屏幕绘制pDC->SetMapMode(MM_TEXT);//pDC->SetWindowOrg(50, 50); //屏幕左上角的坐标设置为(50, 50)pDC->Rectangle(CRect(50, 50, 100, 100));pDC->SetMapMode(MM_TEXT);pDC->SetViewportOrg(50,50);//当前原点位置移动到(50, 50)的位置pDC->Rectangle(CRect(50, 50, 100, 100));pDC->SetMapMode(MM_TEXT);pDC->SetViewportOrg(100,100); // 当前原点位置移动到(50, 50)的位置pDC->SetWindowOrg(50, 50); // 当前的原点坐标设置为(50, 50)pDC->Rectangle(CRect(100, 100, 100, 100));pDC->SetWindowOrg(100,100);*/pDC->SetViewportOrg(150,150); // 当前原点位置移动到(50, 50)的位置COLORREF rgbBkClr = RGB(192,192,192);// 定义灰色pDC->SetBkColor(rgbBkClr);// 背景色为灰色pDC->SetTextColor(RGB(0,0,128));// 文本颜色为兰色pDC->TextOut(250,10,"我*****的*****绘*****图*****板!");CPen *pPenOld, PenNew;int nPenStyle[]= { PS_SOLID, // 实线PS_DOT,// 点线PS_DASH,// 虚线PS_DASHDOT,// 点划线PS_DASHDOTDOT, // 双点划线PS_NULL,// 空的边框PS_INSIDEFRAME,// 边框实线};char *strStyle[]={"Solid", "Dot", "Dash", "DashDot","DashDotDot", "Null", "InsideFrame"};pDC->TextOut(50, 65, "用不同样式的画笔绘图"); pDC->TextOut(400, 40, "***制作者***陈斌生***");for(int i=0; i<7; i++){ // 用不同样式的画笔绘图if(PenNew.CreatePen(nPenStyle[i],1,RGB(0,0, 0))){ //创建画笔pPenOld=pDC->SelectObject(&PenNew);// 选择画笔pDC->TextOut(45,90+20*i,strStyle[i]);pDC->MoveTo(150,100+20*i);pDC->LineTo(200,100+20*i);pDC->SelectObject(pPenOld);// 恢复原来的画笔PenNew.DeleteObject();// 删除底层的GDI对象}else {MessageBox("不能创建画笔!"); } }char *strWidth[]={"1", "2", "3", "4", "5", "6", "7"};pDC->TextOut(260, 65, "用不同宽度的画笔绘图");for(i=0; i<7; i++){ // 用不同宽度的画笔绘图if(PenNew.CreatePen(PS_SOLID,i+1,RGB(0, 0,0))){ // 创建画笔pPenOld=pDC->SelectObject(&PenNew); // 选择画笔pDC->TextOut(260,90+20*i,strWidth[i]);pDC->MoveTo(300, 100+20*i);pDC->LineTo(400, 100+20*i);pDC->SelectObject(pPenOld);// 恢复原来的画笔PenNew.DeleteObject();// 删除底层的GDI对象}else { MessageBox("不能创建画笔!"); }}char *strColor[]={"红","绿","蓝","黄","紫","青","灰"};COLORREF rgbPenClr[]={RGB(255,0,0), RGB(0,255,0),RGB(0,0,255), RGB(255,255,0), RGB(255,0,255),RGB(0,255,255),RGB(192,192,192)};pDC->TextOut(460,65,"用不同颜色的画笔绘图");for(i=0; i<7; i++){ // 用不同颜色的画笔绘图CPen *pPenNew=new CPen(PS_INSIDEFRAME,10,rgbPenClr[i]);// 创建画笔的另一种方法pPenOld=pDC->SelectObject(pPenNew);// 选择创建的画笔pDC->TextOut(460,90+20*i, strColor[i]);pDC->MoveTo(500,100+20*i);pDC->LineTo(600,100+20*i);pDC->SelectObject(pPenOld); // 恢复原来的画笔delete pPenNew; // 自动删除底层的GDI对象}}/////////////////////////////////////////// //////////////////////////////////// CCbsView printingBOOL CCbsView::OnPreparePrinting(CPrintInfo* pInfo){// default preparationreturn DoPreparePrinting(pInfo);} void CCbsView::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/){// TODO: add extra initialization before printing}void CCbsView::OnEndPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/){// TODO: add cleanup after printing}/////////////////////////////////////////// //////////////////////////////////// CCbsView diagnostics#ifdef _DEBUGvoid CCbsView::AssertValid() const{CView::AssertValid();}void CCbsView::Dump(CDumpContext& dc) const {CView::Dump(dc);}CCbsDoc* CCbsView::GetDocument() // non-debug version is inline{ASSERT(m_pDocument->IsKindOf(RUNTIME_CL ASS(CCbsDoc)));return (CCbsDoc*)m_pDocument;}#endif //_DEBUG/////////////////////////////////////////// //////////////////////////////////// CCbsView message handlersvoid CCbsView::OnMouseMove(UINT nFlags, CPoint point)// TODO: Add your message handler code here and/or call defaultif(m_bDragging){CClientDC dc(this);OnPrepareDC(&dc);dc.DPtoLP(&point);dc.MoveTo(this->m_ptOrigin);dc.LineTo(point);m_ptOrigin = point;}CView::OnMouseMove(nFlags, point);}void CCbsView::OnLButtonDown(UINT nFlags, CPoint point){// TODO: Add your message handler code here and/or call defaultCClientDC dc(this);OnPrepareDC(&dc); // 调整设备环境的属性dc.DPtoLP(&point); // 将设备坐标转换为逻辑坐标SetCapture(); // 捕捉鼠标//::SetCursor(m_hCross); // 设置十字光标m_ptOrigin=point;m_bDragging=TRUE; // 设置拖拽标记CView::OnLButtonDown(nFlags, point);}void CCbsView::OnLButtonUp(UINT nFlags, CPoint point){// TODO: Add your message handler code here and/or call defaultm_bDragging = false;ReleaseCapture();CView::OnLButtonUp(nFlags, point);}结果水准网平差结果#include<iostream.h>#include<stdlib.h>#include<iomanip.h>#include<math.h>#define max 50class CMatrix{public:CMatrix(){row=0; column=0;}; // 默认构造函数CMatrix(int i, int j){row=i;column=j;} // 构造函数一CMatrix(const CMatrix& m); // 复制构造函数~CMatrix(void){/*cout<<"谢谢使用,矩阵所占空间以释放!"<<endl;*/} // 默认析构函数CMatrix& operator=(const CMatrix& m); // 赋值运算符bool operator==(const CMatrix& m); // 比括较运算符bool operator!=(const CMatrix& m); // 比括较运算符CMatrix operator+(const CMatrix& m); // 加运算符CMatrix operator-(const CMatrix& m); // 减运算符CMatrix& operator+=(const CMatrix& m); // 自加运算符CMatrix& operator-=(const CMatrix& m); // 自减运算符CMatrix operator-();// 取负数CMatrix& operator*(const CMatrix& 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]; }// 设置行、列的值void setRow(const int i) { row = i; }int getRow() const { return row; }void setColunm(const int i) { column = i; }int getColumn() const { return column; }CMatrix& change(int i, int j);//交换矩阵的行CMatrix& transpose(); // 矩阵转置CMatrix& inverse(); // 矩阵求逆void find(int& f)const;// 判断该矩阵是否可用于迭代求解friend void jocabi(const CMatrix& 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(const CMatrix& 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=(const CMatrix& m) // 赋值运算符{int i,j;for(i=0;i<row;i++){for(j=0;j<column;j++)A[i][j]=m.A[i][j];}return *this;};bool CMatrix::operator ==(const CMatrix& m) // 比括较运算符{ int i,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;};bool CMatrix::operator !=(const CMatrix& m) // 比括较运算符{ int i,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;};CMatrix CMatrix::operator+(const CMatrix& m)// 加运算符{ int i,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-(const CMatrix& m)// 减运算符{ int i,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+=(const CMatrix& m) //自加运算符{ int i,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-=(const CMatrix& 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]-m.A[i][j];return *this;};void CMatrix::find(int& f)const{ int i;for(i=0;i<this->row;i++) if(this->A[i][i]!=0) f=1;else f=0;};CMatrix CMatrix::operator-() // 取负数{ int i,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*(const CMatrix& m) // 乘法运算符{ int i,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()// 输出该矩阵{ int i,j;for(i=1;i<=row;i++){ for(j=1;j<=column;j++)cout<<A[i-1][j-1]<<" ";cout<<endl;}};CMatrix& CMatrix::transpose()// 矩阵转置{ int i,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(const CMatrix& a) //高斯迭代求解{ int f=1;a.find(f);if(f==0) {cerr<<"该矩阵不满足迭代求解条件!"<<endl; exit(1); }else{CMatrix x,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() //列主元素法求解{ CMatrix x,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){ int j,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){ int j,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--;}int j,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() //可逆线性矩阵求解{CMatrix x,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;}}};struct CElvDif{double value; // 观测值double weight; // 权重long startPoint; // 起始点编号long endPoint; // 终点编号};//水准点类的设计struct CLevelPoint{ // 高程平差值=高程值+高程值改正数long index; // 水准点编号double eleValue; // 高程值double dv; // 高程值改正数(初始化为0)bool isKnown; // 是否为已知点};class CElevationNet{public: // 成员函数CElevationNet(){numElvDif=0;numPoints=0 ;numKnPoint=0;}; // 构造函数~CElevationNet(){}; // 析构函数 void input();void output();void jsgc();int getgz(){return numElvDif;} //返回观测值数目int getzs(){return numPoints;} //返回总点数int getys(){return numKnPoint;} //返回已知点数int getws(){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(int i,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: // 成员变量int numElvDif; // 高差总数int numPoints; // 控制网中点的总数目int numKnPoint; //控制网中已知点的数目CElvDif edVec[max]; // 观测值数组CLevelPoint lpVec[max]; // 高程值数组};void CElevationNet::input(){ int a,b,c,i;double z,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;}};void CElevationNet::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(){ int i,j;for(i=0;i<this->numElvDif;i++){ if((lpVec[edVec[i].startPoint-1].eleValu e!=0.0)&&(lpVec[edVec[i].endPoint-1].eleVal ue==0.0)){lpVec[edVec[i].endPoint-1].eleValue=(l pVec[edVec[i].startPoint-1].eleValue)+(edVe c[i].value);}elseif(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)-(edVe c[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].endPoi nt-1].eleValue!=0.0)this->lpVec[j].eleValue=this->lpVec[thi s->edVec[i].endPoint-1].eleValue-this->edVe c[i].value;elseif(this->edVec[i].endPoint==lpVec[j].index& &this->lpVec[this->edVec[i].startPoint-1].e leValue!=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;}};void xishu(CMatrix& B, CMatrix& X,CElevationNet A){int i,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].isKn own==1&&A.lpVec[A.edVec[i].endPoint-1].isKn own==0){B.setValue(i+1,A.edVec[i].endPoint-A.numKn Point,1);}elseif(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);}elseif(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) {int i,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].ele Value+A.edVec[i].value-A.lpVec[A.edVec[i].e ndPoint-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.getgz(),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.getgz(),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));int i,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)/100 0)<<"米"<<endl;vz.operator =(V.transpose());vp=(vz.operator *(Q));p=(vp.operator *(V));cout<<"平差后单位权中误差为:"<<"+_"<<sqrt(p.getValue(1,1)/A.getws())<<"毫米"<<endl;}。
测绘常用程序C语言
测量平差程序设计1.角度(度分秒)到弧度AngleToRadian#define PI 3.14159265double AngleToRadian(double angle){int D,M;double S,radian,degree, angle,MS;D=int(angle+0.3);MS=angle-D;M=int((MS)*100+0.3);S=(MS*100-M)*100;degree=D+M/60.0+S/3600.0;radian=degree*PI/180.0;return radian;}注意:防止数据溢出,要加个微小量,例如0.3.2.弧度换角度(度分秒) RadianToAngle#define PI 3.14159265double RadianToAngle(double radian){int D,M;double S,radian,degree,MS,angle;degree=radian*180/PI;D=int(degree);MS=degree-D;M=int(MS*60);S=(MS*60-M)*60;angle=D+M/100.0+S/10000.0;return angle;}3.已知两点求坐标方位角Azimuth#include <math.h>double Azimuth(double xi,double yi,double xj,double yj) {double Dx,Dy,S,T;Dx=xj-xi;Dy=yj-yi;S=sqrt(Dx*Dx+Dy*Dy);if(S<1e-10) return 0;T=asin(Dy/S);if(Dx<0) T=PI-T;if(Dx>0&&(Dy<0)||T<0) T=2*PI+T;return T;}4.开辟二维数组的动态空间的宏#include <malloc.h>#define NewArray2D(type,A,i,n,m){A=(type**)malloc(n*sizeof(type*));\for(i=0;i<m;i++)\A[i]=(type*)malloc(m*sizeof(type));\}5.释放开辟的二维数组的空间#define FreeSpace(A,i,m){for(i=0;i<m;i++)\free(A[ i]);\free(A);\}注意:释放空间与开辟空间相反,释放空间是先释放列,后释放行.6.矩阵求转置transformmatrixvoid transformmatrix(double **A,double **B,int i,int j){int m,n;for(m=0;m<=i;m++)for(n=0;n<=j;n++){B[n][m]=A[m][n]:}}7.矩阵相乘(mulmatrix)void mulmatrix(double **A,double **B,double **C,int i,int j,int k) {int m,n,p;for(m=0;m<i;m++)for(n=0;n<j;n++){C[m][n]=0;for(p=0;p<k;p++){C[m][n]+=A[m][p]*B[p][n]:}}}8.矩阵求逆(countermatrix)#include <math.h>void countermatrix(double **T, double **s, double **r, double **Q,double **N, double **rt,int n){for(i=0;i<n;i++){s=N[i][i];for(k=0;k<i;k++){s-=T[k][i]*T[k][i];}T[i][i]=sqrt(s)for(j=i+1;j<n;j++){s=N[i][j];for(k=0;k<i;k++){s-=T[k][i]*T[k][j];}T[i][j]=s/T[i][i];}}for(i=0;i<n;i++)for(j=0;j<n;j++){T[i][j]=0;}for(i=n-1;i>=0;i++){r[i][i]=1/T[i][i];for(j=i+1;j<n;j++){s=0;for(k=i;k<j-1;k++){s-=r[i][k]*T[k][j];}r[i][j]=s/T[i][i];}}for(i=0;i<n;i++)for(j=0;j<n;j++){r[i][j]=0;}transformmatrix(r,rt,n,n)mulmatrix(r,rt,Q,n,n)}9.平差主程序之读入数据typedef struct POINT{char name[8];double x,y;int type;}POINT;typedef struct READV ALUE{POINT *begin;POINT *end;double value;}READV ALUE;POINT *GETPOINT(char *name,POINT *pPoint,int nPoint) {int i;for(i=0;i<nPoint;i++){if (strcmp(pPoint[i].name,name)==0)return (pPoint+i)}for(i=0;i<nPoint;i++){if(pPoint[i]=NULL)strcmp(pPoint[i].name,name);pPoint[i].type=0;return(pPoint+i);}}double AngleToRadian(double angle){int D,M;double S,radian,degree, angle,MS;D=int(angle+0.3);MS=angle-D;M=int((MS)*100+0.3);S=(MS*100-M)*100;degree=D+M/60.0+S/3600.0;radian=degree*PI/180.0;return radian;}main(){POINT *pPoint=NULL;READV ALUE *pDirect=NULL;READV ALUE *pDistance=NULL;int nPoint,nKnownPoint,nDirect,nDistance,i;double mo,mf,ms;char begin[8],end[8];FILE *fp=0;fp=fopen(“c:\\dat\\t1.txt”,”r”)fscanf(fp,”%d,%d,%d,%d\n”,&nPoint,&nKnowPoint,&nDirect,&nDistance) if(nPoint>0)pPoint=(POINT*)malloc(nDirect*sizeof(POINT));if(nDirect>0)pDirect=(READV ALUE*)malloc(nDirect*sizeof(READV ALUE));if(nDistance>0)pDistance=(READV ALUE*)malloc(nDistance*sizeof(RAADV ALUE));fscanf(fp,”%lf,%lf,%lf\n”,&mo,&mf,&ms);for(i=0;i<nKnownPoint;i++){fscanf(fp,”%s,%lf,%lf\n”,pPoint[i].name,&pPoint[i].x,&pPoint[i].y);type=1;}for( ;i<nPoint;i++){pPoint[i].name=NULL;pPoint[i].x=0;pPoint[i].y=0;pPoint[i].type=0;}for(i=0;i<nDirect;i++){fscanf(fp,”%s,%s,%lf\n”,begin,end,&pDirect[i].value);pDirect[i].begin=GetPoint(begin,pPoint,nPoint);pDirect[i].end=GetPoint(end,pPoint,nPoint);}for(i=0;i<nDistance;i++){fscanf(fp,”%s,%s,%lf\n”,begin,end,&pDistance[i].value);pDistance[i].begin=GetPoint(begin,pPoint,nPoint);pDistance[i].end=GetPoint(end,pPoint,nPoint);}fclose(fp);}10.角度检验(checkangle)#include <math.h>int checkangle(double angle){int M,S;double MS;if(angle>=0&&angle<360){MS=angle-(int)(angle);if(M<6){S=(int)(MS*1000);if(S%10<6){return 1;}}}return 0;}11.前方交会#define PI=3014159265/***此处调用程序角度换弧度AngleToRadian***/Qianfang(double XE, double YE, double XF, double YF, doubleDEG, double DEF, double DFG, double DFE, double *DFE, double *DFG){double C,A,B;C=DGE-DGF;A=DEF-DEG;B=DFG-DFE;if((C<-PI&&C>-2*PI)||(C>0&&C<PI){XG=(XE/tan(B)+XF/tan(A)-YE+YF)/(1/tan(A)+ 1/tan(B);YG=(YE/tan(B)+YF/tan(A)+XE-XF)/ (1/tan(A)+ 1/tan(B);}if((C>-PI&&C<0)||(C>PI&&C<2*PI)){XG=(XE/tan(B)+XF/tan(A)+YE-YF)/(1/tan(A)+ 1/tan(B);YG=(YE/tan(B)+YF/tan(A)-XE+XF)/ (1/tan(A)+ 1/tan(B);}}12.坐标概算全方向法子函数取出观测方向GetAllDirectint GetAllDirect(char *name,int nDirect,READV ALUE *pDirect, READV ALUE *pStation) {int i,nCount=0;for(i=0;i<nDirect;i++)if(strcmp(pDirect[i].begin->name,name)==0)){pStation[nCount].begin=p(pDirect[nCount].begin;pStation[nCount].end=p(pDirect[nCount].end;pStation[nCount].value=p(pDirect[nCount].value;nCount++;}return nCount;}坐标概算全方向法子程序实现流程(coordinate)coordinate (入口参数设置){READV ALUE pStation[50],pObject[50];int nCount,i,j,k,m,n,p,nobject;for(i=0;i<nPoint;i++){nCount=GetAllDirect(pPoint[i].name,nDirect,pStation)if((nCount>1)||( nCount=1)){for(j=0;j<nCount;j++){if(pStation[j].end->type==1){for(k=0;k<nCount;k++){if(pStation[k].end->type==0)nobject=GetAllDirect(pStation[j].end->name,nDirect,pDirect,pobject)m=-1;n=-1;for(p=0;p<nobject;p++){if(strcmp(pobject[p].end->name,pPoint[i].name)==0){m=p;}if(strcmp(pobject[p].end->name,pStation[k].end->name)==0){n=p;}if(m>=0&&n>=0){pPoint[i]=pStation[k].end-pStation[j].end;pStation[j].end=pObject[m].value-pObject[n].value;{Xe=pPoint[i].x;Ye=pPoint[i].y;Xf=pStation[j].end->x;Yf=pStation[j].end->y;Lef=pStation[j].value;Leg=pStation[k].value;Lfe=pObject[m].value;Lfg=pObject[n].value;Qianfang(Xe,Xf,Ye,Yf,Lef,Leg,Lfe,Lfg,*Xg,*Yg;)pStation[k].end->x=*xg;pStation[k].end->y=*yg;pStation[k].end.type=2;}}}}}}}}13.坐标增量法(calcoordinate)子函数由端点名称得边长值的函数GetDistancedouble GetDistance(char *begin,char *end,int nDistance,READV ALUE *pDistance){int i;for(i=0;i<nDistance;i++){if((strcmp(pDistance[i].begin->name,begin)==0&&strcmp(pDistance[i].end->name,end==0)| |(strcmp(pDistance[i].begin->name,end)==0&&strcmp(pDistance[i].end,begin)==0))return pDistance[i].value;}return -1;}/***函数取出观测方向GetAllDirect***/void calcoordinate(int nDirect,READV ALUE *pDirect,int nDistace,READV ALUE *pDistance,int nPoint,POINT *pPoint){int nPoint,nCount,nDirect,nDistance;int m=-1,i,j,k;double x1,y1,x2,y2,A0,A,S,dx,dy;READV ALUE*pDirect=NULL;READV ALUE pStation[50];for(i=0;i<nPoint;i++){if(pPoint[i].type>0){nCount=GetAllDirect(pPoint[i].name,nDirect,pDirect,pStation[50]);for(j=0;j<nCount;j++){if(pStation[j].end->type>0)m=j;if(m!=-1){for(k=0;k<nCount;k++){if(pStation[k].end->type==0){x1=pPoint[i].x;y1=pPoint[i].y;x2=pStation[j].end->x;y2=pStation[j].end->y;A0=Bearing(x1,y1,x2,y2);A=A0-(DMSToRAD(pStation[m].value)-DMSToRAD(pStation[k].value));if(A<0)A=A+2*PI;if(A>2*PI)A=A-2*PI;S=GetDistance(pPoint[i],pStation[k].end,nDistance,pDistance);if(S<0)continue;else{dx=S*cos(A);dy=S*sin(A);pStation[k].end->x=pPoint[i].x+dx;pStation[k].end->y=pPoint[i].y+dy;pStation[k].end->type=2;}}}}}}}}14.高斯正反算高斯正算:#include <math.h>#include <stdio.h>#define PI 3.14159265double DMSToRAD(double dDMS){int L1,L2;double T,L3;L1=(int)(dDMS+0.3);L2=(int)((dDMS-L1)*100+0.3);L3=((dDMS-L1)*100-L2)*100;T=(L1+L2/60.0+L3/3600.0)*PI/180.0;return T;}void PreGausePositive(double B,double L,double L0, double a, double b, double *N, double *l, double *c, double *t, double *X,double *B1){double a0,a2,a4,a6,a8,m0,m2,m4,m6,m8;double e,e1;e=(sqrt(a*a-b*b))/a;e1=(sqrt(a*a-b*b))/b;B1=DMSToRAD(B);t=tanB1;c=sqrt(e1*e1*cosB1*cos*B1);l=L-L0;N=a/(sqrt(1-e*e*sinB1*sinB1));m0=a*(1-e*e);m2=3/2*e*e*m0;m4=5/4*e*e*m2;m6=7/6*e*e*m4;m8=9/8*e*e*m6;a0=m0+m2/2+3*m4/8+5*m6/16+35*m8/128;a2=m2/2+m4/2+15*m6/32+7/16*m8;a4=m4/8+3*m6/16+7*m8/32;a6=m6/32+m8/16;a8=m8/128;X=a0*B1-a2*(sin(2*B1))/2+a4*(sin(4*B1))/4-a6*(sin(6*B1))/6+a8*(sin(8*B1))/8;}Void BLToXY(double *x,double *y,double N,double l,double c,double t,double B1,double X){x=X+N*l*l*t*cosB1*cosB1*((3+l*l*cosB1*cosB1*(5-t*t+9*c*c+4*c*c*c*c)/4+l*l*\cosB1*cosB1*(61-58*t*t+t*t*t*t)/30))/6;y=N*l*cosB1(1+l*l*cosB1*((1+c*c-t*t)+l*l*cosB1*cosB1(5-18*t*t+t*t*t*t+14*c*c\-58*t*t*c*c)));}高斯反算void XYToBL(double x,double y,double L0,double a,double b,double q,double *B,\double *L){double Bf,c,t,y,N,e1,ee=(sqrt(a*a-b*b))/(a*a);e1=(sqrt(a*a-b*b))/(b*b);for(Bf=0;;){t=tanBf;c=e1*e1*cosBf;N=a/(sqrt(1-e*e*sinBf*sinBf));B=Bf-(1+c*c)*t*y*y/(2*N*N)*(1-y*y)/(12*N*N)*(15+3*t*t+c*c-9*t*t*c*c)-y*y/ (30*N*N)*(61+90*t*t+45*t*t*t*t);if(fabs(B-Bf)<q)break;Bf=B;}L=L0+y/(N*cosBf)*(1-y*y/(6*N*N))*((1+2*t*t+c*c)-y*y/(20*N*N)*(5+6*c*c)+28*t* t+24*t*t*t*t+8*t*t*c*c)-y*y/(42*N*N)*(61+662*t*t+1320*t*t*t*t+720*t*t*t*t*t*t));B=RADToDMS(B);L=RADToDMS(L);}- 11 -。
测量平差程序设计5
测量平差程序设计
For j = nl(i - 1) + 1 To nl(i) '再次遍访i测站 '再次遍访i 上的所有方向值 h = seqn(lb(j)) '查目标点对应的序号 '查目标点对应的序号 If x(h) < 1 Then '该照准方向点坐标未 '该照准方向点坐标未 求出
测量平差程序设计
测量平差程序设计
For j = nl(i - 1) + 1 To nl(i) '遍访i测站上所有方 '遍访i 向值 h = seqn(lb(j)) '查目标点对应的序号 '查目标点对应的序号 If x(h) > 1 Then '目标点坐标已解出 '目标点坐标已解出 t = azimuth(x(i), y(i), x(h), y(h)): p1 = j '反 '反 算方Next循环 '退出For-Next循环 End If Next j
测量平差程序设计
t1 = azimuth(x(m(1)), y(m(1)), x(m(2)), y(m(2))) '利用一对公共点 '利用一对公共点 求旋转角 t2 = azimuth(xo(1), yo(1), xo(2), yo(2)) dt = t2 - t1: x1 = x(m(1)): y1 = y(m(1))
测量平差程序设计
For i = 1 To yds '置入已知点坐标 '置入已知点坐标 x(m(i)) = xo(i): y(m(i)) = yo(i) Next I c = Space(11) & "近 似 坐 标 值" & "近 Chr(13) & Chr(13) & Chr(10) c = c & Space(3) & "点 名" & Space(4) "点 & "X-坐标" & Space(6) & "Y-坐标" & "X-坐标" "Y-坐标" Chr(13) & Chr(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, intcolumn)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{double value; // 观测值double weight; // 权重long startPoint; // 起始点编号long endPoint; // 终点编号};//水准点类的设计structCLevelPoint{ // 高程平差值=高程值+高程值改正数long index; // 水准点编号double eleValue; // 高程值double dv; // 高程值改正数(初始化为0)bool isKnown; // 是否为已知点};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].eleValu e!=0.0)&&(lpVec[edVec[i].endPoint-1].eleVal ue==0.0)){lpVec[edVec[i].endPoint-1].eleValue=(l pVec[edVec[i].startPoint-1].eleValue)+(edVec[i].value);}elseif(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)-(edVe c[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].endPoi nt-1].eleValue!=0.0)this->lpVec[j].eleValue=this->lpVec[thi s->edVec[i].endPoint-1].eleValue-this->edVe c[i].value;elseif(this->edVec[i].endPoint==lpVec[j].index& &this->lpVec[this->edVec[i].startPoint-1].e leValue!=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].isKn own==1&&A.lpVec[A.edVec[i].endPoint-1].isKn own==0){B.setValue(i+1,A.edVec[i].endPoint-A.numKn Point,1);}elseif(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);}elseif(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].ele Value+A.edVec[i].value-A.lpVec[A.edVec[i].e ndPoint-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.getgz(),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.getgz(),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)/100 0)<<"米"<<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、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
getchar();}
if(xiugai==1)
{printf("请输入从%s到%s的距离:",xg1,xg2);
scanf("%lf",&a[cd]);
printf("请输入从%s到%s的高差:",xg1,xg2);
{cd1=i;break;}
if((abs(cd1-cd)!=1)&&(abs(cd1-cd)!=count-2))printf("\n请输入两个相连的点!\n");
if((abs(cd1-cd)==1)||(abs(cd1-cd)==count-2))break;
if(count>1)
{
ioop2:
{ if(xiugai==0)
{printf("请输入从%s到%s的距离:",cz[i-1].a,cz[i].a);
scanf("%lf",&a[i-1]);
printf("请输入从%s到%s的高差:",cz[i-1].a,cz[i].a);
}
continue;
}
if((xg=='n')||(xg=='N'))break;
}
sum1+=a[i-1];
sum2+=b[i-1];
}
} /*点号输入*/
scanf("%lf",&b[cd]);
getchar();} }
if(xiugai==1)
{
printf("/继续修改 y/修改完成 n/ y :");
scanf("%c",&xg);
printf("请确认:");
getchar();
{
printf("\n请正确输入命令!!!!");
printf("/要平差 y/修改数据后再平差 x/不平差并重新输入数据 n/ y :");
scanf("%c",&pcf);
printf("\n请确认:\n");
getchar();
}
printf("请输入第二点点号:");
scanf("%s",xg2);
for(i=1;i<=count;i++)
if(strcmp(xg1,cz[i].a)==0)
{cd=i;break;}
for(i=1;i<=count;i++)
if((pcf=='\n')||(pcf=='y')||(pcf=='Y'))printf("");
if((pcf=='x')||(pcf=='X'))
{
while(1)
{
printf("请输入第一点点号:");
scanf("%s",xg1);
printf("请输入第%d点点号(当输入完时请输入off) :",i+1);
scanf("%s",cz[i+1].a);
getchar();
count++;
i++;}
if((strcmp(wc,cz[i].a)==0)||(strcmp(wc1,cz[i].a)==0))break;
if(abs(su)<=bh)printf("闭合差为| %.1lf |<=%lf(限差)符合标准,可以平差!\n",su,bh);
if(abs(su)>bh)
{printf("\n闭合差为| %.1lf |>%lf超限!\n",su,bh);
printf("你这组数据有问题!\n");
while((xg!='\n')&&(xg!='y')&&(xg!='Y')&&(xg!='n')&&(xg!='N'))
{
printf("\n请正确输入命令!!!!");
printf("/继续修改 y/修改完成 n/ y :");
scanf("%c",&xg);
}
s=sum1/1000;
if(xiugai==0)
{
printf("请输入闭合差超限点数值:");
scanf("%lf",&bh);
getchar();
}
su=sum2*1000;
printf("请输入第二点点号:");
scanf("%s",xg2);
for(i=1;i<=count;i++)
if(strcmp(xg1,cz[i].a)==0)
{cd=i;break;}
for(i=1;i<=count;i++)
if((gs=='y')||(gs=='Y')) /*修改部分*/
{
while(1)
{
printf("请输入第一点点号:");
scanf("%s",xg1);
printf("请输入第二点点号:");
if(xiugBiblioteka i==0) { printf("\n要修改数据吗?");
printf("\n/不需要修改 n/要修改 y/ n :");
scanf("%c",&gs);
printf("请确认:\n");
getchar();
while((gs!='\n')&&(gs!='y')&&(gs!='Y')&&(gs!='n')&&(gs!='N'))
for(i=0;i<100;i++)
{
a[i]=0;
b[i]=0;
}
i=0;
while(1) /*点号输入*/
{
if(xiugai==0)
{
{
printf("\n请正确输入命令!!!!");
printf("\n/不需要修改 n/要修改 y/ n :");
scanf("%c",&gs);
printf("\n请确认:\n");
getchar();
}
if((gs=='\n')||(gs=='n')||(gs=='N'))printf("");
if(xiugai==0)
count--;
if(xiugai==1)
{
sum1=0;
sum2=0;
for(i=1;i<count;i++)
{sum1+=a[i];
sum2+=b[i];}
}
ioop3:
if(xiugai==0)
{
printf("\n请确认:\n");
getchar();
}
if((xg=='\n')||(xg=='y')||(xg=='Y'))
{
while(1)
{
printf("请输入第一点点号:");
scanf("%s",xg1);
if(strcmp(xg2,cz[i].a)==0)
{cd1=i;break;}
if((abs(cd1-cd)!=1)&&(abs(cd1-cd)!=count-2))printf("\n请输入两个相连的点!\n");
if((abs(cd1-cd)==1)||(abs(cd1-cd)==count-2))break;
double a[100],b[100],c[100],d[100],e[100],sum1=0,sum2=0,danwei1,yizhi1,yizhi2,bh,su,s;/*a距离b高差c改正数d平差高差e平查结果*/
static char wc[20]="off",wc1[20]="OFF";
printf("/要平差 y/修改数据后再平差 x/不平差并重新输入数据 n/ y :");
scanf("%c",&pcf);
printf("\n请确认:\n");
getchar();
while((pcf!='\n')&&(pcf!='y')&&(pcf!='Y')&&(pcf!='n')&&(pcf!='N')&&(pcf!='x')&&(pcf!='X'))