高程拟合的方法和原理(二次曲面拟合代码)

高程拟合的方法和原理(二次曲面拟合代码) By Kiseigo

kiseigo /lvyeqish 2011-01-06 22:37:14

'原理是用方程 h=b0+b1*x+b2*y+b3*x*x+b4*y*y+b5*x*y 来表达曲面,h指的是高程异常值，比如WGS84到bj54的高程差，然后根据6或者6个以上的公共点求出b0,b1……b5,然后如果要求某点的高程值，输入它的x,y就可以得到高程异常值h，然后利用WGS84的BLH中的H加上高程异常值就可以得到54的高程.

'这个程序经过2011年01月上旬的实战精度比较高，不过存在一个弱点，就是如果北坐标比较大，如2333444.555，应该先人为的去掉最高位，这样矩阵运算才不会出异常。这是因为矩阵运算的算法不够完善。有空再解决它。

'Code By Kiseigo 2011.01.06

Option Explicit

Private Sub cmdCalc_Click()

Dim matA() As Double

Dim matB() As Double

ReDim matA(6, 5) As Double '7个已知点

ReDim matB(6, 0) As Double

Call SetKnownValueAB(matA, matB)

Dim arrPara() As Double 'b0,b1,b2……b6这6个参数

Call CalcB0toB6(matA, matB, arrPara) '计算b0,b1,b2……b6这6个参数

Dim Hout As Double

Hout = calcHfit(11, 3, arrPara) '计算某位置的高程，这里刚好取已知点来验算

FrmMain.Caption = Format(Hout, "0.000") '结果得93.7,说明结果正确

End Sub

'求高程拟合(二次曲面拟合)的参数B0,B1,B2,B3,B4,B5,B6 By Kiseigo

2011.01.06 21:53 Helped by BluePan

'输入matA(5,5) 最少6行，也就是最少6个已知高程点

'输入matB(5, 0) 最少6个点,这里是高程值,matB(0)是第一个点

'输出：B0toB6Out, 下标从0取起，一维数组，下标0-5

Public Function CalcB0toB6(matA() As Double, matB() As Double, B0toB6Out() As Double)

'假设方程是 h=b0+b1*x+b2*y+b3*x*x+b4*y*y+b5*x*y; 方程由

BluePan提供

Dim maxPt As Integer '公共点个数,要求>=6个.6表示6个点。

maxPt = UBound(matA, 1) + 1

'步骤1：加1空行,加1空列.因为矩阵运算是从1开始，麻烦

Call RedimMatrisAFrom1Nor0(matA)

Call RedimMatrisAFrom1Nor0(matB)

'步骤2：计算 AT * A 矩阵

Dim matAT() As Double 'A的转置矩阵

ReDim matAT(UBound(matA, 2), UBound(matA, 1))

Call MTrans(UBound(matAT, 1), UBound(matAT, 2), matA, matAT) '求A的转置矩阵

Dim ATA() As Double 'A的转置*A

ReDim ATA(UBound(matAT, 1), UBound(matA, 2)) '方阵

Call MMul(UBound(matAT, 1), UBound(matAT, 2), UBound(matA, 2), matAT, matA, ATA) '计算ATA(A的转置*A )

'步骤3：计算(A的转置*A) 的逆矩阵

Dim ATAinv() As Double 'A的转置*A 的逆矩阵

ReDim ATAinv(UBound(ATA, 1), UBound(ATA, 2))

Dim i As Integer

Dim j As Integer

For i = 0 To UBound(ATA, 1)

For j = 0 To UBound(ATA, 2)

ATAinv(i, j) = ATA(i, j)

Next j

Next i

Call MRinv(UBound(ATAinv), ATAinv) '输出ATAinv

'步骤4：计算ATB(A的转置*B )

Dim ATB() As Double 'A的转置*B

ReDim ATB(UBound(matAT, 1), UBound(matB, 2))

Call MMul(UBound(matAT, 1), UBound(matAT, 2), UBound(matB, 2), matAT, matB, ATB) '计算ATB(A的转置*B )

'步骤5: 计算 (A的转置 * A的逆矩阵) * (A的转置 * B)

Dim B0toB6OutWithZero() As Double

ReDim B0toB6OutWithZero(UBound(ATAinv, 1), UBound(ATB, 2)) As Double

Call MMul(UBound(ATAinv, 1), UBound(ATAinv, 2), UBound(ATB, 2), ATAinv, ATB, B0toB6OutWithZero)

'把结果的第一行空行和第一列空列去掉

ReDim B0toB6Out(5)

B0toB6Out(0) = B0toB6OutWithZero(1, 1)

B0toB6Out(1) = B0toB6OutWithZero(2, 1)

B0toB6Out(2) = B0toB6OutWithZero(3, 1)

B0toB6Out(3) = B0toB6OutWithZero(4, 1)

B0toB6Out(4) = B0toB6OutWithZero(5, 1)

B0toB6Out(5) = B0toB6OutWithZero(6, 1)

End Function

'下标都变成从1开始。0下标的数据都不用，置0.

'相对于在最前加1行和加1列，都是0的行和列

Private Function RedimMatrisAFrom1Nor0(mtsA() As Double)

Dim a()

ReDim a(UBound(mtsA(), 1) + 1, UBound(mtsA(), 2) + 1)

Dim i As Integer

Dim j As Integer

For i = 0 To UBound(mtsA(), 1)

For j = 0 To UBound(mtsA(), 2)

a(i + 1, j + 1) = mtsA(i, j)

Next

Next

ReDim mtsA(UBound(a, 1), UBound(a, 2))