导线网平差.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
升级会员 