水准网间接平差程序设计(C++).doc

/／//／／／//////////／/／//／／////／//／/／//////／////／/／/／/／
// 　 ｖisual　C++6、0 　 编译通过 //
／/／/／//／//／//／/／/／／／／/／////////／//／／／//／/／/／/／／//／/／
//////／//////／/////／／//／//／／///／////////／／/////////
// 　 　参考资料 　 　　 ／/
// 　 　部分网络资料 　 　 　 　//
// 宋力杰《测量平差程序设计》　　　 //
／/姚连壁《基于maｔlａb得控制网平差程序设计》 //
/／/／////／///／////////／//／//／//／////／/／//／／／//／///／/
＃iｎclude〈ｉosｔream>
#incluｄe＜fstｒeaｍ〉
＃iｎclｕdｅ ＜stdlib、h＞
＃inclｕｄe<math、ｈ＞
#iｎｃlude <ioｍaniｐ>
usiｎg namesｐacｅ　ｓｔd；
//／////////////／/／/////／／///／//／///／／//／／//／//／///////／/////／／/／／/////////clａsｓ
ｃlaｓs ＳZWPＣ
{
ｐｒivatｅ：
ｉnt gcz＿zs； 　 　／／高差总数
iｎt　sｚｄ_ｚｓ； 　 //总点数
ﻩiｎt yz＿ｓzd_ｚs；　　 //已知点数
dｏｕｂle　m_pvｖ； 　 　　 //［ｐｖv]
ﻩｉｎt　＊qｓｄ_dh;　 　 //高差起点号
int ＊zd_ｄh; 　 /／高差终点号
cｈａr 　＊＊ｄm； 　　／/点名地址数组
ﻩｄouble　 *gcｚ；　 //观测值数组
ﻩdoｕblｅ　 *ｓｚd＿gc; //高程值数组
ﻩdｏuble　　*P； 　　 　／/观测值得权
double 　*AＴPＡ,*ATＰL; ／/法方程系数矩阵与自由项
dｏuble　 ＊dX;　 //高程改正数、平差值
double *V; 　 　//残差
ﻩdouble　 ｍ_mu; 　 /／单位权中误差
public:
ﻩSZWPC()；
ﻩ~ＳZWＰＣ（）；
ﻩinｔ　ij(ｉｎt　i，inｔ　ｊ);/／对称矩阵下标计算函数
boｏl ｉｎverse（ｄouｂle a[],ｉnt ｎ);/／对称正定矩阵求逆（仅存下三角元素）(参考她人）
ﻩvoｉd 　inpｕtdａta（ｃhａｒ　＊dａta输入原始数据函数
inｔ dｍ_ｄｈ(cｈａr ＊namｅ);　／／点名转点号
ﻩvoｉd 　 ｃa_H０(）; /／近似高程计算函数
void cａ_ATPA（）; //法方程组成函数
ﻩvoｉｄ cａ_dX（); //高程平差值计算函数
voｉｄ 　 prｉｎtrｅsulｔ(char *resulｔfiｌe）; //精度估计与平差值输出函数
doubｌe　ca_V（）； //残差计算函数
　voｉd ｚxeｃpc(ｃhar ＊resulｔ最小二乘平差函数
}；
///////／////／//／／/／////／//////／//////////／///／/／／/////／//／//／/／//////／
／/ 构造函数
SZWPC：:SZWPＣ（）
｛
ｇｃｚ_ｚs=0；
ｓzd＿ｚs＝0;
ｙz_szd_zｓ=0；
}
/／/////／／／／/／/／//////／／/／／//／///／///////／／//／//////////////／///／/／//／/
// 析构函数
SZＷPC：:~SＺＷPＣ()
｛
if(gｃz_zｓ>0)
{
ﻩ deletｅ []qsd_ｄh；
ﻩdｅlete　[］zd_ｄh；
ﻩﻩdelｅte ［]ｇcz；
ﻩdeletｅ [］P；
ﻩﻩdelete [］V；
ﻩ}
ﻩｉf(szｄ_zs〉０)
{
ﻩﻩdｅlｅｔe [］szd_gc;
dｅｌete ［］ATPA；
ﻩ dｅlete　［]AＴPＬ;
deｌete [］dX；
ﻩﻩfor（iｎt　ｉ=0; i〈ｓｚd_ｚs；i+＋)
ﻩﻩ if(dｍ[ｉ］！=NULL）ｄelete[］（dｍ［i]）；
delete ［］ｄm;
ﻩ｝
｝
／////／/////／/////////／/／/／/／///／//／//／／//／////／／//／／/／／/／/////／／／/////////
// 　 对称矩阵下标计算函数
inｔ　SZＷPC：:ij(ｉｎt i，int j)
{
rｅtuｒn （i〉=j）？ ｉ*（i+１）／2+ｊ ：ｊ*（j+1）／2+i;
}
///////////////////／/／/／/／/／////／////／／//／／//////／／／////／///／／／///////／///
// 　对称正定矩阵求逆(仅存下三角元素)（参考她人）
bｏol SZＷPC:：inverｓｅ(dｏuｂｌe a［]，iｎt n）
{
　 double *a0＝new dｏｕble［n］;
　for(iｎt ｋ＝0；k<n;k+＋)
ﻩ｛
ﻩ dｏｕｂle a0０=a［0］；
ﻩﻩif（a０0+1、0==1、0)
ﻩ ｛
ﻩﻩ ｄelete ［]a0；
ﻩ ｒｅｔｕrn ｆａlｓe;
ﻩ｝
ﻩ ｆoｒ(int i=1;ｉ＜n；i++)
｛
ﻩ doｕbｌe ai0 = a[i＊（i+１)／2];
if(ｉ＜=n－k-1）ａ０[i］=　－ai0/a0０；
ﻩﻩ ｅlｓｅ 　 　 a0［i］= aｉ0／ａ０0；
ﻩ for（iｎt j=1;j<＝i；j＋＋）
ﻩ ｛
ﻩﻩ ﻩa［（ｉ—1）*i/2＋j－１]=a[ｉ＊（ｉ+1）/２+j］＋aｉ０＊ａ0[j］；
ﻩ}
ﻩﻩ}
ﻩﻩfoｒ（i=1；i〈n；i++)
ﻩﻩ｛
ﻩ ﻩa［（ｎ—1)＊n/2+i－１］=a0［i］;
ﻩ }
ａ[n*（n+1）/2-1］=１、０/a00；
ﻩ}
ﻩdelｅｔe ［］a０;
return　trｕｅ;
}
/／///////////////／／／／///／／／///／/／／//／／//////／／///／/／／////／//／//／/／//／／／ 原始数据输入函数
voｉｄ SZWPC::inputdａtａ（ｃhaｒ ＊datafile）
{
ifsｔrｅａｍ in）；
ﻩｉｆ（！ infile)
ﻩ｛
ﻩﻩcerr〈＜＂ Oｐｅn errｏr！"〈〈enｄl；
ﻩ}
ﻩinfiｌe＞>ｇｃｚ_zｓ>＞sｚｄ_zｓ〉>ｙｚ_szd＿zs;
ﻩint unPnumber=szd_zs—yz_sｚd_zs；ﻩ
ｓzd_ｇc＝ｎｅw ｄouble ［szd_zs］;
dX=new　doubｌe [ｓzd_zs］；
ﻩＡTPA=ｎｅw douｂle [szd＿zs*（szd＿zｓ+1)／2];
ﻩATPL＝new doｕblｅ　［sｚd_zs］；
qｓd_ｄｈ=new　iｎｔ [gcz_zs]；
ﻩzd_dｈ=new iｎt　[gｃｚ_ｚs］；
ﻩgcｚ＝neｗ douｂｌｅ ［gｃｚ_zs］；
ﻩV=nｅw　doublｅ ［gcz_ｚs］;
P＝nｅw ｄouｂｌe [ｇcz_ｚs］;
dｍ=new　char＊ [szｄ＿zs]；
fｏr(inｔ i=0;i<ｓzd_zs；i+＋）
ﻩ{
ﻩﻩdm［i］ = NＵＬＬ；// dm_dh函数根据dm［ｉ］就是否为NULＬ确定dm［ｉ］就是否为点名地址
}
ﻩchaｒ buffeｒ［12８]； ／/临时数组，保存从文件中读到得点名
for（ i=０;i＜=yz_szd_zs—1；i+＋)／／ 读取已知高程数据
ﻩ{
infile>＞bufｆｅｒ；
ﻩ ｉnt　c=dm＿ｄｈ（bufｆer）;
ﻩ infｉle〉＞szｄ_gc［ｉ］;
ﻩ}
ﻩｆor（i=0;ｉ〈ｇｃz_zs；i＋＋)// 读取观测数据
ﻩ{
ﻩ ｉnfｉle>>buffeｒ； //读取高程起点名
ﻩﻩｑsd_dｈ[i］=dm_ｄh（bｕfｆer）;
ﻩﻩinｆiｌｅ＞＞ｂuffｅｒ;/／读取高程终点
ﻩzｄ＿dｈ[ｉ]=dm＿ｄh（ｂuffer)；ﻩ
ﻩﻩinfile〉＞ｇｃｚ[i]〉〉P［i］； /／读取高差值与路线长度
ﻩP［i]=１、0/P[i]；／/线路长转化为观测值得权
｝
ﻩin();
｝
///////／/／//／／/／////／／//////／//////////／／／///////////////／／/／////／/／ 点名转点号,返回点名对应得点号
int ＳZWPＣ：：ｄm_ｄh（char ＊ｎaｍe）
｛
　 　for(inｔ　i=0；　i〈szd_zs； i+＋）
ﻩ｛
if(dｍ[i］！=ＮULL）
｛
ﻩ ﻩif（ｓtrｃmp（ｎame,dm[i］)==0）rｅｔuｒn i;／/将待查点名与已经存入点名数组得点名比较，若存在返回点号
ﻩ }
ｅlｓe
ﻩ{
ﻩ ﻩｉｎt lｅn =　sｔｒlen(ｎａｍe）;//判断点名长度
ﻩﻩdm[i] = new　ｃhar[len＋１]；//为点名申请存储空间
ﻩﻩﻩstrｃpｙ（dm［i],　ｎaｍe);/／待查点名就是一个新得点名,将新点名得地址放到ｄm数组中
ﻩﻩﻩrｅtｕrn　ｉ；//返回点号
ﻩ｝
ﻩ}
ｒeｔurn －1; //ｄｍ数组已经存满，且没有待查点名
}
//／///／／／/／／／//////////／/////////／//////／//／///／/////／///////////////／//////高程近似值计算
vｏｉｄ ＳZＷPC:：ｃa_H0（）
｛
foｒ（int i=ｙz_szｄ＿zs;ｉ＜szd＿zs；i+＋)ｓｚｄ_gc［i]=—10０00、9;//为计算机设置辨别未知高程点得标志
ﻩfｏｒ(int j=1;；j＋＋)
ﻩ{
ｉnt k=0; //计算出近似高程得点数
ｆoｒ(i=0;ｉ<gcz＿zs；i＋+）
ﻩ {
ﻩﻩﻩｉnｔ ｋ1＝qsd_dh［i]；　 //高差起点号
ﻩﻩint　k2=zd_ｄh［i]； 　 /／高差终点号
ﻩﻩ if（ｓzｄ＿gc［k1]＞—１0000、0　&& szd_gc［k２]〈—1０000、0)//ｋ１点高程或高程近似值已知,k２点高程或高程近似值未知
ﻩ ｛
ﻩﻩ ﻩsｚｄ_gc[ｋ2］=szd_ｇc[k1］+gcz[i］;//计算近似高程
ﻩﻩ k+＋;
ﻩ }
ﻩﻩﻩelsｅ
ﻩ iｆ(ｓzd_gc［k1]<—１0000、0 ＆＆　szd_ｇc[k2］>-10０0０、0）／/ｋ2点高程或高程近似值已知，ｋ１点高程或高程近似值未知
ﻩ ﻩ {
ﻩ ﻩ ｓzｄ_ｇc[k1］＝ｓｚd＿ｇｃ［k２］-gcｚ［i］；//计算近似高程
ﻩｋ＋+;
ﻩﻩ｝
ﻩ ｝
ﻩ ｉf（k==（szd_zｓ－yz_szd_ｚｓ））ｂreak;//所有得近似高程计算完成，退出
}
｝
////////／/////////////／／/／/／/／／////／///／/////／/////////////／/／／／//／／／/////
/／　 组成法方程
ｖoiｄ ＳZWPC::ca_ATPA()
{
　//int t=szd＿ｚs；
fｏr(int i=0; ｉ〈ｓzd_ｚs*(szd_zs+1)/2； i++) ATPA［ｉ］=０、0;／/赋初值
　 fｏr（ｉ=０; i〈szd_zs; i++） AＴPL[i]=０、0；//赋初值
fｏr（iｎｔ ｋ＝0;　k＜ｇcz_zs；　k++)
　 {
　　 int i=qsｄ_dｈ［k］;//获取点号
ﻩﻩiｎt　j=zd_ｄh［k］；//获取点号
　　 　　ﻩdoｕｂle Pk=P［ｋ］;//获取权值
　 ﻩｄouble lk=gcｚ［k]－(ｓzd_ｇc[j］—ｓｚd_gc［i])；／/获得第k个自由项
　 ﻩATＰL［i]—=Ｐk＊ｌk；／/获得法方程自由项
ﻩﻩAＴPL［j]+=Pｋ*ｌk;
ATPA[ij（ｉ,i）]+=Ｐk;//获得法方程系数矩阵
ﻩＡTPA［ij(j,j）］+=Pk；
ﻩＡTＰA［ｉj(ｉ，j)］—=Pｋ；
ﻩ｝
}
/／／//／//／/////／//////／//／//／/／／／///／/／／///／//／/／//／///////／///////／//／////
// 高程平差值计算
ｖoiｄ ＳＺWPＣ：：ca_dX(）
｛
ﻩfｏr（iｎt i=0;i〈yｚ_ｓzd_zs;i＋＋） AＴPA［iｊ（i，i)]=1、０e3０；/／处理已知点
ﻩiｆ（!ｉnvｅrsｅ（ＡＴPA，szd_ｚｓ）)／/矩阵求逆
ﻩ{
ﻩ cerr〈<”法方程系数矩阵降秩！”<<enｄl；/／矩阵为奇异矩阵，无法求逆
exit(0）;//退出程序
}
foｒ(i＝0； ｉ〈szｄ_zs； i+＋）//计算高程改正数
｛
ﻩdouble　ｘi=０、0；
ﻩ ｆｏr（int　j＝0; ｊ〈ｓｚd＿ｚs; j+＋)
｛
ﻩｘｉ+=ＡTPＡ［ij(i,ｊ)］＊ATPL［j]；
ﻩ ｝
ﻩﻩdＸ［ｉ]=xi；
szd＿ｇｃ［i]+=ｘｉ;／/计算高程平差值
｝
｝
／／/／//／/／／////／/////／//／//////／/／/／/／／////／//／///／///／／////／///／////／／////
／/ 　 残差计算
double SＺWＰC:：ca_Ｖ(）
{
ﻩdoubｌe ｐｖv=０、0；
ﻩfor（int i＝0；i〈=gcz_zｓ—1;i+＋)
｛
ﻩ int　ｋ1＝qsd_ｄh［i］；
ﻩｉｎt k2=ｚd＿dｈ[i]；
ﻩ V[i］=szd_gc［ｋ2］-sｚd_gc[k1］－gｃｚ［i］；
ﻩ pｖv+=V［ｉ]*V[i]＊P[ｉ］;
｝
returｎ(pvv）;
｝
//／／/////／//////／/////／////／／／/／//////／/／/／//／/／//／////／／/////／/////／／/／/／
／/　 原始数据与平差值输出
ｖｏｉd SZWＰC：：printresult(ｃhａｒ *ｒｅｓultfiｌｅ）
{
ﻩdouｂｌe pvv=cａ_V(）；　 ／/ 残差计算
ofｓtream out)；//以输出方式打开文件,若文件不存在，创建文件
ﻩ//输出原始观测数据
ｏutfｉle＜<endl〈<＂观测总数：”＜<ｇcz＿ｚs〈<＂ 　　"＜＜＂总点数:＂〈〈ｓzd_ｚｓ;
ouｔｆile<〈" 　 ＂<<"已知点数：”＜<yz_sｚｄ_ｚs〈〈ｅndl；
ﻩoutfile＜<eｎdl<〈＂====＝=========＝=＝===＝ 已知高程 =＝===＝==＝============”<＜enｄl；//输出原始观测数据已知点点号、高程
ｆｏr（iｎt i＝0；i＜=yz＿sｚd_zs—1；i＋+）
{
outfile〈＜" 　＂<<ｄｍ［ｉ］；
outfiｌｅ<<sｅｔｉoｓｆｌags（ioｓ::ｆｉｘeｄ);
ﻩ oｕｔfile＜<seｔw（10)<＜setprecｉsｉon（4)〈＜sｚd_gｃ[i］〈<ｅndｌ;
}
oｕtｆile＜<ｅndl＜＜eｎdl＜〈”＝==＝====＝===＝＝======= 高差观测值=======＝＝=====＝======”<<ｅndl<<ｅnｄl；//输出原始观测数据高程观测值与路线长
ﻩｏｕｔfiｌe<〈”起始点名”＜<” ＂<<"终点点名”＜<＂ 　 ”＜＜"高差观测值（ｍ）”＜<＂ ”<<”两点间距离(km)＂<＜ｅnｄl；
ﻩfｏr（i＝0；i＜=gｃｚ_ｚs-1；i+＋）
ﻩ{
ﻩﻩouｔｆiｌe＜<"　 　"<〈dm［ｑsd_dh[i]］<＜setw（9）〈＜dm［zd＿dｈ[i]]；
ﻩﻩoutfｉlｅ<〈setiosflagｓ（ｉoｓ：：fixeｄ）；
ﻩ outfile<〈ｓeｔｗ（16）〈〈setprｅcision(4）〈<gｃz［i];
ﻩﻩoutfｉle<<ｓetiosfｌags(ｉos：:fixｅd）；
ﻩﻩouｔfiｌe〈<setw（１6）〈〈ｓｅtprecisｉon（４)<<1、0/P［ｉ］<＜ｅndl;
ﻩ}
m_mu=sqrｔ(pvv/(gcz_zｓ－(szｄ_zs-yz_ｓｚd_zs）））；//计算单位权中误差
ﻩoutfiｌe<〈endl<<"===＝===＝==＝======＝＝== 单位权中误差=====＝＝==＝===＝==＝===＝”〈〈enｄl;／/输出单位权中误差
outｆiｌe〈〈ｅndｌ〈〈"σ0="〈<ｍ＿mu〈＜endl;ﻩ
ﻩoｕtfile＜<eｎdｌ<〈”======＝=＝＝=========＝= 高程平差值及其精度 ====＝===＝========＝=＝=”<＜ｅndｌ＜<endl；//输出高程平差值及其精度
outfｉｌe＜〈”点名 近似高程 改正数 高程平差值 中误差＂〈＜endl；
ｆoｒ（ i＝0； ｉ〈szd_ｚs； i＋＋）
{
　 　　　 outfile〈＜setw（2）＜〈dm［i];
ﻩdｏuble dx=dＸ[i]；
ﻩ ｄｏｕｂle ｑiｉ=ＡTPA［ij（ｉ，ｉ）］;
ﻩﻩouｔfｉle＜〈sｅtiｏｓflａgs(ｉos：:fiｘed);
ﻩouｔfｉle〈<ｓetｗ(12）〈＜ｓetｐrecｉsioｎ（4)＜〈sｚd_gc[i］-ｄx；
ﻩﻩouｔfile＜<sｅtｉosｆｌaｇｓ(ｉos::fiｘed）；
ﻩﻩoutfile＜＜setw（10）〈〈seｔｐｒeciｓiｏn（4）＜<dx；
ｏuｔfilｅ<<sｅtiｏsfｌags(ios:：ｆixｅd）；
outｆile<〈setw（１1）〈<setprecisｉon（4）<＜szｄ_ｇc[i];
ﻩﻩoutfilｅ<<setioｓｆｌaｇs（ｉｏs:：fixed)；
ﻩﻩouｔfｉle<〈sｅtw（１0）〈<ｓetprｅcisｉｏn(4)<〈sqrｔ（qｉi)＊m＿mu＜〈endl;
ﻩ}
//输出观测值平差值及其精度
ﻩoｕtfｉle〈<endl<〈ｅndl＜〈”＝＝＝＝＝====＝======＝＝=== 观测值平差值及其精度 ====＝=＝＝==＝＝==＝======＂<<endｌ＜〈endl;
outfile<＜＂起 点 终 点 观测高差 　 v”＜<＂ 　　高差平差值 观测权 　　中误差＂〈<eｎdｌ;
for(ｉ=0；ｉ<=ｇcｚ_zs－1;i＋+）
ﻩ{
ﻩ iｎt k1＝qsd_dh[i］；
ﻩ ｉnt　k2＝zd＿dh［i］；
dｏuｂle qii=ATＰA[ｉj（k1,k１）];
ﻩdouble　qjj=　ATPA［ij(k2,k2）] ；
doubｌｅ ｑij=AＴPＡ［ｉj（ｋ１,k２）］;
ﻩﻩｄoｕble ｍl=ｓｑrｔ(qii+ｑjj－２、０*qｉj）＊m_mu；
ﻩ ouｔ(2)；
ﻩﻩoutfiｌe<<dm［k1]；
ﻩﻩoｕｔ（7）；
outfile〈〈dm[k2]；
ｏutfile＜〈setiosｆlａgｓ(ios：:ｆixｅd）;
ﻩﻩoutfile〈＜ｓｅtw(12）〈〈setpreｃｉｓion（4)<＜ｇｃz［i］；
ﻩﻩoutfｉle<<sｅtiosflａgs(ios:：ｆixed）;
ﻩﻩoutfilｅ<<setw（10）〈＜sｅtprｅcisiｏｎ(4）<＜V[i］;
ouｔｆile＜＜sｅｔioｓｆlags(ioｓ：：fixed)；
ﻩ ouｔfｉle〈〈seｔw（10)〈<setprecision（4）〈〈gｃz[i]+V[i];
outfiｌe<〈sｅtiｏsｆlags（ｉos：：ｆiｘed）；
ﻩoｕtfiｌe〈〈setｗ（1０）<<setpｒeｃisｉｏn（4）〈＜P[i]；
ﻩ ｏutfile〈<ｓeｔiosｆｌags（iｏs:：fixeｄ)；
ﻩ ouｔｆｉle〈＜sｅtw（10）＜＜setｐreciｓion(４）<＜ｍl<<endｌ；
ﻩ}
out（）；
｝
/／//／//／／／//////////////／///／/／////／///／／/////////／///／／////／/／////／////／／
／/ 　水准网最小二乘平差
void SZWPC：:ｚxecｐc(char *resｕｌｔｆile)
｛
ﻩca_H0（）;　//近似高程计算
ca＿ＡTPＡ();　／/ 组成法方程
ｃa_ｄＸ();　／／ 高程平差值计算ﻩ
｝
///////／／///／/／/／//／///／///／//////／／///／///／/／/／/／///／/／////／/／/////／/／
iｎｔ main（)
｛
ﻩcｈar ＊datａｆｉle　＝”算例＼＼Ｄata、txt”；//原始数据文件存储地址指针
ｃhaｒ ＊resｕltfｉｌe =＂算例＼＼Resｕlt、txt”;／/平差结果输出地址指针
coｕt<＜endl〈<endｌ〈〈"水准网经典间接平差"<<ｅnｄl〈〈endｌ；
cout<<＂原数据文件位置：＂〈<dａtafｉlｅ<＜eｎdl;
ﻩcｏｕt＜〈”平差结果文件位置：”〈＜resultfile〈＜eｎdｌ〈＜enｄｌ；
ＳZＷPC new_ｎet;//定义新得对象
ﻩｎew_net、inpｕｔdaｔa(ｄatａ输入原始数据
ﻩｎew＿ｎet、zｘecpc（resulｔ最小二乘平差计算
ｎew_neｔ、priｎｔresuｌt（resｕｌt输出平差结果
　 reｔurn 0;
}