导线网平差.doc

导线网平差.txt老子忽悠孩子叫教育，孩子忽悠老子叫欺骗，互相忽悠叫代沟。▲ 男人 这花花世界，我要用什么颜色来吸引你。 '***********************************************************************************
' 本程序用于单一附（闭）合导线严密平差计算，采用按角度条件平差法。以方向观测值中
'误差的先验值作为单位权中误差。计算结果可求得各待定点的坐标平差值及其点位精度Mx，My，
'及M，并计算出各待定点误差椭圆元素E，F，Z
'
'参考文献： 郭久训．《控制网平差程序设计》　北京：原子能出版社，2004．8
'
'平差数据来源：潘正风等．《数字测图原理与方法》　武汉大学出版社　　186页表6-5
'
'等级：所用平差数据为首级图根导线数据（精度很不高）。本程序中写出了方位角和导线全长相对闭合差的判
'别，但考虑到程序的通用性，将这些限差判断当作了注释处理，而不实际运行。
'
'程序不足：没有导线网的图形表达。
'***********************************************************************************

Private i%, j%, n%, tc#, tb#, B_x!, B_y!, e1!, e2!, m!, m0#, z#, aa As Boolean, bb As Boolean, cc As Boolean ',dd As Boolean
Private Naa#(2, 2), Naa逆#(2, 2), W#(2), K#(2), qq#(2), fx#(2), fy#(2)
Private A#(), Q#(), V#(), C#(), mx#(), my#(), mk#(), e#(), f#(), zz#()
Private Po() As Point

'文件格式说明：
'文件格式详见文件"平差数据.txt"
'
Private Sub 打开文件_Click() '打开文件
Dim ff$, 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
'出错处理
On Error GoTo FileErr
ff = CommonDialog1.FileName 'ff是文件路径名
Open ff For Input As #1 '以顺序文件方式打开文件，使用input

Line Input #1, temp '读取文件中的说明语句
Line Input #1, temp '读取文件中的说明语句
Line Input #1, temp '读取文件中的说明语句
Input #1, n '读取n，n为 （测站数-1）
ReDim Po(n) As Point '定义Po(n)，其中Po(0)存A点数据，Po(n)存B点数据，Po(1)到Po(n-1)存n-1个未知点数据。
Line Input #1, temp '读取文件中的说明语句
Line Input #1, temp '读取文件中的说明语句
Input #1, m, e1, e2 '读取先验方向观测值中误差m，测距仪固定误差e1，比例误差e2
Line Input #1, temp '读取文件中的说明语句
Line Input #1, temp '读取文件中的说明语句
Input #1, C_name, C_x, C_y '读取已知点C
Input #1, A_name, A_x, A_y, A_l, A_s '读取已知点A

For i = 1 To n - 1 '读取n-1个未知点
Input #1, Po(i).name, Po(i).l, Po(i).s
Next i

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点数据
For i = 0 To n
Po(i).l = deg(Po(i).l) '将观测方向左角的单位度分秒化作度
Next i

ReDim Q#(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
Exit Sub
FileErr:
MsgBox "您的文件未打开或打开的文件格式有误！注意：请重新运行本程序！！", , "提示"
End Sub

'
'
'
Private Sub 开始计算_Click()

If bb = True Then
推算方位角和坐标 '求近似坐标
求条件式
'If dd = True Then '若平差数据没有超限，则进行下面的计算
组法方程式
求逆求K
求改正数和平差值
推算方位角和坐标 '求平差后坐标
精度评定
'Else
'MsgBox "限差超限！", , "提示"
'Exit Sub
'End If
Else
MsgBox "数据文件未打开！！", , "提示"
Exit Sub
End If
aa = True
MsgBox "计算完毕！", , "提示"
显示平差结果
End Sub

'
'
'
Private Sub 保存_Click()
Dim ff$

If aa = True Then
CommonDialog1.DialogTitle = "保存平差结果"
CommonDialog1.Filter = "(*txt)|*.txt"
CommonDialog1.FileName = "平差结果.txt"
CommonDialog1.ShowSave
ff = CommonDialog1.FileName
If ff = "" Then
MsgBox "文件名不能为空", , "警告"
Exit Sub
End If

Open ff For Output As #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(度分秒)"
For i = 1 To n - 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), "###.#####")
Next i
Print #2, "-------------------------------------------------------------------------------------"
Close #2

cc = True
MsgBox "保存完毕！", , "提示"
Else
MsgBox "没有需要保存平差结果！", , "提示"
Exit Sub
End If
End Sub

'
'
'
Private Sub 退出_Click()
If aa = True And cc = False Then
提示.Show '若平差结果未保存则进行提示
Else
End '若没有平差结果或平差结果已保存则退出
End If
End Sub

'
'
'
Private Sub 推算方位角和坐标()
Po(0).t = tc + Po(0).l - 180
If Po(0).t > 360 Then Po(0).t = Po(0).t - 360
If Po(0).t < 0 Then Po(0).t = Po(0).t + 360

For i = 1 To n
'求方位角
Po(i).t = Po(i - 1).t + Po(i).l - 180
If Po(i).t > 360 Then Po(i).t = Po(i).t - 360
If Po(i).t < 0 Then Po(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)
Next i
End Sub

'
'角度改正数的单位是秒，边长改正数的单位是厘米
'方向观测值中误差为单位权
'方向观测值中误差的先验值为 m(秒)
'角度观测值的权为 0.5
'距离观测值中误差的先验值为 ms=e1+e2*s*0.0001(厘米)
'
Private Sub 求条件式()
Dim sums!

'求A
For i = 0 To n '求角度改正数系数
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 '求角度观测值的权倒数
Next i
For i = 0 To n - 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 '求距离观测值的权倒数
Next i
'求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
For i = 0 To n - 1
sums = sums + Po(i).s
Next i

'限差判断
'这是首级图根导线的精度要求，其中 n+1 为测站数
'If Abs(W(0)) > 40 * Sqr(n + 1) Or Sqr((W(1) / 100) ^ 2 + (W(2) / 100) ^ 2) / sums > 0.00025 Then '方位角闭合差和导线全长相对闭合差
'Exit Sub
'End If
'dd = True
End Sub

'
'
'
Private Sub 组法方程式()
Dim g%
For i = 0 To 2
For j = i To 2 '由于Naa是对称的，只求其上三角即可。
For g 0 To 2 * n
Naa(i, j) = Naa(i, j) + A(i, g) * A(j, g) * Q(g)
Next g
Next j
Next i
End Sub

'
'
'对于3*3的矩阵来说，利用公式　Naa逆=Naa*/|Naa| 　去求Naa的逆即可，其中Naa*是伴随矩阵。
Private Sub 求逆求K()
Dim ee#, ii#, jj#, kk#, det# 'det是|Naa|的值
ee = 1
For i = 0 To 2
ee = ee * Naa(i, i)
Next i
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逆
For i = 0 To 2
For j = 0 To 2
Naa逆(i, j) = Naa逆(i, j) / det
Next j
Next i
'求K
For i = 0 To 2
For j = 0 To 2
K(i) = K(i) - Naa逆(i, j) * W(j)
Next j
Next i
End Sub

'
'
'
'
Private Sub 求改正数和平差值()
Dim pvv#
pvv = 0
For i = 0 To 2 * n
For j = 0 To 2
V(i) = V(i) + K(j) * A(j, i) * Q(i) '求改正数V
Next j
pvv = pvv + V(i) * V(i) / Q(i) '求 pvv
Next i
m0 = Sqr(pvv / 3) '求单位权中误差 m0
For i = 0 To n - 1
Po(i).s = Po(i).s + V(i) / 100 '求导线边长的平差值，其中边长改正数的单位是厘米，边长的单位是米
Next i
For i = 0 To n
Po(i).l = Po(i).l + V(n + i) / 3600 '求导线左角的平差值，其中角度改正数的单位是秒，角度的单位是度
Next i
End Sub

'
'列出待定点坐标的权函数式，据此进行精度评定。
'
Private Sub 精度评定()
Dim uu#, vv#, ww#, ss#, r#, kk!
For j = 1 To n - 1
uu = 0: vv = 0: ww = 0
For i = 0 To j - 1
C(0, i) = Cos(Po(i).t / p0)
C(1, i) = Sin(Po(i).t / p0)
C(0, n + i) = (Poi).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)
Next i
For kk = 0 To 2
ss = 0: r = 0
For i = 0 To j - 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)
Next i
fx(kk) = ss: fy(kk) = r
Next kk
For kk = 0 To 2
ss = 0: r = 0
For i = 0 To 2
ss = ss - Naa逆(kk, i) * fx(i)
r = r - Naa逆(kk, i) * fy(i)
Next i
K(kk) = ss: qq(kk) = r
Next kk
For i = 0 To 2
uu = uu + fx(i) * K(i)
vv = vv + fy(i) * qq(i)
ww = ww + fx(i) * qq(i)
Next i
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
If z < 0 Then z = z + 90
If ww < 0 Then z = z + 180
zz(j) = dms(z) '点位误差椭圆元素Z
Next j
End Sub

'
'
'
Private Function 显示平差数据(ByVal fname$)
Dim linesfromFile$, Nextline$

Open fname For Input As #1
Do Until EOF(1)
Line Input #1, Nextline
linesfromFile = linesfromFile + Nextline + Chr(13)
Loop
Close #1
Print linesfromFile
End Function

'
'
'
Private Sub 显示平差结果()
Print "平差结果："
Print "--------------------------------------------------------------------------------------"
Print " 单位权中误差 m0=" & Format(m0, "####.##") & " 秒"
Print "---------------------------------------------------------------------------------------------"
Print "点名", "x(m)", "y(m)", "点位误差(cm)", "E(cm)", "F(cm)", "Z(度分秒)"
For i = 1 To n - 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), "###.#####")
Next i
Print "---------------------------------------------------------------------------------------------"
'
'
'
'可以将平差所得的B点坐标和B点的实际坐标输出进行比较，以验证平差结果是否正确
'Print Format(Po(n).x, "########.######"), Format(Po(n).y, "########.######")
'Print B_x, B_y
End Sub