WGS84转BJ54的高程拟合算法公式.doc

原理是用方程 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