////////////////////////////////////////////////////
// visual C++6.0 编译通过//
////////////////////////////////////////////////////
///////////////////////////////////////////////////
// 参考资料//
// 部分网络资料//
// 宋力杰《测量平差程序设计》//
//姚连壁《基于matlab的控制网平差程序设计》//
///////////////////////////////////////////////////
#include
#include
#include
#include
#include
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
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
double a00=a[0];
if(a00+1.0==1.0)
{
delete []a0;
return false;
}
for(int i=1;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
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!"<
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
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
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
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
{
int k=0; //计算出近似高程的点数
for(i=0;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++;
}
else
if(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
+) ATPA[i]=0.0;//赋初值
for(i=0; i
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
{
cerr<<"法方程系数矩阵降秩!"<
}
for(i=0; i
double xi=0.0;
for(int j=0; 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<
{
outfile<<" "<
outfile<
{
outfile<<" "<
m_mu=sqrt(pvv/(gcz_zs-(szd_zs-yz_szd_zs)));//计算单位权中误差
outfile<
r( i=0; i
outfile<
double qii=ATPA[ij(i,i)];
outfile<
//输出观测值平差值及其精度
outfile<
{
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<
outfile<
outfile.close();
}
//////////////////////////////////////////////////////////////////////////
// 水准网最小二乘平差
void
SZWPC::zxecpc(char *resultfile)
{
ca_H0(); //近似高程计算
ca_ATPA(); //组成法方程
ca_dX(); // 高程平差值计算
}
///////////////////////////////////////////////////////////////////////
int main()
{
char *datafile ="算例\\Data.txt";
//原始数据文件存储地址指针
char *resultfile ="算例\\Result.txt";
//平差结果输出地址指针
cout<
new_net.inputdata(datafile);//输入原始数据
new_net.zxecpc(resultfile);//最小二乘平差计算
new_net.printresult(resultfile);//输出平差结果
return 0;
}