导线网平差.docx
《导线网平差.docx》由会员分享,可在线阅读,更多相关《导线网平差.docx(16页珍藏版)》请在冰豆网上搜索。
导线网平差
导线网平差
'***********************************************************************************
'本程序用于单一附(闭)合导线严密平差计算,采用按角度条件平差法。
以方向观测值中
'误差的先验值作为单位权中误差。
计算结果可求得各待定点的坐标平差值及其点位精度Mx,My,
'及M,并计算出各待定点误差椭圆元素E,F,Z
'
'参考文献:
郭久训.《控制网平差程序设计》 北京:
原子能出版社,2004.8
'
'平差数据来源:
潘正风等.《数字测图原理与方法》 武汉大学出版社 186页表6-5
'
'等级:
所用平差数据为首级图根导线数据(精度很不高)。
本程序中写出了方位角和导线全长相对闭合差的判
'别,但考虑到程序的通用性,将这些限差判断当作了注释处理,而不实际运行。
'
'程序不足:
没有导线网的图形表达。
'***********************************************************************************
Privatei%,j%,n%,tc#,tb#,B_x!
B_y!
e1!
e2!
m!
m0#,z#,aaAsBoolean,bbAsBoolean,ccAsBoolean',ddAsBoolean
PrivateNaa#(2,2),Naa逆#(2,2),W#
(2),K#
(2),qq#
(2),fx#
(2),fy#
(2)
PrivateA#(),Q#(),V#(),C#(),mx#(),my#(),mk#(),e#(),f#(),zz#()
PrivatePo()AsPoint
'文件格式说明:
'文件格式详见文件"平差数据.txt"
'
PrivateSub打开文件_Click()'打开文件
Dimff$,temp$,A_name$,A_x!
A_y!
A_l#,A_s!
B_name$,B_l#,C_name$,C_x!
C_y!
D_name$,D_x!
D_y!
Form1.Cls'清屏
CommonDialog1.DialogTitle="打开数据文件"
CommonDialog1.FileName=""
CommonDialog1.ShowOpen
'出错处理
OnErrorGoToFileErr
ff=CommonDialog1.FileName'ff是文件路径名
OpenffForInputAs#1'以顺序文件方式打开文件,使用input
LineInput#1,temp'读取文件中的说明语句
LineInput#1,temp'读取文件中的说明语句
LineInput#1,temp'读取文件中的说明语句
Input#1,n'读取n,n为(测站数-1)
ReDimPo(n)AsPoint'定义Po(n),其中Po(0)存A点数据,Po(n)存B点数据,Po
(1)到Po(n-1)存n-1个未知点数据。
LineInput#1,temp'读取文件中的说明语句
LineInput#1,temp'读取文件中的说明语句
Input#1,m,e1,e2'读取先验方向观测值中误差m,测距仪固定误差e1,比例误差e2
LineInput#1,temp'读取文件中的说明语句
LineInput#1,temp'读取文件中的说明语句
Input#1,C_name,C_x,C_y'读取已知点C
Input#1,A_name,A_x,A_y,A_l,A_s'读取已知点A
Fori=1Ton-1'读取n-1个未知点
Input#1,Po(i).name,Po(i).l,Po(i).s
Nexti
Input#1,B_name,B_x,B_y,B_l'读取已知点B
Input#1,D_name,D_x,D_y'读取已知点D
Close#1
tc=ZBiaoFSuan(C_x,C_y,A_x,A_y)'坐标反算,求点C到A的坐标方位角,并记作tc,单位是度
tb=ZBiaoFSuan(B_x,B_y,D_x,D_y)'坐标反算,求点B到D的坐标方位角,并记作tb,单位是度
Po(0).name=A_name:
Po(0).x=A_x:
Po(0).y=A_y:
Po(0).l=A_l:
Po(0).s=A_s'Po(0)存A点数据
Po(n).name=B_name:
Po(n).x=B_x:
Po(n).y=B_y:
Po(n).l=B_l'Po(n)存B点数据
Fori=0Ton
Po(i).l=deg(Po(i).l)'将观测方向左角的单位度分秒化作度
Nexti
ReDimQ#(2*n),A#(2,2*n),V#(2*n),C#(1,2*n-2),mx#(n-1),my#(n-1),_
mk#(n-1),e#(n-1),f#(n-1),zz#(n-1)'变量重新定义
bb=True
MsgBox"文件已成功打开",,"提示"
显示平差数据ff
ExitSub
FileErr:
MsgBox"您的文件未打开或打开的文件格式有误!
注意:
请重新运行本程序!
!
",,"提示"
EndSub
'
'
'
PrivateSub开始计算_Click()
Ifbb=TrueThen
推算方位角和坐标'求近似坐标
求条件式
'Ifdd=TrueThen'若平差数据没有超限,则进行下面的计算
组法方程式
求逆求K
求改正数和平差值
推算方位角和坐标'求平差后坐标
精度评定
'Else
'MsgBox"限差超限!
",,"提示"
'ExitSub
'EndIf
Else
MsgBox"数据文件未打开!
!
",,"提示"
ExitSub
EndIf
aa=True
MsgBox"计算完毕!
",,"提示"
显示平差结果
EndSub
'
'
'
PrivateSub保存_Click()
Dimff$
Ifaa=TrueThen
CommonDialog1.DialogTitle="保存平差结果"
CommonDialog1.Filter="(*.txt)|*.txt"
CommonDialog1.FileName="平差结果.txt"
CommonDialog1.ShowSave
ff=CommonDialog1.FileName
Ifff=""Then
MsgBox"文件名不能为空",,"警告"
ExitSub
EndIf
OpenffForOutputAs#2'以顺序文件方式保存文件,使用output
Print#2,"平差结果:
"
Print#2,"---------------------------------------------------------------------------------------------"
Print#2,"单位权中误差m0="&Format(m0,"####.##")&"秒"
Print#2,"---------------------------------------------------------------------------------------------"
Print#2,"点名","x(m)","y(m)","点位误差(cm)","E(cm)","F(cm)","Z(度分秒)"
Fori=1Ton-1
Print#2,""&Po(i).name,Format(Po(i).x,"########.####"),Format(Po(i).y,"########.####"),""&Format(mk(i),"#####.##"),_
""&Format(e(i),"#####.##"),""&Format(f(i),"#####.##"),""&Format(zz(i),"###.#####")
Nexti
Print#2,"---------------------------------------------------------------------------------------------"
Close#2
cc=True
MsgBox"保存完毕!
",,"提示"
Else
MsgBox"没有需要保存平差结果!
",,"提示"
ExitSub
EndIf
EndSub
'
'
'
PrivateSub退出_Click()
Ifaa=TrueAndcc=FalseThen
提示.Show'若平差结果未保存则进行提示
Else
End'若没有平差结果或平差结果已保存则退出
EndIf
EndSub
'
'
'
PrivateSub推算方位角和坐标()
Po(0).t=tc+Po(0).l-180
IfPo(0).t>360ThenPo(0).t=Po(0).t-360
IfPo(0).t<0ThenPo(0).t=Po(0).t+360
Fori=1Ton
'求方位角
Po(i).t=Po(i-1).t+Po(i).l-180
IfPo(i).t>360ThenPo(i).t=Po(i).t-360
IfPo(i).t<0ThenPo(i).t=Po(i).t+360
'求坐标
Po(i).x=Po(i-1).x+Po(i-1).s*Cos(Po(i-1).t/p0)'t要以度为单位
Po(i).y=Po(i-1).y+Po(i-1).s*Sin(Po(i-1).t/p0)
Nexti
EndSub
'
'角度改正数的单位是秒,边长改正数的单位是厘米
'方向观测值中误差为单位权
'方向观测值中误差的先验值为m(秒)
'角度观测值的权为0.5
'距离观测值中误差的先验值为ms=e1+e2*s*0.0001(厘米)
'
PrivateSub求条件式()
Dimsums!
'求A
Fori=0Ton'求角度改正数系数
A(0,n+i)=1
A(1,n+i)=(Po(i).y-Po(n).y)/2062.65
A(2,n+i)=(Po(n).x-Po(i).x)/2062.65
Q(n+i)=2'求角度观测值的权倒数
Nexti
Fori=0Ton-1'求边长改正数系数
A(1,i)=Cos(Po(i).t/p0)
A(2,i)=Sin(Po(i).t/p0)
Q(i)=(e1+e2*Po(i).s*0.0001)^2/m^2'求距离观测值的权倒数
Nexti
'求W
W(0)=(Po(n).t-tb)*3600'单位是s
W
(1)=(Po(n).x-B_x)*100'单位是cm
W
(2)=(Po(n).y-B_y)*100'单位是cm
Fori=0Ton-1
sums=sums+Po(i).s
Nexti
'限差判断
'这是首级图根导线的精度要求,其中n+1为测站数
'IfAbs(W(0))>40*Sqr(n+1)OrSqr((W
(1)/100)^2+(W
(2)/100)^2)/sums>0.00025Then'方位角闭合差和导线全长相对闭合差
'ExitSub
'EndIf
'dd=True
EndSub
'
'
'
PrivateSub组法方程式()
Dimg%
Fori=0To2
Forj=iTo2'由于Naa是对称的,只求其上三角即可。
Forg=0To2*n
Naa(i,j)=Naa(i,j)+A(i,g)*A(j,g)*Q(g)
Nextg
Nextj
Nexti
EndSub
'
'
'对于3*3的矩阵来说,利用公式 Naa逆=Naa*/|Naa| 去求Naa的逆即可,其中Naa*是伴随矩阵。
PrivateSub求逆求K()
Dimee#,ii#,jj#,kk#,det#'det是|Naa|的值
ee=1
Fori=0To2
ee=ee*Naa(i,i)
Nexti
ee=ee+2*Naa(0,1)*Naa(1,2)*Naa(0,2)
ii=Naa(0,1)^2
jj=Naa(0,2)^2
kk=Naa(1,2)^2
det=ee-(ii*Naa(2,2)+jj*Naa(1,1)+kk*Naa(0,0))
'求伴随矩阵Naa*,用Naa逆存储
Naa逆(0,0)=Naa(1,1)*Naa(2,2)-kk
Naa逆(0,1)=Naa(0,2)*Naa(1,2)-Naa(0,1)*Naa(2,2)
Naa逆(1,0)=Naa逆(0,1)
Naa逆(1,1)=Naa(0,0)*Naa(2,2)-jj
Naa逆(0,2)=Naa(0,1)*Naa(1,2)-Naa(0,2)*Naa(1,1)
Naa逆(2,0)=Naa逆(0,2)
Naa逆(2,2)=Naa(0,0)*Naa(1,1)-ii
Naa逆(1,2)=Naa(0,2)*Naa(0,1)-Naa(0,0)*Naa(1,2)
Naa逆(2,1)=Naa逆(1,2)
'求Naa逆
Fori=0To2
Forj=0To2
Naa逆(i,j)=Naa逆(i,j)/det
Nextj
Nexti
'求K
Fori=0To2
Forj=0To2
K(i)=K(i)-Naa逆(i,j)*W(j)
Nextj
Nexti
EndSub
'
'
'
'
PrivateSub求改正数和平差值()
Dimpvv#
pvv=0
Fori=0To2*n
Forj=0To2
V(i)=V(i)+K(j)*A(j,i)*Q(i)'求改正数V
Nextj
pvv=pvv+V(i)*V(i)/Q(i)'求pvv
Nexti
m0=Sqr(pvv/3)'求单位权中误差m0
Fori=0Ton-1
Po(i).s=Po(i).s+V(i)/100'求导线边长的平差值,其中边长改正数的单位是厘米,边长的单位是米
Nexti
Fori=0Ton
Po(i).l=Po(i).l+V(n+i)/3600'求导线左角的平差值,其中角度改正数的单位是秒,角度的单位是度
Nexti
EndSub
'
'列出待定点坐标的权函数式,据此进行精度评定。
'
PrivateSub精度评定()
Dimuu#,vv#,ww#,ss#,r#,kk!
Forj=1Ton-1
uu=0:
vv=0:
ww=0
Fori=0Toj-1
C(0,i)=Cos(Po(i).t/p0)
C(1,i)=Sin(Po(i).t/p0)
C(0,n+i)=(Po(i).y-Po(j).y)/2062.65
C(1,n+i)=(Po(j).x-Po(i).x)/2062.65
uu=uu+C(0,i)^2*Q(i)+C(0,n+i)^2*Q(n+i)
vv=vv+C(1,i)^2*Q(i)+C(1,n+i)^2*Q(n+i)
ww=ww+C(0,i)*C(1,i)*Q(i)+C(0,n+i)*C(1,n+i)*Q(n+i)
Nexti
Forkk=0To2
ss=0:
r=0
Fori=0Toj-1
ss=ss+A(kk,i)*C(0,i)*Q(i)+A(kk,n+i)*C(0,n+i)*Q(n+i)
r=r+A(kk,i)*C(1,i)*Q(i)+A(kk,i+n)*C(1,n+i)*Q(n+i)
Nexti
fx(kk)=ss:
fy(kk)=r
Nextkk
Forkk=0To2
ss=0:
r=0
Fori=0To2
ss=ss-Naa逆(kk,i)*fx(i)
r=r-Naa逆(kk,i)*fy(i)
Nexti
K(kk)=ss:
qq(kk)=r
Nextkk
Fori=0To2
uu=uu+fx(i)*K(i)
vv=vv+fy(i)*qq(i)
ww=ww+fx(i)*qq(i)
Nexti
mx(j)=m0*Sqr(Abs(uu))
my(j)=m0*Sqr(Abs(vv))
mk(j)=Sqr((mx(j))^2+(my(j))^2)'每个点的点位中误差,单位cm
r=uu+vv
ss=Sqr((uu-vv)^2+4*ww^2)
e(j)=m0*Sqr((r+ss)/2)'点位误差椭圆元素E,单位cm
f(j)=m0*Sqr(Abs(r-ss)/2)'点位误差椭圆元素F,单位cm
z=Atn(2*ww/(uu-vv))
z=z/2*p0
Ifz<0Thenz=z+90
Ifww<0Thenz=z+180
zz(j)=dms(z)'点位误差椭圆元素Z
Nextj
EndSub
'
'
'
PrivateFunction显示平差数据(ByValfname$)
DimlinesfromFile$,Nextline$
OpenfnameForInputAs#1
DoUntilEOF
(1)
LineInput#1,Nextline
linesfromFile=linesfromFile+Nextline+Chr(13)
Loop
Close#1
PrintlinesfromFile
EndFunction
'
'
'
PrivateSub显示平差结果()
Print"平差结果:
"
Print"---------------------------------------------------------------------------------------------"
Print"单位权中误差m0="&Format(m0,"####.##")&"秒"
Print"---------------------------------------------------------------------------------------------"
Print"点名","x(m)","y(m)","点位误差(cm)","E(cm)","F(cm)","Z(度分秒)"
Fori=1Ton-1
Print""&Po(i).name,Format(Po(i).x,"########.####"),Format(Po(i).y,"########.####"),""&Format(mk(i),"#####.##"),_
""&Format(e(i),"#####.##"),Format(f(i),"#####.##"),Format(zz(i),"###.#####")
Nexti
Print"---------------------------------------------------------------------------------------------"
'
'
'
'可以将平差所得的B点坐标和B点的实际坐标输出进行比较,以验证平差结果是否正确
'PrintFormat(Po(n).x,"########.######"),Format(Po(n).y,"########.######")
'PrintB_x,B_y
EndSub