'水準網平差類 LevelNet
Option Explicit
Dim qsd As New DCollection '起算點串列
Dim wzd As New DCollection '未知點串列
Dim sections() As New LevelSection '測段
Dim kNum As Integer '已知點數
Dim sNum As Integer '測段數
Dim msNum As Integer '必要觀測數
Dim fDegree As Integer '自由度
Dim uWeight As Double '單位權
Dim uWeight0 As Double '先驗單位權
Dim tDistance As Double '總距離數
Dim tPoints As Integer '總點數
Dim pvv As Double '用于精度評定
Const eps As Double = 0.00000001 '預設精度
Private Function readQsd(dlist() As String) As DCollection '讀取起算點資料
Dim bm As Benchmark, ds() As String, d As Variant, bs As New DCollection
For Each d In dlist
Set bm = New Benchmark
ds = Split(d, ",")
bm.id = UCase(ds(0))
bm.appH = Val(ds(1))
bm.known = True
If bs.ExistsKey(bm.id) = False Then
bs.Add bm, bm.id
End If
Next
Set readQsd = bs
End Function
Private Function readObs(dlist() As String) As LevelSection() '讀取起算點資料
Dim ss() As New LevelSection, ds() As String, i As Integer, j As Integer
j = UBound(dlist)
ReDim ss(j) As New LevelSection
For i = 0 To j
ds = Split(dlist(1), ",")
ss(i).fromID = UCase(ds(0))
ss(i).toID = UCase(ds(1))
ss(i).Observe = Val(ds(2))
ss(i).Distance = Val(ds(3))
Next i
End Function
Private Function SplitLine(ByVal s As String) As String() '分解格式資料
s = Replace(s, Chr(9), " ")
s = Replace(s, " ", " ")
s = Replace(s, Chr(13), " ")
s = Replace(s, Chr(10), " ")
SplitLine Split(Trim(s))
End Function
Public Sub readData(ByVal obsText As String, qsdText As String, ByVal uw As String) '讀取資料
Dim bs0 As DCollection '起算點串列
Set bs0 = readQsd(SplitLine(qsdText)) '讀取起算點,必須先呼叫
sections = readObs(SplitLine(obsText)) '讀取觀測資料
uWeight0 = Val(uw)
Dim bm As Benchmark, s As Variant
For Each s In sections '根據觀測資料構建水準點集合
Set bm = New Benchmark
bm.id = s.fromID
If bs0.ExistsKey(bm.id) And qsd.ExistsKey(bm.id) = False Then
bm , appH = bs0.Item(bm.id).appH
bm.adjH = bm.appH
bm.known = True
qsd.Add bm, bm.id
ElseIf bs0.ExistsKey(bm.id) = False And wzd.ExistsKey(bm.id) = False Then
bm.appH = -99999
bm.known = False
wzd.Add bm, bm.id
End If
Set bm = New Benchmark
bm.id = s.tolD
If bs0.ExistsKey(bm.id) And qsd.ExistsKey(bm.id) = False Then
bm.appH = bs0.Item(bm.id).appH
bm.adjH = bm.appH
bm.known = True
qsd.Add bm, bm.id
ElseIf bs0.ExistsKey(bm.id) = False And wzd.ExistsKey(bm.id) = False Then
bm.appH = -99999
bm.known = False
wzd.Add bm, bm.id
End If
Next
sNum = UBound(sections) + 1 '測段數
msNum = wzd.Count '必要觀測數=總點數-已知點數
tPoints = kNum + msNum
fDegree = sNum - msNum '自由度=總測段數-必要觀測數
End Sub
Public Property Get sectionNum()
sectionNum = UBound(sections) + 1 '測段數
End Property
Public Property Get mustNum()
mustNum = msNum '必要觀測數
End Property
Public Property Get freeDegree()
freeDegree = fDegree '自由度=總測段數-必要觀測數
End Property
Public Property Get unitWeight()
unitWeight = uWeight '單位權
End Property
Public Property Get totalDistance()
totalDistance = tDistance '總距離
End Property
Public Property Get totalPoints()
totalPoints = tPoints '單位權
End Property
Public Function calcApproximateAltitude() As Boolean '計算近似高程
Dim i As Integer
Dim bm1 As Benchmark, bm2 As Benchmark, s As Variant
For i = 0 To UBound(sections)
For Each s In sections
If qsd.ExistsKey(s.fromID) Then
Set bm1 = qsd.Item(s.fromID)
If wzd.ExistsKey(s.toID) Then
Set bm2 = wzd.Item(s.toID)
bm2.appH = bm1.appH + s.Observe
End If
ElseIf qsd.ExistsKey(s.toID) Then
Set bm2 = qsd.Item(s.toID)
If wzd.ExistsKey(s.fromID) Then
Set bm1 = wzd.Item(s.fromID)
bml.appH = bm2.appH - s.Observe
End If
Set bm1 = wzd.Item(s.fromID)
ElseIf wzd.ExistsKey(s.fromID) And wzd.ExistsKey(s.toID) Then
Set bm1 = wzd.Item(s.fromID)
Set bm2 = wzd.Item(s, .toID)
If (bm1.appH + 9999) > eps And (bm2.appH + 9999) < eps Then
bm2.appH = bm1.appH + s.Observe
ElseIf (bm1.appH + 9999) < eps And (bm2.appH + 9999) > eps Then
bml.appH = bm2.appH - s.Observe
End If
End If
Next
Next i '檢測是否有水準點不能推算
For Each s In sections
If wzd.ExistsKey(s.fromID) Then
Set bm1 = wzd.Item(s.fromID)
If Abs(bml.appH + 9999) < eps Then
MsgBox s.fromID + s.toID + "測段中," + bm1.id + "水準點的近似高程不能推算,請檢查觀測資料。"
calcApproximateAltitude = False
Exit Function
End If
End If
If wzd.ExistsKey(s.toID) Then
Set bm2 = wzd.Item(s.toID)
If Abs(bm2.appH + 9999) < eps Then
MsgBox s.fromID + s.toID + "測段中," + bm2.id + "水準點的近似高程不能推算,請檢查觀測資料。"
calcApproximateAltitude = False
Exit Function
End If
End If
Next
calcApproximateAltitude = True
End Function
Public Sub buildBlP(mB As Matrix, mL As Matrix, mP As Matrix) '建立誤差方程
mB.init sNum, msNum '誤差方程系數矩陣
mL.init sNum, 1 'L系數矩陣
mP.init sNum, sNum '權矩陣
Dim bm1 As Benchmark, bm2 As Benchmark, s As LevelSection
Dim row As Integer
For row = 0 To sNum - 1
Set s = sections(row)
If qsd.ExistsKey(s.fromID) Then
Set bm1 = qsd.Item(s.fromID)
ElseIf wzd.ExistsKey(s.fromID) Then
Set bm1 = wzd.Item(s.fromID)
mB.Element(row, wzd.getindex(bm1, id) - 1) = -1
End If
If qsd.ExistsKey(s.toID) Then
Set bm2 = qsd.Item(s, toID)
ElseIf .wzd.ExistsKey(s.tolD) Then
Set bm2 = wzd.Item(s.toID)
mB.Element(row, 0) = (s.Observe - bm2.appH + bm1.appH) * 1000#
mP.Element(row, row) = uWeight0 / s.Distance
Next row
End Sub
Public Sub Adj(B As Matrix, L As Matrix, P As Matrix)
Dim BT As Matrix, w As Matrix, Nbb As Matrix, Nbb1 As Matrix, i As Integer
Set BT = B.transpose()
Set w = BT.multiply(P).multiply(L) '法方程W
Set Nbb = BT.multiply(P).multiply(B) '法方程N
Set Nbb1 = Nbb.invert() '法方程N的逆陣
Dim x As Matrix, v As Matrix
Set x = Nbb1.multiply(w) '高程改正數
Set v = B.multiply(x).subtract(L) '高差改正數
'精度評定
pvv = v.invert().multiply(P).multiply(v).Element(0, 0)
uWeight = Sqr(pvv / fDegree) '計算單位權
Dim LL As Matrix, bm As Benchmark, s As LevelSection
Set LL = B.multiply(Nbb1).multiply(BT)
For i = 0 To msNum - 1 'miNum為必要觀測數
Set bm = wzd.Item(i + 1)
bm.adjH = bm.appH + x.Element(i, 0) / 1000#
bm.msError = uWeight * Sqr(Nbb1.Element(i, i))
Next i
tDistance = 0
For i = 0 To sNum - 1
Set s = sections(i)
s.CorrectedValue = v.Element(i, 0)
s.AdjustedObserve = s.Observe + s.CorrectedValue / 1000#
s.meanSquareError = uWeight * Sqr(LL.Element(i, i))
tDistance = tDistance + s.Distance
Next i
End Sub
Public Function outResult()
Dim s As String
Dim newline As String, i As Integer, bm As Variant
Dim j As Integer, ts As Variant
newline = Chr(13) + Chr(10)
s = String(90, "-") + newline
s = s + " APPROXIMATE HEIGHT" + newline '輸出近似高程
s = s + String(90, "-") + newline
s = s + " No. Name Height(m) " + newline
s = s + String(90, "-") + newline
j = 0
For Each bm In qsd.Values ' '輸出已知點
j = j + 1
s = s + Format(i, String(6, "@")) + Format(bm.id, String(17, "@"))
s = s + Format(Format(bm.appH, "0.0000"), String(12, "@")) + newline
Next
For Each bm In wzd.Values '輸出未知點
j = j + 1
s = s + Format(j, String(6, "@")) + Format(bm.id, String(17, "@"))
s = s + Format(Format(bm.appH, "0.0000"), String(12, "@")) + newline
Next
s = s + String(90, "-") + newline '輸出已知點高程
s = s + " KNOWN HEIGHT" + newline
s = s + String(90, "-") + newline
For i = 1 To qsd.Count '輸出已知點
Set bm = qsd.Item(i)
s = s + Format(i, String(6, "@")) + Format(bm.id, String(17, "@"))
s = s + Format(Format(bm.appH, "0.0000"), String(12, "@")) + newline
Next i
s = s + String(90, "-") + newline '輸出測段觀測高差
s = s + " MEASURING DATA OF HEIGHT DIFFERENCE" + newline
s = s + String(90, "-") + newline
s = s + " NO. From TO "
s = s + "Observe(m) Distance(km) Weight" + newline
s = s + String(90, "-") + newline
For i = 0 To UBound(sections)
Set ts = sections(i)
s = s + Format(i + 1, String(6, "@")) + Format(ts.fromID, String(18, "@"))
s = s + Format(ts.toID, String(19, "@"))
s = s + Format(Format(ts.Observe, "0.0000"), String(14, "(@"))
s = s + Format(Format(ts.Distance, "0.00"), String(12, "(@"))
s = s + Format(Format(1# / ts.Distance, " 0.000000"), String(15, "(@"))
s = s + newline
Next i
s = s + String(90, "-") + newline '輸出平差高程
s = s + " ADJUSTED HEIGHT" + newline
s = s + String(90, "-") + newline
s = s + "No Name Height(m) Mh(mm)" + newline
s = s + String(90, "-") + newline
j = 0
For Each bm In qsd.Values '輸出已知點
j = j + 1
s = s + Format(j, String(6, "@")) + Format(bm.id, String(18, "@"))
s = s + Format(Format(bm.adjH, " 0.0000"), String(13, "@"))
s = s + Format(Format(bm.msError, "0.0000"), String(13, "@"))
s = s + newline
Next
For Each bm In wzd.Values '輸出未知點
j = j + 1
s = s + Format(j, String(6, "(@")) + Format(bm.id, String(18, "(@"))
s = s + Format(Format(bm.adjH, " 0.0000"), String(13, "(@"))
s = s + Format(Format(bm.msError, " 0.0000"), String(13, "(@"))
s = s + newline
Next
s = s + String(90, "-") + newline
s = s + " ADJUSTED HEIGHT DIFFERENCE" + newline
s = s + String(90, "-") + newline
s = s + " No. From "
s = s + "To Adjusted_dh(m) V(mm) Mdh(mm) " + newline
s = s + String(90, "-") + newline
For i = 0 To UBound(sections)
Set ts = sections(i)
s = s + Format(i + 1, String(6, "@")) + Format(ts.fromID, String(21, "@"))
s = s + Format(ts.toID, String(19, "(@"))
s = s + Format(Format(ts.AdjustedObserve, " 0.0000"), String(20, "@"))
s = s + Format(Format(ts.CorrectedValue, "0.00"), String(10, "@"))
s = s + Format(Format(ts.meanSquareError, "0.00"), String(13, "@")) + newline
Next i
s = s + String(90, "-") + newline
s = s + " UNIT WEIGHT AND PVV " + newline
s = s + String(90, "-") + newline
s = s + " PVV= " + Format(pvv, "0.0000 ") + newline
s = s + " Free Degree=" + Format(fDegree, "####")
uj5u.com熱心網友回復:
代碼不完整,你需要在工程中加入Dcollection和Levelsection這兩個類模塊的代碼。
轉載請註明出處,本文鏈接:https://www.uj5u.com/gongcheng/283578.html
標籤:VB基礎類
上一篇:如何隨機宣告位元組陣列變數?