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

1、.'原理是用方程 h=b0+b1*x+b2*y+b3*x*x+b4*y*y+b5*x*y 来表达曲面,h指的是高程异常值，比如WGS84到bj54的高程差，然后根据6或者6个以上的公共点求出b0,b1b5,然后如果要求某点的高程值，输入它的x,y就可以得到高程异常值h，然后利用WGS84的BLH中的H加上高程异常值就可以得到54的高程.'这个程序经过2011年01月上旬的实战精度比较高，不过存在一个弱点，就是如果北坐标比较大，如2333444.555，应该先人为的去掉最高位，这样矩阵运算才不会出异常。这是因为矩阵运算的算法不够完善。有空再解决它。'Code By Kis

2、eigo 2011.01.06 Option ExplicitPrivate 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,b2b6这6个参数 Call CalcB0toB6(matA, matB, arrPara) '计算b0,b1,b2b6这6

3、个参数 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个

4、点,这里是高程值,matB(0)是第一个点'输出：B0toB6Out, 下标从0取起，一维数组，下标0-5Public 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空列

5、.因为矩阵运算是从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(U

6、Bound(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 UB

7、ound(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, m

8、atB, 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 B0

9、toB6Out(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行和

10、加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,

11、2) For i = 0 To UBound(a(), 1) For j = 0 To UBound(a(), 2) mtsA(i, j) = a(i, j) Next NextEnd Function'设置测试数据Private Sub SetKnownValueAB(matA() As Double, matB() As Double) '测试的数据 Dim x1 As Double Dim y1 As Double Dim x2 As Double Dim y2 As Double Dim x3 As Double Dim y3 As Double Dim x4 As D

12、ouble Dim y4 As Double Dim x5 As Double Dim y5 As Double Dim x6 As Double Dim y6 As Double Dim x7 As Double Dim y7 As Double Dim z1 As Double Dim z2 As Double Dim z3 As Double Dim z4 As Double Dim z5 As Double Dim z6 As Double Dim z7 As Double x1 = 55 y1 = 66 x2 = 11 y2 = 2 x3 = 22 y3 = 9 x4 = 3 y4

13、= 4 x5 = 13 y5 = 11 x6 = 8 y6 = 17 x7 = 11 y7 = 3 z1 = 6710.2 z2 = 82.6 z3 = 439.6 z4 = 25.2 z5 = 265.7 z6 = 310# z7 = 93.7 matA(0, 0) = 1 matA(1, 0) = 1 matA(2, 0) = 1 matA(3, 0) = 1 matA(4, 0) = 1 matA(5, 0) = 1 matA(6, 0) = 1 matA(0, 1) = x1 matA(1, 1) = x2 matA(2, 1) = x3 matA(3, 1) = x4 matA(4,

14、 1) = x5 matA(5, 1) = x6 matA(6, 1) = x7 matA(0, 2) = y1 matA(1, 2) = y2 matA(2, 2) = y3 matA(3, 2) = y4 matA(4, 2) = y5 matA(5, 2) = y6 matA(6, 2) = y7 matA(0, 3) = x1 * x1 matA(1, 3) = x2 * x2 matA(2, 3) = x3 * x3 matA(3, 3) = x4 * x4 matA(4, 3) = x5 * x5 matA(5, 3) = x6 * x6 matA(6, 3) = x7 * x7

15、matA(0, 4) = y1 * y1 matA(1, 4) = y2 * y2 matA(2, 4) = y3 * y3 matA(3, 4) = y4 * y4 matA(4, 4) = y5 * y5 matA(5, 4) = y6 * y6 matA(6, 4) = y7 * y7 matA(0, 5) = x1 * y1 matA(1, 5) = x2 * y2 matA(2, 5) = x3 * y3 matA(3, 5) = x4 * y4 matA(4, 5) = x5 * y5 matA(5, 5) = x6 * y6 matA(6, 5) = x7 * y7 matB(0

16、, 0) = z1 matB(1, 0) = z2 matB(2, 0) = z3 matB(3, 0) = z4 matB(4, 0) = z5 matB(5, 0) = z6 matB(6, 0) = z7End Sub'计算某位置(a,b)的高程值'输入: a'输入: b'输入: arrParaB0toB6'输出：calcHfitPublic Function calcHfit(ByVal a As Double, ByVal b As Double, arrParaB0toB6() As Double) As Double calcHfit =

17、arrParaB0toB6(0) + arrParaB0toB6(1) * a + arrParaB0toB6(2) * b + _ arrParaB0toB6(3) * a * a + arrParaB0toB6(4) * b * b + arrParaB0toB6(5) * a * b 'calcHfit = 0.2 + 0.3 * a + 0.4 * b + 0.5 * a * a + 0.6 * b * b + 0.7 * a * b 'x=55, y=66 'x=11, y = 2 'x=22, y = 9 'x=3, y = 4 'x

18、=13, y = 11 'x=8, y = 17 'h=6710.200,82.600,439.600,25.200,265.700,310.000End Function还有矩阵运算的代码没有贴，因为sina不让贴那么多代码。 vb矩阵运算的代码：Option Explicit''''''''''''''''''''''''''''''&#

19、39;''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' 模

20、块名：MatrixModule.bas' 功能： 矩阵运算'''''''''''''''''''''''''''''''''''''''''''''''''''''&#

21、39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&#

22、39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&#

23、39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&#

24、39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&#

25、39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''' 功能： 将矩阵转换为显示字符串&

26、#39; 模块名：MatrixModule.bas' 函数名：MatrixToString' 参数： m - Integer型变量，矩阵的行数' n - Integer型变量，矩阵的列数' mtxA - Double型m x n二维数组，存放相加的左边矩阵' sFormat - 显示矩阵各元素的格式控制字符串' 返回值：String型，显示矩阵的字符串''''''''''''''''''''

27、;''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

28、;''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

29、;''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

30、;''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

31、;Function MatrixToString(M As Integer, N As Integer, mtxA() As Double, sFormat As String) As String Dim i As Integer, j As Integer Dim s As String s = "" For i = 1 To M For j = 1 To N s = s + Format(mtxA(i, j), sFormat) + " " Next j s = s + Chr(13) Next i MatrixToString = sEnd Fu

32、nction'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

33、;''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

34、;''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

35、;''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

36、;''''''''''''''''''''''' 功能： 将矩阵的指定行转换为显示字符串' 模块名：MatrixModule.bas' 函数名：MatrixRowToString' 参数： n - Integer型变量，矩阵的列数' r - Integer型变量，要显示的矩阵的行数' mtxA - Double型m x n二维数组，存放相加的左边矩阵' sForma

37、t - 显示矩阵各元素的格式控制字符串' 返回值：String型，显示矩阵指定的行向量''''''''''''''''''''''''''''''''''''''''''''''''''

38、9;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

39、9;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

40、9;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

41、9;'''''''''''''''''''''''''''''Function MatrixRowToString(N As Integer, r As Integer, mtxA() As Double, sFormat As String) As String Dim i As Integer, j As Integer Dim s As String s =

42、"" For j = 1 To N s = s + Format(mtxA(r, j), sFormat) + " " Next j MatrixRowToString = sEnd Function''''''''''''''''''''''''''''''''''''

43、''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

44、''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

45、''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

46、'''''''''''''''''''''''''''''''''''''''''''''' 功能： 将矩阵的指定列转换为显示字符串' 模块名：MatrixModule.bas' 函数名：MatrixColToStr

47、ing' 参数： m - Integer型变量，矩阵的行数' c - Integer型变量，要显示的矩阵的列数' mtxA - Double型m x n二维数组，存放相加的左边矩阵' sFormat - 显示矩阵各元素的格式控制字符串' 返回值：String型，显示矩阵指定的列向量''''''''''''''''''''''''''''

48、;''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

49、;''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

50、;''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

51、;''''''''''''''''''''''''''''''''''''''''''''''''''''Function MatrixColToString(M As Integer

52、, c As Integer, mtxA() As Double, sFormat As String) As String Dim i As Integer, j As Integer Dim s As String s = "" For i = 1 To M s = s + Format(mtxA(i, c), sFormat) + " " Next i MatrixColToString = sEnd Function'''''''''''''&

53、#39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&

54、#39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&

55、#39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&

56、#39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&

57、#39;'''''''' 功能： 计算矩阵的加法' 模块名：MatrixModule.bas' 函数名：MAdd' 参数： m - Integer型变量，矩阵的行数' n - Integer型变量，矩阵的列数' mtxA - Double型m x n二维数组，存放相加的左边矩阵' mtxB - Double型m x n二维数组，存放相加的右边矩阵' mtxC - Double型m x n二维数组，返回和矩阵'''''''

58、''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

59、''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

60、''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

61、''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

62、'''''''''''''Sub MAdd(M As Integer, N As Integer, mtxA() As Double, mtxB() As Double, mtxC() As Double) Dim i As Integer, j As Integer For i = 1 To M For j = 1 To N mtxC(i, j) = mtxA(i, j) + mtxB(i, j) Next j Next iEnd Sub'''''&

63、#39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&

64、#39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&

65、#39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&

66、#39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&

67、#39;'''''''''''''''' 功能： 计算矩阵的减法' 模块名：MatrixModule.bas' 函数名：MSub' 参数： m - Integer型变量，矩阵的行数' n - Integer型变量，矩阵的列数' mtxA - Double型m x n二维数组，存放相减的左边矩阵' mtxB - Double型m x n二维数组，存放相减的右边矩阵' mtxC - Double型m x n二维数组，

68、返回差矩阵'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

69、''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

70、''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

71、''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

72、'''''''''''''''''''''Sub MSub(M As Integer, N As Integer, mtxA() As Double, mtxB() As Double, mtxC() As Double) Dim i As Integer, j As Integer For i = 1 To M For j = 1 To N mtxC(i, j) = mtxA(i, j) + mtxB(i, j) Next j

73、 Next iEnd Sub'''''''''''''''''''''''''''''''''''''''''''''''''''''''''&

74、#39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&

75、#39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&

76、#39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&

77、#39;'''''''''''''''''''''''' 功能： 计算矩阵的数乘' 模块名：MatrixModule.bas' 函数名：MNmul' 参数： m - Integer型变量，矩阵的行数' n - Integer型变量，矩阵的列数' dblNum - Double型变量，乘数' mtxA - Double型m x n二维数组，存放乘数矩阵'

78、; mtxB - Double型m x n二维数组，存放数乘的结果矩阵'''''''''''''''''''''''''''''''''''''''''''''''''''''

79、''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

80、''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

81、''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

82、'''''''''''''''''''''''''''Sub MNmul(M As Integer, N As Integer, dblNum As Double, mtxA() As Double, mtxB() As Double) Dim i As Integer, j As Integer For i = 1 To M For j = 1 To N mtxB(i, j) =

83、 dblNum * mtxA(i, j) Next j Next iEnd Sub''''''''''''''''''''''''''''''''''''''''''''''''''''

84、;''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

85、;''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

86、;''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

87、;'''''''''''''''''''''''''''''' 功能： 计算矩阵的转置' 模块名：MatrixModule.bas' 函数名：MTrans' 参数： m - Integer型变量，矩阵的行数 2009.06.30 Kiseigo 注释 m为输出矩阵mtxAT的行数' n - Integer型变量，矩阵的列数&#

88、39; mtxA - Double型m x n二维数组，存放原矩阵' mtxAT - Double型n x m二维数组，返回转置矩阵''''''''''''''''''''''''''''''''''''''''''''''

89、''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

90、''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

91、''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

92、''''''''''''''''''''''''''''''''''Sub MTrans(M As Integer, N As Integer, mtxA() As Double, mtxAT() As Double) Dim i As Integer, j As Integer For i = 1 To M For j =

93、1 To N mtxAT(i, j) = mtxA(j, i) '2009.06.30 Kiseigo 把右边的mtxAT修改为mtxA Next j Next iEnd Sub'''''''''''''''''''''''''''''''''''''''''&#

94、39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&#

95、39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&#

96、39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&#

97、39;'''''''''''''''''''''''''''''''''''''''' 功能： 计算矩阵的乘法' 模块名：MatrixModule.bas' 函数名：MMul' 参数： m - Integer型变量，相乘的左边矩阵的行数' n -

98、Integer型变量，相乘的左边矩阵的列数和右边矩阵的行数' l - Integer型变量，相乘的右边矩阵的列数' mtxA - Double型m x n二维数组，存放相乘的左边矩阵' mtxB - Double型n x l二维数组，存放相乘的右边矩阵' mtxC - Double型m x l二维数组，返回矩阵乘积矩阵'''''''''''''''''''''''''

99、;''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

100、;''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

101、;''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

102、;'''''''''''''''''''''''''''''''''''''''''''''''''''''''Sub MMul(M As Integer, N

103、 As Integer, l As Integer, mtxA() As Double, mtxB() As Double, mtxC() As Double) Dim i As Integer, j As Integer, k As Integer For i = 1 To M For j = 1 To l mtxC(i, j) = 0# For k = 1 To N mtxC(i, j) = mtxC(i, j) + mtxA(i, k) * mtxB(k, j) Next k Next j Next iEnd Sub'''''''&

104、#39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&

105、#39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&

106、#39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&

107、#39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&

108、#39;'''''''''''''' 功能： 计算复矩阵乘法' 模块名：MatrixModule.bas' 函数名：MCmul' 参数： m - Integer型变量，相乘的左边矩阵的行数' n - Integer型变量，相乘的左边矩阵的列数和右边矩阵的行数' l - Integer型变量，相乘的右边矩阵的列数' mtxAR - Double型m x n二维数组，存放相乘的左边矩阵的实部' mtxAI - Double型m x

109、n二维数组，存放相乘的左边矩阵的虚部' mtxBR - Double型n x l二维数组，存放相乘的右边矩阵的实部' mtxBI - Double型n x l二维数组，存放相乘的右边矩阵的虚部' mtxCR - Double型m x l二维数组，返回矩阵乘积矩阵的实部' mtxCI - Double型m x l二维数组，返回矩阵乘积矩阵的虚部''''''''''''''''''''''&#

110、39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&#

111、39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&#

112、39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''&#

113、39;'''''''''''''''''''''''''''''''''''''''''''''''''''''''''Sub MCmul(M

114、As Integer, N As Integer, l As Integer, mtxAR() As Double, mtxAI() As Double, _mtxBR() As Double, mtxBI() As Double, mtxCR() As Double, mtxCI() As Double) Dim i As Integer, j As Integer, k As Integer Dim p As Double, q As Double, s As Double For i = 1 To M For j = 1 To l mtxCR(i, j) = 0# mtxCI(i, j)

115、 = 0# For k = 1 To N p = mtxAR(i, k) * mtxBR(k, j) q = mtxAI(i, k) * mtxBI(k, j) s = (mtxAR(i, k) + mtxAI(i, k) * (mtxBR(k, j) + mtxBI(k, j) mtxCR(i, j) = mtxCR(i, j) + p - q mtxCI(i, j) = mtxCI(i, j) + s - p - q Next k Next j Next iEnd Sub''''''''''''

116、''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

117、''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

118、''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

119、'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' 功能： 矩阵求逆' 模块名：MatrixModule.bas' 函数名：MRinv' 参数： n - Integer型变量，矩阵的阶数' mtxA - Double型二维数组，体积为n x n。存放原矩阵A；返回时存放其逆矩阵A-1。' 返回值：Boolean型，失败为False，成功为True''''''''''''''