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

1、高程拟合的方法和原理（二次曲面拟合代码）By Kiseigokiseigo 2011-01-06 22:37:14原理是用方程h=b0+b1*x+b2*y+b3*x*x+b4*y*y+b5*x*y 来表达曲面,h指的 是高程异常值，比如 WGS8到bj54的高程差，然后根据6或者6个以上的公共 点求出b0,b1b5,然后如果要求某点的高程值，输入它的x,y就可以得到高程异常值h,然后利用WGS8的 BLH中的H加上高程异常值就可以得到54的高 程这个程序经过2011年01月上旬的实战精度比较高，不过存在一个弱点，就是 如果北坐标比较大，如2333444.555，应该先人为的去掉最高位，这样矩阵

2、运算 才不会出异常。这是因为矩阵运算的算法不够完善。有空再解决它。Code By Kiseigo 2011.01.06Opti on ExplicitPrivate Sub cmdCalc_Click()Dim matA() As DoubleDim matB() As DoubleReDim matA(6, 5) As Double 7 个已知点ReDim matB(6, 0) As DoubleCall SetK nownV alueAB(matA, matB)Dim arrPara() As Double b0,b1,b2b6 这 6 个参数Call CalcB0toB6(matA, m

3、atB, arrPara)计算 b0,b1,b2 b6 这 6 个参数Dim Hout As DoubleHout = 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,

4、 0) 最少 6 个点, 这里是高程值 ,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空行,加

5、 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 的转置 *AReDim ATA(UBound(matAT, 1),

6、 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 IntegerDim j As IntegerFor i = 0 To UBound(ATA, 1)For j = 0 To UBound(ATA, 2) AT

7、Ainv(i, j) = ATA(i, j)Next j输出 ATAinvNext iCall MRinv(UBound(ATAinv), ATAinv) 步骤 4：计算 ATB(A 的转置 *B )Dim ATB() As Double A 的转置 *BReDim 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 的转

8、置 * B)Dim B0toB6OutWithZero() As DoubleReDim B0toB6OutWithZero(UBound(ATAinv, 1), UBound(ATB,2) As DoubleCall MMul(UBound(ATAinv, 1), UBound(ATAinv, 2), UBound(ATB, 2), ATAinv, ATB, B0toB6OutWithZero) 把结果的第一行空行和第一列空列去掉ReDim B0toB6Out(5)B0toB6Out(0) = B0toB6OutWithZero(1, 1)B0toB6Out(1) = B0toB6OutWi

9、thZero(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()ReDi

10、m a(UBound(mtsA(), 1) + 1, UBound(mtsA(), 2) + 1)Dim i As IntegerDim j As IntegerFor i = 0 To UBound(mtsA(), 1)For j = 0 To UBound(mtsA(), 2)a(i + 1, j + 1) = mtsA(i, j)NextNextReDim mtsA(UBound(a, 1), UBound(a, 2)For i = 0 To UBound(a(), 1)For j = 0 To UBound(a(), 2) mtsA(i, j) = a(i, j)NextNextEnd

11、 Function 设置测试数据Private Sub SetKnownValueAB(matA() As Double, matB() As Double) 测试的数据Dim x1 As DoubleDim y1 As DoubleDim x2 As DoubleDim y2 As DoubleDim x3 As DoubleDim y3 As DoubleDim x4 As DoubleDim y4 As DoubleDim x5 As DoubleDim y5 As DoubleDim x6 As DoubleDim y6 As DoubleDim x7 As DoubleDim y7

12、As DoubleDim z1 As DoubleDim z2 As DoubleDim z3 As DoubleDim z4 As DoubleDim z5 As DoubleDim z6 As DoubleDim z7 As Doublex1 =55y1 =66x2 =11y2 =2x3 =22y3 =9x4 =3y4 =4x5 =13y5 =11x6 =8y6 =17x7 =11y7 =3z1 = 6710.2 z2 = 82.6 z3 = 439.6 z4 = 25.2 z5 = 265.7 z6 = 310# z7 = 93.7matA(0, 0) = 1matA(1, 0) = 1

13、matA(2, 0) = 1matA(3, 0) = 1matA(4, 0) = 1matA(5, 0) = 1matA(6, 0) = 1matA(0, 1) = x1 matA(1, 1) = x2 matA(2, 1) = x3 matA(3, 1) = x4 matA(4, 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

14、, 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 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) =

15、 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 * y7matB(0, 0) = z1matB(1, 0) = z2matB(2, 0) = z3matB(3, 0) = z4matB(4, 0) = z5matB(5, 0) = z6matB(6, 0) = z7End Sub 计算某位置 (a,b) 的高程值 输入 : a 输入: b 输入 : arrParaB0toB6 输出：

16、calcHfitPublic Function calcHfit(ByVal a As Double, ByVal b As Double, arrParaB0toB6() As Double) As DoublecalcHfit = arrParaB0toB6(0) + arrParaB0toB6(1) * a + arrParaB0toB6(2) * b + _arrParaB0toB6(3) * a * a + arrParaB0toB6(4) * b * b + arrParaB0toB6(5) * a * bcalcHfit = 0.2 + 0.3 * a + 0.4 * b + 0

17、.5 * a * a + 0.6 * b * b + 0.7* a * bx=55, y=66x=11, y = 2x=22, y = 9x=3, y = 4x=13, y = 11x=8, y = 17 h=6710.200,82.600,439.600,25.200,265.700,310.000End Function还有矩阵运算的代码没有贴，因为 sina 不让贴那么多代码vb 矩阵运算的代码：Option ExplicitII III III III III III Illi III III III III III III III III III III III III Illi I

18、II III II 模块名：MatrixModule.bas 功能：矩阵运算II III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi

19、 III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III功能：将矩阵转换为显示字符串模块名：MatrixModule.bas函数名：MatrixToString 参数： m - Integer 型变量，矩阵的行数 n - Integer 型变量，矩阵的列数 mtxA - Double 型 m x n 二维数组

20、，存放相加的左边矩阵 sFormat - 显示矩阵各元素的格式控制字符串 返回值： String 型，显示矩阵的字符串II III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III

21、 III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III IIFunction MatrixToString(M As Integer, N As Integer, mtxA() As Double, sFormat As String) As StringDim i As Integer, j As IntegerDim s As Strings = For i = 1 To MFor j = 1 To Ns = s + Format(mtxA(i,

22、j), sFormat) + Next j s = s + Chr(13)Next iMatrixToString = sII III III III III III Illi III III III III III III III III III III III III Illi III III IIEnd FunctionII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III I

23、II III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III功能：将矩阵的指定行转换为显示字符串模块名：MatrixModule.bas函数名：MatrixRowToString 参数： n - Integer 型变量，矩阵的列数 r - Integer 型变量，要显示的矩阵的行数 mtxA - Double 型 m x n 二维数组，存放相加的左边矩阵 sFormat - 显示矩阵各元素的格式控制字符串 返回值： Str

24、ing 型，显示矩阵指定的行向量II III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III

25、III III III Illi III III III III III III III III III III IIFunction MatrixRowToString(N As Integer, r As Integer, mtxA() As Double, sFormat As String) As StringDim i As Integer, j As IntegerDim s As Strings = For j = 1 To Ns = s + Format(mtxA(r, j), sFormat) + Next jMatrixRowToString = sEnd Function

26、II III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi

27、 III III III III III III III III III III III功能： 将矩阵的指定列转换为显示字符串模块名： MatrixModule.bas函数名： MatrixColToString 参数： m - Integer 型变量，矩阵的行数c - Integer 型变量，要显示的矩阵的列数mtxA - Double型 m x n 二维数组，存放相加的左边矩阵sFormat -显示矩阵各元素的格式控制字符串返回值： String 型，显示矩阵指定的列向量II III III III III III Illi III III III III III III III III

28、III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III IIFunction MatrixCo

29、lToString(M As Integer, c As Integer, mtxA() As Double, sFormat As String) As StringDim i As Integer, j As IntegerDim s As Strings = For i = 1 To Ms = s + Format(mtxA(i, c), sFormat) + Next iMatrixColToString = sEnd FunctionII III III III III III Illi III III III III III III III III III III III III

30、Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III函数名： MAdd 参数： m - Integer n - In

31、teger mtxA - Double mtxB - Double mtxC - Double功能： 计算矩阵的加法 模块名： MatrixModule.bas型变量，矩阵的行数型变量，矩阵的列数型m x n二维数组，存放相加的左边矩阵 型m x n二维数组，存放相加的右边矩阵 型mxn二维数组，返回和矩阵II III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III

32、 III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III IISub MAdd(M As Integer, N As Integer, mtxA() As Double, mtxB() As Double, mtxC() As Double)Dim i As Inte

33、ger, j As IntegerFor i = 1 To MFor j = 1 To NmtxC(i, j) = mtxA(i, j) + mtxB(i, j)Next jNext iEnd SubII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III

34、 III Illi III III III III III III III III III III III功能： 计算矩阵的减法 模块名： MatrixModule.bas 函数名： MSub参数： m - Integer n - Integer型变量，矩阵的行数型变量，矩阵的列数mtxA - Double mtxB - DoublemtxC - Double型m x n二维数组，存放相减的左边矩阵 型m x n二维数组，存放相减的右边矩阵 型mxn二维数组，返回差矩阵II III III III III III Illi III III III III III III III III III

35、 III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III IISub MSub(M As Intege

36、r, N As Integer, mtxA() As Double, mtxB() As Double, mtxC() As Double)Dim i As Integer, j As IntegerFor i = 1 To MFor j = 1 To NmtxC(i, j) = mtxA(i, j) + mtxB(i, j)Next jNext iEnd SubII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Ill

37、i III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III功能：计算矩阵的数乘模块名：MatrixModule.bas函数名：MNmul参数： m - Integer型变量，矩阵的行数n - Intege

38、r型变量，矩阵的列数dblNum - Double mtxA - Double mtxB - Double型变量，乘数型 m x n 二维数组，存放乘数矩阵型m x n二维数组，存放数乘的结果矩阵II III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III I

39、II Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III IISub MNmul(MAs Integer, N As Integer, dblNum As Double, mtxA() As Double,mtxB() As Double)Dim i As Integer, j As IntegerFor i = 1 To MFor j = 1 To N mtxB(i, j)

40、= dblNum * mtxA(i, j)Next jNext iEnd SubII III III III III III Illi III III III III III III III III III III III III Illi III III II II III III III III III Illi III III III III III III III III III III III III Illi III III II II III III III III III Illi III III III III III III III III III III III III

41、Illi III III II II III III III III III Illi III III III III III III III III III III III 功能：计算矩阵的转置 模块名：MatrixModule.bas 函数名：MTrans 参数：m - Integer 型变量，矩阵的行数 2009.06.30 Kiseigo 注释 m为输出矩阵mtxAT的行数 n - Integer 型变量，矩阵的列数 mtxA - Double 型 m x n 二维数组，存放原矩阵 mtxAT - Double 型 n x m 二维数组，返回转置矩阵II III III III III

42、 III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III I

43、II III III III III III IISub MTrans(MAs Integer, N As Integer, mtxA() As Double, mtxAT() As Double) Dim i As Integer, j As IntegerFor i = 1 To MFor j = 1 To N把右边的 mtxATmtxAT(i, j) = mtxA(j, i) 2009.06.30 Kiseigo 修改为 mtxANext jNext iEnd SubII III III III III III Illi III III III III III III III III I

44、II III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III功能： 计算矩阵的乘法 模块名：

45、MatrixModule.bas函数名： MMul 参数： m - Integer型变量，相乘的左边矩阵的行数n - Integer 型变量，相乘的左边矩阵的列数和右边矩阵的行数l - IntegermtxA - DoublemtxB - DoublemtxC - Double型变量，相乘的右边矩阵的列数型m x n二维数组，存放相乘的左边矩阵 型 n x l 二维数组，存放相乘的右边矩阵 型 m x l 二维数组，返回矩阵乘积矩阵II III III III III III Illi III III III III III III III III III III III III Illi I

46、II III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III IISub MMul(MAs Integer, N As Integer, l A

47、s Integer, mtxA() As Double, mtxB() As Double, mtxC() As Double)Dim i As Integer, j As Integer, k As IntegerFor i = 1 To MFor j = 1 To l mtxC(i, j) = 0#For k = 1 To NmtxC(i, j) = mtxC(i, j) + mtxA(i, k) * mtxB(k, j) Next kNext jNext iEnd SubII III III III III III Illi III III III III III III III III

48、 III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III功能：计算复矩阵乘法模块名：M

49、atrixModule.bas函数名：MCmul参数：m - Integer型变量，相乘的左边矩阵的行数n - Integer型变量，相乘的左边矩阵的列数和右边矩阵的行数l - Integer 型变量，相乘的右边矩阵的列数mtxAR- Double型m x nmtxAI -Double型m x nmtxBR- Double型n x lmtxBI -Double型n x lmtxCR- Double型m x lmtxCI -Double型m x l二维数组，存放相乘的左边矩阵的实部 二维数组，存放相乘的左边矩阵的虚部 二维数组，存放相乘的右边矩阵的实部 二维数组，存放相乘的右边矩阵的虚部 二维数

50、组，返回矩阵乘积矩阵的实部 二维数组，返回矩阵乘积矩阵的虚部II III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III

51、 IIII III III III III III Illi III III III III III III III III III III IISub MCmul(M 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 IntegerDim p As Double,

52、 q As Double, s As DoubleFor i = 1 To MFor j = 1 To l mtxCR(i, j) = 0# mtxCI(i, j) = 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 - qNext kN

53、ext jNext iEnd SubII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III II

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

55、II III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III III III Illi III III IIII III III III III III Illi III III III III III III III III III III IIFunction MRinv(N As Integer, mtxA() As Double) A

56、s Boolean 局部变量ReDim nIs(N) As Integer, nJs(N) As IntegerDim i As Integer, j As Integer, k As IntegerDim d As Double, p As Double 全选主元，消元For k = 1 To Nd = 0#For i = k To NFor j = k To N p = Abs(mtxA(i, j) If (p d) Then d = p nIs(k) = i nJs(k) = jEnd IfNext jNext i 求解失败If (d + 1# = 1#) Then MRinv = False Exit FunctionEnd IfIf (nIs(k) k)