当前位置：文档之家› 高程拟合的方法和原理(二次曲面拟合代码)

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

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

kiseigo https://www.doczj.com/doc/bd5410669.html,/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)

ReDim 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)

End 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 Double

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 = 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, 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

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, 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) = z7

End Sub

'计算某位置(a,b)的高程值

'输入: a

'输入: b

'输入: arrParaB0toB6

'输出：calcHfit

Public Function calcHfit(ByVal a As Double, ByVal b As Double, arrParaB0toB6() As Double) As Double

calcHfit = 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=13, y = 11

'x=8, y = 17

'h=6710.200,82.600,439.600,25.200,265.700,310.000

End Function

还有矩阵运算的代码没有贴，因为sina不让贴那么多代码。

vb矩阵运算的代码：

Option Explicit

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' ''''''''''''''''''''

' 模块名：MatrixModule.bas

' 功能：矩阵运算

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' '''''''''''''''''''

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' ''''''''''''''''''''''''''''''''''''''''''''''''''''''

' 功能：将矩阵转换为显示字符串

' 模块名：MatrixModule.bas

' 函数名：MatrixToString

' 参数： m - Integer型变量，矩阵的行数

' n - Integer型变量，矩阵的列数

' mtxA - Double型m x n二维数组，存放相加的左边矩阵

' sFormat - 显示矩阵各元素的格式控制字符串

' 返回值：String型，显示矩阵的字符串

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 = s

End Function

' 功能：将矩阵的指定行转换为显示字符串

' 模块名：MatrixModule.bas

' 函数名：MatrixRowToString

' 参数： n - Integer型变量，矩阵的列数

' r - Integer型变量，要显示的矩阵的行数

' mtxA - Double型m x n二维数组，存放相加的左边矩阵

' sFormat - 显示矩阵各元素的格式控制字符串

' 返回值：String型，显示矩阵指定的行向量

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 = ""

For j = 1 To N

s = s + Format(mtxA(r, j), sFormat) + " "

Next j

MatrixRowToString = s

End Function

' 功能：将矩阵的指定列转换为显示字符串

' 模块名：MatrixModule.bas

' 函数名：MatrixColToString

' 参数： m - Integer型变量，矩阵的行数

' c - Integer型变量，要显示的矩阵的列数

' mtxA - Double型m x n二维数组，存放相加的左边矩阵

' sFormat - 显示矩阵各元素的格式控制字符串

' 返回值：String型，显示矩阵指定的列向量

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''''''''''''''''''''''''''

Function MatrixColToString(M As Integer, 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 = s

End Function

' 功能：计算矩阵的加法

' 模块名：MatrixModule.bas

' 函数名：MAdd

' 参数： m - Integer型变量，矩阵的行数

' n - Integer型变量，矩阵的列数

' mtxA - Double型m x n二维数组，存放相加的左边矩阵

' mtxB - Double型m x n二维数组，存放相加的右边矩阵

' mtxC - Double型m x n二维数组，返回和矩阵

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 i

End Sub

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' ''''''''''''''''''''''''''''''''''''''''''''''''''''''

' 功能：计算矩阵的减法

' 模块名：MatrixModule.bas

' 函数名：MSub

' 参数： m - Integer型变量，矩阵的行数

' n - Integer型变量，矩阵的列数

' mtxA - Double型m x n二维数组，存放相减的左边矩阵

' mtxB - Double型m x n二维数组，存放相减的右边矩阵

' mtxC - Double型m x n二维数组，返回差矩阵

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

Next i

End Sub

' 功能：计算矩阵的数乘

' 模块名：MatrixModule.bas

' 函数名：MNmul

' 参数： m - Integer型变量，矩阵的行数

' n - Integer型变量，矩阵的列数

' dblNum - Double型变量，乘数

' mtxA - Double型m x n二维数组，存放乘数矩阵

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

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) = dblNum * mtxA(i, j)

Next j

Next i

End Sub

' 功能：计算矩阵的转置

' 模块名：MatrixModule.bas

' 函数名：MTrans

' 参数： m - Integer型变量，矩阵的行数 2009.06.30 Kiseigo 注释 m 为输出矩阵mtxAT的行数

' n - Integer型变量，矩阵的列数

' mtxA - Double型m x n二维数组，存放原矩阵

' mtxAT - Double型n x m二维数组，返回转置矩阵

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 = 1 To N

mtxAT(i, j) = mtxA(j, i) '2009.06.30 Kiseigo 把右边的mtxAT 修改为mtxA

Next j

Next i

End Sub

' 功能：计算矩阵的乘法

' 模块名：MatrixModule.bas

' 函数名：MMul

' 参数： m - Integer型变量，相乘的左边矩阵的行数

' n - Integer型变量，相乘的左边矩阵的列数和右边矩阵的行数' l - Integer型变量，相乘的右边矩阵的列数

' mtxA - Double型m x n二维数组，存放相乘的左边矩阵

' mtxB - Double型n x l二维数组，存放相乘的右边矩阵

' mtxC - Double型m x l二维数组，返回矩阵乘积矩阵

Sub MMul(M As Integer, N 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 i

End Sub

' 功能：计算复矩阵乘法

' 模块名：MatrixModule.bas

' 函数名：MCmul

' 参数： m - Integer型变量，相乘的左边矩阵的行数

' n - Integer型变量，相乘的左边矩阵的列数和右边矩阵的行数' l - Integer型变量，相乘的右边矩阵的列数

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

' mtxCI - Double型m x l二维数组，返回矩阵乘积矩阵的虚部

'''''''''''''''''''''''''''''''''''''''''''''''''''''

Sub 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 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) = 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 i

End Sub

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

''''''''''''''''''''''''''''''''''''''''''''''''''''''

' 功能：矩阵求逆

' 模块名：MatrixModule.bas

' 函数名：MRinv

' 参数： n - Integer型变量，矩阵的阶数

' mtxA - Double型二维数组，体积为n x n。存放原矩阵A；返回时存放其逆矩阵A-1。

' 返回值：Boolean型，失败为False，成功为True

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''''''''''''''''''''''''''

Function MRinv(N As Integer, mtxA() As Double) As Boolean

' 局部变量

ReDim nIs(N) As Integer, nJs(N) As Integer

Dim i As Integer, j As Integer, k As Integer

Dim d As Double, p As Double

' 全选主元，消元

For k = 1 To N

d = 0#

For i = k To N

For j = k To N

p = Abs(mtxA(i, j))

If (p > d) Then

d = p

nIs(k) = i

nJs(k) = j

End If

Next j

Next i

' 求解失败

If (d + 1# = 1#) Then

MRinv = False

Exit Function

End If

If (nIs(k) <> k) Then

For j = 1 To N

p = mtxA(k, j)

mtxA(k, j) = mtxA(nIs(k), j)

mtxA(nIs(k), j) = p

Next j

End If

If (nJs(k) <> k) Then

For i = 1 To N

p = mtxA(i, k)

mtxA(i, k) = mtxA(i, nJs(k))

mtxA(i, nJs(k)) = p

Next i

End If

mtxA(k, k) = 1# / mtxA(k, k)

For j = 1 To N

If (j <> k) Then mtxA(k, j) = mtxA(k, j) * mtxA(k, k)

Next j

For i = 1 To N

If (i <> k) Then

For j = 1 To N

If (j <> k) Then mtxA(i, j) = mtxA(i, j) - mtxA(i, k) * mtxA(k, j)

Next j

End If

Next i

For i = 1 To N

If (i <> k) Then mtxA(i, k) = -mtxA(i, k) * mtxA(k, k) Next i

Next k

' 调整恢复行列次序

For k = N To 1 Step -1

If (nJs(k) <> k) Then

For j = 1 To N

p = mtxA(k, j)

mtxA(k, j) = mtxA(nJs(k), j)

mtxA(nJs(k), j) = p

Next j

End If

If (nIs(k) <> k) Then

For i = 1 To N

p = mtxA(i, k)

mtxA(i, k) = mtxA(i, nIs(k))

mtxA(i, nIs(k)) = p

Next i

End If

Next k

' 求解成功

MRinv = True

End Function

' 功能：复矩阵求逆

' 模块名：MatrixModule.bas

' 函数名：MCinv

' 参数： n - Integer型变量，矩阵的阶数

' mtxAR - Double型二维数组，体积为n x n。存放原矩阵A的实部；返回时存放其逆矩阵A-的实部。

' mtxAI - Double型二维数组，体积为n x n。存放原矩阵A的虚部；返回时存放其逆矩阵A-的虚部。

' 返回值：Boolean型，失败为False，成功为True

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' ''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''''''''''''''''''''''''''

Function MCinv(N As Integer, mtxAR() As Double, mtxAI() As Double) As Boolean

' 局部变量

ReDim nIs(N) As Integer, nJs(N) As Integer

Dim i As Integer, j As Integer, k As Integer

Dim d As Double, p As Double, s As Double, t As Double, q As Double, b As Double

' 全选主元，消元

For k = 1 To N

d = 0#

For i = k To N

For j = k To N

p = mtxAR(i, j) * mtxAR(i, j) + mtxAI(i, j) * mtxAI(i, j)

If (p > d) Then

d = p

nIs(k) = i

nJs(k) = j

End If

Next j

Next i

' 求解失败

If (d + 1# = 1#) Then

MCinv = False

Exit Function

End If

If (nIs(k) <> k) Then

For j = 1 To N

t = mtxAR(k, j)

mtxAR(k, j) = mtxAR(nIs(k), j)

mtxAR(nIs(k), j) = t

t = mtxAI(k, j)

mtxAI(k, j) = mtxAI(nIs(k), j)

mtxAI(nIs(k), j) = t

Next j

End If

If (nJs(k) <> k) Then

For i = 1 To N

t = mtxAR(i, k)

mtxAR(i, k) = mtxAR(i, nJs(k))

mtxAR(i, nJs(k)) = t

t = mtxAI(i, k)

mtxAI(i, k) = mtxAI(i, nJs(k))

mtxAI(i, nJs(k)) = t

Next i

End If

mtxAR(k, k) = mtxAR(k, k) / d

mtxAI(k, k) = -mtxAI(k, k) / d

For j = 1 To N

If (j <> k) Then

p = mtxAR(k, j) * mtxAR(k, k)

q = mtxAI(k, j) * mtxAI(k, k)

s = (mtxAR(k, j) + mtxAI(k, j)) * (mtxAR(k, k) + mtxAI(k, k))

mtxAR(k, j) = p - q

mtxAI(k, j) = s - p - q

End If

Next j

For i = 1 To N

If (i <> k) Then

For j = 1 To N

If (j <> k) Then

p = mtxAR(k, j) * mtxAR(i, k)

q = mtxAI(k, j) * mtxAI(i, k)

s = (mtxAR(k, j) + mtxAI(k, j)) * (mtxAR(i, k) + mtxAI(i, k))

t = p - q

b = s - p - q

mtxAR(i, j) = mtxAR(i, j) - t

mtxAI(i, j) = mtxAI(i, j) - b

End If

Next j

End If

Next i

For i = 1 To N

If (i <> k) Then

p = mtxAR(i, k) * mtxAR(k, k)

q = mtxAI(i, k) * mtxAI(k, k)

s = (mtxAR(i, k) + mtxAI(i, k)) * (mtxAR(k, k) + mtxAI(k, k))

mtxAR(i, k) = q - p

mtxAI(i, k) = p + q - s

End If

Next i

Next k

' 调整恢复行列次序

For k = N To 1 Step -1

If (nJs(k) <> k) Then

For j = 1 To N

t = mtxAR(k, j)

mtxAR(k, j) = mtxAR(nJs(k), j)

mtxAR(nJs(k), j) = t

t = mtxAI(k, j)

mtxAI(k, j) = mtxAI(nJs(k), j)

mtxAI(nJs(k), j) = t

Next j

End If

If (nIs(k) <> k) Then

For i = 1 To N

t = mtxAR(i, k)

mtxAR(i, k) = mtxAR(i, nIs(k))

mtxAR(i, nIs(k)) = t

t = mtxAI(i, k)

mtxAI(i, k) = mtxAI(i, nIs(k))

mtxAI(i, nIs(k)) = t

Next i

End If

Next k

' 求解成功

MCinv = True

End Function

' 模块名：MatrixModule.bas

' 函数名：MSsgj

' 功能：计算对称正定矩阵的逆

' 参数： n - Integer型变量，矩阵阶数。

' mtxA - Double型二维数组，体积为n x n。存放实对称正定矩阵A；返回时存放逆矩阵A-。

' 返回值：Boolean型，成功为True，失败为False。

Function MSsgj(N As Integer, a() As Double) As Boolean

' 局部变量

Dim i As Integer, j, k, M

Dim W As Double, g As Double

ReDim b(N) As Double

' 循环重新编号求解

For k = 1 To N

W = a(1, 1)

' 求解失败

If (Abs(W) + 1# = 1#) Then

MSsgj = False

Exit Function

End If

M = N - k - 1

For i = 2 To N

g = a(i, 1)

b(i) = g / W

If (i <= M) Then b(i) = -b(i)

For j = 2 To i

a(i - 1, j - 1) = a(i, j) + g * b(j)

Next j

Next i

a(N, N) = 1# / W

For i = 2 To N

a(N, i - 1) = b(i)

Next i

Next k

For i = 1 To N - 1

For j = i + 1 To N

a(i, j) = a(j, i)

Next j

Next i

' 求解成功

MSsgj = True

End Function

' 模块名：MatrixModule.bas

' 函数名：MTinv

' 功能：用特兰持（Trench）方法求托伯利兹（Toeplitz）矩阵的逆矩阵

' 参数： n - Integer型变量，T型矩阵阶数。

' dblT - Double型一维数组，长度为n。存放n阶T型矩阵中的上三角元素t0,t1,…tn-1。

' dblTT - Double型一维数组，长度为n。其中后n-1个元素

tt(1),…,tt(n-1)依次存放n阶T型矩阵中的元素。

' dblB - Double型二维数组，体积为n x n。返回n阶T型矩阵的逆矩阵。

' 返回值：Boolean型，成功为True，失败为False。

Function MTinv(N As Integer, dblT() As Double, dblTT() As Double, dblB() As Double) As Boolean

' 局部变量

Dim i As Integer, j As Integer, k As Integer

Dim a As Double, s As Double

ReDim c(N) As Double, r(N) As Double, p(N) As Double

' 矩阵非T型矩阵

If (Abs(dblT(1)) + 1# = 1#) Then

MTinv = False

Exit Function

End If

' 取初值

a = dblT(1)

c(1) = dblTT(2) / dblT(1)

r(1) = dblT(2) / dblT(1)

' 循环计算

For k = 1 To N - 2

s = 0#

For j = 2 To k + 1

s = s + c(k + 1 - j + 1) * dblTT(j)

Next j

s = (s - dblTT(k + 2)) / a

GPS水准高程拟合报告

GPS水准高程拟合报告实验目的： 1掌握GPS水准高程拟合的基本原理，了解高精度GPS水准的研究意义； 2能够利用Matlab编程实现几何内插法拟合GPS水准高程；实验内容：利用Matlab编程实现几何内插法拟合GPS水准高程，并作内插结果分析实验原理： 1大地水准面，参考椭球面，正高，大地高之间的几何关系 A 正高的定义是：由地面点沿通过该点的铅垂线到大地水准面的距离。一般用符号Hg表示。 B 大地高的定义是：由地面点沿通过该点的椭球面法线到椭球面的距离。也称为椭球高，一般用符号H表示。大地高是一个纯几何量，不具有物理意义。同一个点，在不同的基准下，具有不同的大地高。利用GPS，可以测定地面点的WGS-84中的大地高。 C 大地水准面差距：大地水准面到椭球面的距离，称为大地水准面差距，记为hg （或N）。如上图可以看出大地高和正高之间的关系：H=Hg+ hg 2几何内插法原理几何内插法是通过一些既进行了GPS观测又具有水准资料的点上的大地水准面差距，采用平面或曲面拟合，配置三次样条等内插方法，得到其他点上的大地水准面差距从而反算这些点上的正高。 3二多项式拟合 N=a0+a1*dB+a2*dL+a3*dB2+a4*dL2+a5*dL*dB 公式一

式中dB=B-B0;dL=L-L0;B0=1/n∑B;L0=1/n∑L,n为GPS观测点的数量。利用其中一些具有水准观测资料的公共点上的的大地高和正高可以计算出这些点的大地水准面差距。利用这些公共点的观测资料求得公式一的参数，再利用求得的公式进行其他点的大地水准面差距内插，和正高的拟合；实验步骤： 1输入已知点的GPS观测值和相应的正常高构成矩阵B,L,H,h,分别是纬度矩阵，经度矩阵，大地高矩阵，正高矩阵； 2计算dB=B-B0;dL=L-L0;B0=1/n∑B;L0=1/n∑L，构成矩阵矩阵dB,dL和大地水准面差距矩阵N=H-h; 3将以上计算得到的矩阵代入公式一经过间接平差求得相应的参数a i,这样就能构成一个确定的多项式二； 4输入未知待求点的GPS观测值构成矩阵BB,LL,HH，计算相应的dBB,dLL； 5将dBB,dLL矩阵代入多项式二，解算出对应点的大地水准面差距NN矩阵； 6反算各点的正高h=H-NN; 7对计算得到的正高，大地水准面差距做对比分析；实验分析： 1本实验中可以选择两种差值公式算法（1）N=a0+a1*dB+a2*dL+a3*dB2+a4*dL2+a5*dL*dB （2）N=a0+a1*B+a2*L+a3*B2+a4*L2+a5*L*B 采用公式（1）的插值结果如下： Δh(dB)散点图注：Δh(dB)是插值点的水准资料与插值结果的差值采用公式（2）的插值结果如下：

曲面拟合实例教程总结

例7.2.1试用最小二乘法求拟合曲线，并估计其误差，做出拟合曲线。 (1)做散点图 x=[-2.5,-1.7,-1.1,-0.8,0,0.1,1.5,2.7,3.6]; y=[-192.9,-85.50,-36.15,-26.52,-9.10,-8.43,-13.12,6.50,68.04]; plot(x,y,'r*') legend('实验数据(xi,yi)') xlabel('x'),ylabel('y') title('例7.2.1的数据点(xi,yi)的散点图') 2.CFTOOL工具箱使用 Shift+enter:换行输入 Gaussian：高斯曲线 Interpolant：最小二乘法差值 Polynomial：多项式 3.y1=polyfit(x,y,3) 拟合多项式的阶数为3 4.matlab绘制三维曲面图已知曲线关系方程以二元函数图z = xexp(-x^2-y^2) 为例讲解基本操作，（1）首先需要利用meshgrid函数生成X-Y平面的网格数据，如下所示： % 生成二维网格数据 xa = [-2,0.2,2]; ya =[-1,0.15,1.5]; [x,y] = meshgrid(xa,ya); （2）此外，需要计算纵轴数据(z轴），如下所示： % calculate z data z = x.*exp(-x.^2 - y.^2); （3）在计算出(x,y,z)数据后，就可以使用三维绘图函数mesh绘制三维曲面图，如下所示：mesh(x,y,z); 4（2）、另一种方法： [x,y] = meshgrid(-2:0.2:2,-1:0.15:1.5); z = x.*exp(-x.^2 - y.^2); mesh(x,y,z); 5.由三组散点图绘制曲面（网格划分） xyz=[40 2 1.4 40 5 2.5 40 7 1.4 40 9 0.9 70 8 5.6 ]; tri = delaunay(xyz(:,1), xyz(:,2));

GPS高程拟合方法及其应用

GPS高程拟合方法及其应用论文介绍了GPS高程拟合的原理。介绍了多种拟合模型的拟合原理、模型参数的优化选择，给出了利用地表拟合求解较高精度高程异常的方法，将各种模型进行应用对比。标签：大地高GPS水准高程异常拟合模型 1 GPS高程异常当前GPS技术在平面控制测量工作中已经得到了广泛的应用，但在高程控制测量中却未能得到广泛应用。原因是GPS高程测量得到的是建立在WGS-84坐标系上的大地高H，而我国测量工作中采用的是正常高H。GPS高程测量可以获得厘米级精度的大地高，但在GPS大地高转换为正常高过程中，由于未能获得同等精度的高程异常ζ，导致转换所得的GPS正常高达不到精度要求。 2高程拟合常用方法拟合法是对GPS观测点进行几何水准联测，同一点的大地高减去正常高得到该点的高程异常，再把测区的似大地水准面假定为多项式曲面或者其他数学曲面去拟合已知高程异常的点，根据拟合的曲面内插其他GPS点的高程异常值。拟合法进行GPS高程转换的数学模型很多，如多项式曲线拟合、最小二乘平面拟合、二次多项式曲面拟合等，归纳起来可以分为线状拟合模型、平面拟合模型和曲面线状拟合模型三类。 3高程拟合实例分析一测区，选取其中32个GPS水准高程点进行拟合，将32个水准点的X与Y值通过AutoCAD一个简短的VB加载程序展绘成图：方案一：16个起算点均匀分布选取点2，4，8，10，11，13，16，17，19，20，24，25，26，30，31，32十六个点均匀分布于分布已知水准点，经由GPS拟合程序拟合后，计算成果中得拟合高程与水准成果的互差中误差为11.820480毫米。方案二：16个起算点分布在一侧（非均匀分布）选取点位集中于右下侧，分别为1，2，3，5，9，10，11，14，18，21， 22，23，25，27，28，29十六个点。经由GPS拟合程序拟合后，计算成果中得拟合高程与水准成果的互差中误差为14.631518毫米。

曲面拟合原理与实例

问题：给定一组坐标(,,)g g g x y z ，1,2,g =…,n ，表示有n 个点。要求用以下二元多项式函数对所给的坐标进行拟合： ,1 1 111,1 11 (,)p q p q i j i j ij ij ij i j f x y a x y a x y ----==== =∑∑∑ 即 2111121312121222321 112 111231 112 11 123(,)q q q q i i i i q i i i iq p p p p q p p p pq a a y a y a y f x y a x a xy a xy a xy a x a x y a x y a x y a x a x y a x y a x y ------------++++= +++++++++++++++L L M L M L 设 1112121222221211,,q q p p pq p q a a a x y a a a x y a a a x y ???? ?????? ????????????=== ?????? ?????? ?? ?????????? x y A L L M M O M M M L 则函数又可表示为(,)T f =x y x Ay ，拟合的目标就是求出系数矩阵A 。最小二乘法：构造关于系数ij a 的多元函数： 2 112111 1 11 (,,)[(,)]()p q n n i j pq g g g g g ij g g g i j s a a f x y z a x y z ωω--=====-=-∑∑∑∑L 点（11a ，…，pq a ）是多元函数11(,,)pq s a a L 的极小点，其中g ω为权函数，默认为1，所以点（11a ，…，pq a ）必须满足方程组 0ij s a ?=? 在1g ω=的情况下，有

高程拟合

作业： 1.高程异常是如何产生的？请从实际角度谈谈如何有效地解决这一问题？答：高程异常是由地下物质及其密度分布不均匀产生的重力异常导致的。大地高与正常高之间的关系式：Hr= H84-ξ 其中ξ表示似大地水准面至椭球面间的高差，叫做高程异常。地面点的正常高Hr是地面点沿铅垂线至似大地水准面的距离。大地高是由地面点沿通过该点的椭球面法线到参考椭球面的距离，是一个几何量，不具有物理上的意义。实际上，很难获得高精度的高程异常，而GPS单点定位误差又较大，一般测区内缺少高精度的GPS基准点，GPS网平差后，很难得到高精度的大地高H84。所以很难应用上式精确的计算各GPS点正常高Hr。实际应用中解决高程异常问题，精确计算各GPS点的正常高Hr，目前主要有GPS水准高程，GPS重力高程，GPS三角高程等方法。 1 GPS水准高程目前，国内外用于GPS水准计算的各种方法主要有：绘等值线图法；解析内插法（包括曲线内插法、样条函数法和Akima法）；曲面拟和法（包括平面拟合法、多项式曲面拟合法、多面函数拟合法、曲面样条拟合法、非参数回归曲面拟和法和移动曲面法）。 1、绘等值线图法这是最早的GPS水准方法。其原理是：设在某一测区，有m个GPS点，用几何水准联测其中n个点的正常高（联测水准的点称为已知点），根据GPS观测获得的点的大地高，可以求出n个已知点的高程异常。然后，选定适合的比例尺，按n个已知点的平面坐标（平面坐标经GPS网平差后获得），展绘在图纸上，并标注上相应的高程异常，再用1～5cm的等高距，绘出测区的等高异常图。在图上内插出未联测几何水准的(m-n)个点（未联测几何水准GPS 的称为待求点），从而求出这些待求点的正常高。 2、解析内插法当GPS点布设成测线时，可应用曲线内插法，求定待求点的正常高。其原理是：根据测线上已知点的平面坐标和高程异常，用数值拟合的方法，拟合出测线方向的似大地水准面曲线，再内插出待求点的高程异常，从而求出点的正常高。

曲面拟合原理与实例

多项式函数对所给的坐标进行拟合: 构造关于系数a j 的多元函数: n 2 s( a i 1，L ， a pq ) g [ f(x g ，y g ) Z g ] g 1 点（311，…，a pq ）是多元函数s （a 11,L ,a pq ）的极小点，其中 g 为权函数，默认为1,所以点（811，…，a pq ）必须满足方程组 s 3 ij f(x,y) i 1 j 3j x y ij 1,1 1 a i i 1 j 1 i 1 j 1 i X y f (x, y) a 11 a 12y 2 a 13y L q 1 a 21x a 22xy 2 a 23xy L q 1 a 2q xy M 1 i 1 i 1 2 L i 1 q 1 a i1x y 33X y a iq X y q M p 1 p 1 p 1 2 L p 1 q a p1X a p2X y a p3X y a pq X y p,q p q 即 1 x 2 x x M x p ,y y 2 y M y q ,A a 12 L a 1q a 22 L a 2q M O M a p2 L a pq a ii a 21 M a p1 则函数又可表示为 f (x, y) x T Ay ,拟合的目标就是求出系数矩阵 A 。给定一组坐标(x g ’y g ’Z g ) , g 1,2,…,n , 表示有 n 个点。要求用以下二元 p q / i 1 j 1 g ( a j X y i 1 j 1 z g )2

在g 1的情况下，有

2 [f (X g ,y g ) Z g ] g i 2[f (X g ,y g ) i 1 j 1 Z g ]X g y g g i n 2 g i 因此可得 n n i 1 j 1 i 1 j 1 X g y g f(X g ,y g ) X g y g Z g g 1 g 1 n p q n i X g 1 y g 1 1 a X g y g 1 i 1 j X g Y g 1 Z g g 1 1 1 g 1 n p,q n i X g 1 y g 1 1 1 a X g y g i 1 j 1 X g y g z g g 1 1,1 g 1 p,q n n a i 1 j (x g y g X g y g 1 、 i 1 X y g 1z g 1,1 g 1 g 1 p,q a u (i, j) v(i, j) (i, j) (1,1),…,(p,q) 1,1 上式实际共有p q 个等式，可将这比1(1,1) L U pq (1,1) an M O M M Un(p,q) L U pq (p,q) a pq 也就是U*a=V 的形式，其中 Un(1,1) L U pq (1,1) U M O M Un(P,q) L U pq (p,q) p q 个等式写成矩阵的形式有: v(1,1) M v( p, q) an v(1,1) a M ，V M a pq v(p,q) 2[f(X g ,y g ) 叩石If")] a ij a ij x g 1 y g 1 f (x g ,y g ) x g+g z u (i, j) n (X g g 1 1 y g 1 i 1 X g y g 1 ), v(i,j) n i 1 j 1 X g y g Z g g 1

GPS静态控制网高程拟合作业指导书 8

GNSS静态控制网高程拟合作业指导书 1总则 1.1目的为了指导我院GNSS控制网高程拟合工作，统一技术标准，明确操作规程，规范技术管理，保证GNSS控制网高程拟合成果质量，特制定本技术要求。 1.2适用范围 ⑴本指导书适用于所有水环境设计院测量项目。 ⑵工作阶段为指导GNSS静态控制网高程拟合由外业到内业的一系列工作。 1.3编制依据 ①GB/T23709-2009《区域似大地水准面精化基本技术规定》; ②GB/T 12898-2009《国家三、四等水准测量规范》; ③GB/T 18314-2009《全球定位系统（GPS）测量规范》. ④SL197-97《水利水电工程测量规范》（规划设计阶段），以下简称《水电规范》。 1.4术语与定义 1.4.1似大地水准面 guasi-geoid 从地面点沿正常重力方向至正常高端点所构成的曲面。 1.4.2高程异常 Height anomaly 似大地水准面相对于地面椭球面的高度。 1.4.3高程异常控制点 control point of height anomaly GNSS/水准点 GNSS/levelling point 大地高由GNSS测定，正常高由水准测量测定的大地点。 2基本要求 2.1、GNSS静态控制网高程拟合作业流程准备工作→技术设计（含高程联测）→网图设计→选点埋石→GNSS网静态观测→高程联测→高程网平常→数据处理→平面网平差计算（高程拟合）→质量检查与自检报告→技术报告→成果整理与提交 2.2、GNSS高程拟合作业流程图 2.3、GNSS高程转换的基本原理

2.3.1、在一个测区内若有若干个既进行了GNSS 测量又联测了水准高程的GNSS 点，这样的点称为水准重合点，后面简称已知点。 2.3.2、利用大地高和水准高之间的关系，推算出各水准重合点上的高程异常，利用这些离散点上的异常值，可以拟合出测区所在局部区域的似大地水准面，进而可以内插出未知点上的高程异常，实现椭球高向正常高的转换。 2.3.3、GNSS 高程转换的基本公式 ⑴首先进行GNSS 相对定位测量，由GNSS 网三维平差可以求得各埋石点的GNSS 大地高Hi ，各点正常高可以通过下式求得： ξ i i i H h - = 式（2-1）由于GNSS 网中起算点的大地高精度不高，控制网平差时一般只取一个点的大地高作为起算数据。由此求得的网中各点的大地高与实际大地高之间偏移了一个平移量△H O ，此时(2-1)式变为 (2-2) 如果在GNSS 网中某点Pk 处联测了水准，则 (2-3) 由式（2-2）．式(2-3)得 (2-4) 即 (2-5) 从(2-5)式可以看出，虽然由GNSS 确定的大地高(H j 、H k )的精度不高，但由于GNSS 基线向量的精度很高，因而各点之间的大地高差(H i -H k )仍可达到很高的精度，所以在GNSS 高程转换时，最好采用高程异常差建立模型。根据重合点的高程异常与水平位置，可以建立测区的高程异常（差）面，从而采用内插法可以获得网中其它各点的高程异常，根据网中待定点的大地高和上面所求得的高程异常便可求得工程需要的正常高，实现高程转换。 3、GNSS 网设计、观测与平差 3.1、GNSS 网设计表3.1 水利工程GNSS 网的精度标准等级项目二三四五等一级五等二级固定误差/mm ≤10 ≤10 ≤10 ≤10 ≤10 比例误差系数 ≤5 ≤10 ≤20 ≤40 ≤80 相邻点最小距离/km 2～4 1～2 0.5～1 0.2～0.5 0.1～0.2 相邻点平均距离/km 8～13 4～8 2～4 1～2 0.5～1 最弱相邻点边长相对中误差 1/150000 1/80000 1/40000 1/20000 1/10000

基于曲面拟合的图像分割算法

基于曲面拟合的图像分割算法作者：禹建东孔月萍来源：《现代电子技术》2008年第22期摘要：传统分割方法，在光照不均匀情况下，很难得到理想的分割结果。针对这种情况，提出一种基于曲面拟合的阈值曲面分割方法。首先利用偏离项和光顺项构造拟合曲面方程，然后使用在统一的大光顺项因子条件下求解的拟合结果，来构造自适应的偏离因子与光顺因子，最后利用这些自适应因子第二次精确拟合阈值曲面。实验表明，该方法对于照度不均匀的图像，分割结果明显优于传统全局阈值法。关键词：图像分割；拟合；阈值曲面；B样条中图分类号：TP391文献标识码：B 文章编号：1004-373X(2008)22-106-02 Image Segmentation Based on Surface Fitting YU Jiandong,KONG Yueping (School of Information and Control,Xi′an University of Architecture,Xi′an,710055,Chi na) Abstract：The traditional segmentation method cannot get the ideal result under the non-uniform illumination.With regard to this situation,a segmentation method using threshold surface based on surface fitting is proposed.Firstly,using deflection item and fairing item to construct the equation of fitting surface.Secondly,In order to construct the adaptive weight factors,using the result of fitting surface with more heavy fairing item,then fit the threshold surface exactly again.The experimental result shows that the method has better than traditional method. Keywords：image segmentation;fitting;threshold surface;B-spline 1 引言图像分割是计算机视觉领域中极为重要的一环，是实现图像内容识别之前首先要完成的工作。分割效果的好坏，决定了识别正确率的高低。传统的全局阈值法，只有在对双峰特征的图像时才有较好的效果。而当图像中存在照度不均匀、或者背景灰度变化等情况，则往往达不到令人满意的分割结果。因此自适应阈值分割技术应运而生，它主要利用图像的局部特征，根据不同的区域自适应地选取相应阈值，构造一个用于分割的阈值曲面。其中分割效果比较好的如近些年出现的基于变分的图像分割，它利用图

RTK高程拟合

工程之星3.0 特色功能之一：控制点测量介绍 S730手簿蓝牙传输文件过程 RTK测量高程精度简析 2011-05-26 13:26:55| 分类：RTK测量资料| 标签：|字号大中小订阅石家庄南方测绘导航产品部郭晓辉使用RTK做地形图测量，既能快速的获得平面坐标又能快速的获得高程，大家都很容易接受，可是当谈论到使用RTK 是否可以做水准测量时，不少朋友都在心里打了一个问号。到底RTK 测得的高程和水准测量差多少呢？能不能满足工程的要求。其实这方面的问题已经被专家论证了多次，答案是在严格控制及选用合理的作业方法下，RTK 测量高程可以满足四等水准测量及等外的水准测量。毫无疑问，使用RTK 进行水准测量将会大大降低工作强度，同时提高作业效率。下面就介绍一下，如何使用RTK达到如上所述的效果。首先分析下GPS测得高程和水准测量求高程的区别，GPS 测量求得的原始坐标是WGS-84坐标（B,L,H）大地纬度，大地精度，大地高。而我国水准测量是采用1985国家高程基准，以似大地水准面为起算面，最后是以正常高作为使用的高程。因为测量原理不同，两种测量的起算面不同，所以两种高程值之间存在高程异常，即大地高= 正常高+高程异常。所以如果使用GPS要达到水准测量要求的正常高的值，必须要求提高得的大地高和高程异常值的精度。大地高的精度如南方灵锐S86RTK的精度指标垂直精度±2cm+1ppm ，静态，快速静态高程精度±5mm+1ppm，而精确的求出高程异常就是关键所在。南方GPS，RTK 用高程拟合的方法精确求得高程异常，从而可以实时的得到控制范围内的正常高。 GPS 水准高程拟合方法是: 在GPS 网中联测一些水准点, 利用这些点上的正常高和大地高求出它们的高程异常值, 再根据这些点上的高程异常值与坐标的关系,用最小二乘的方法拟合出测区的似大地水准面,利用拟合出的似大地水准面,内插出其他GPS 点的高程异常, 从而求出各个未知点的正常高。用于GPS 水准拟合的数学模型很多, 不同的数学模型对不同地形条件具有不同的拟合精度, 因此GPS 水准拟合模型拟合精度的探讨一直是GPS 应用研究领域的热点问题。其中多项式就是GPS 水准拟合模型的一种,其模型可表述为 ζ= f ( x , y ) + ε 当GPS 点布设成网状时,一般采用曲面拟合的方法。设测站点的高程异常ζ与坐标之间存在以下函数关系ζ i = f ( xi , y i ) + ε i其中, f ( xi , y i ) 为ζ的趋势值, ε i 为误差。选用空间曲面函数 f ( x i , yi ) = a0 + a1x i + a2y i + a3x2i + a4x iyi + a5 y2i + a6 x3i + a7 x2iy i + a8x iy2i + a9y3i ( 4)进行拟合,式中ai 为待定参数。在已知点个数大于等于参数个数求出参数ai ,进而求出测区内任意点的高程异常。根据测区的不同情况,也可以选用不同的参数进行拟合。选用的参数不同,拟合出的曲面的形式也不相同。 1多项式拟合模型分型

二次曲面的一般理论

第六章二次曲面的一般理论教学目的 : 本章讨论了一般二次曲面的渐近方向、中心、切线、切平面、径面奇向、主径面与主方向等重要概念 ,从不同角度对二次曲面进行了分类 . 研究了二次曲面的几何性质 , 并通过坐标变换和不变量、半不变量两种形式化二次曲面的一般方程为规范方程 , 对二次曲面进行了分类和判定 , 是二次曲面理论的推广和扩充 . 教学重难点 : 通过坐标变换和运用不变量、半不变量化二次曲面的一般方程为规范方程 , 既是重点又是难点 . 基本概念二次曲面 : 在空间 , 由三元二次方程 2 2 2 a 11x a 22 y a 33z 2a 12 xy 2a 13 xz 2a 23 yz 2a 14 x 2a 24 y 2a 34z a 44 0 (1) 所表示的曲面 . 虚元素 :空间中，有序三复数组 (x,y,z) 叫做空间复点的坐标，如果三坐标全是实数，那么它对应的点是实点，否则叫做虚点二次曲面的一些记号 F(x,y,z) F 1(x,y,z) a 11x a 12y a 13z a 14 F 2(x,y,z) a 12x a 23y a 23z a 24 F 3( x, y, z) a 13x a 23y a 33z a 34 F 4 (x,y,z) a 14x a 24y a 34z a 44 2 2 2 (x, y,z) a 11x 2 a 22 y 2 a 33z 2 2a 12 xy 2a 13 xz 2a 23 yz 1 (x,y,z) a 11x a 12 y a 13z 2 (x,y,z) a 12 x a 22 y a 23z 2 a 11 x 22 a 22 y a 33 z 2a 12 xy 2a 13 xz 2a 23 yz 2a 14 x 2a 24 y 2a 34 z a 44

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

高程拟合的方法和原理(二次曲面拟合代码) By Kiseigo kiseigo https://www.doczj.com/doc/bd5410669.html,/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

二次曲面的一般理论

第六章二次曲面的一般理论教学目的: 本章讨论了一般二次曲面的渐近方向、中心、切线、切平面、径面奇向、主径面与主方向等重要概念,从不同角度对二次曲面进行了分类. 研究了二次曲面的几何性质,并通过坐标变换和不变量、半不变量两种形式,化二次曲面的一般方程为规范方程,对二次曲面进行了分类和判定,是二次曲面理论的推广和扩充. 教学重难点: 通过坐标变换和运用不变量、半不变量化二次曲面的一般方程为规范方程,既是重点又是难点. 基本概念二次曲面: 在空间,由三元二次方程 022222244342414231312233222211=+++++++++a z a y a x a yz a xz a xy a z a y a x a （1）所表示的曲面. 虚元素:空间中，有序三复数组),,(z y x 叫做空间复点的坐标，如果三坐标全是实数，那么它对应的点是实点，否则叫做虚点二次曲面的一些记号 ≡ ),,(z y x F 44 342414231312233222211222222a z a y a x a yz a xz a xy a z a y a x a +++++++++ 141312111),,(a z a y a x a z y x F +++≡ 242323122),,(a z a y a x a z y x F +++≡ 343323133),,(a z a y a x a z y x F +++≡ 443424144),,(a z a y a x a z y x F +++≡ yz a xz a xy a z a y a x a z y x 231312233222211222),,(+++++≡Φ z a y a x a z y x 1312111),,(++≡Φ z a y a x a z y x 2322122),,(++≡Φ

移动二次曲面拟合内插DEM程序

似水无痕一、二次曲面移动拟合法内插DEM的原理 DEM内插就是根据参考点上的高程求出其他待定点上的高程,在数学上属于插值问题。任意一种内插方法都是基于邻近数据点之间存在很大的相关性,从而由邻近点的数据内插出待定点上的数据。移动曲面拟合法内插,是以每一待定点为中心,定义一个局部函数去拟合周围的数据点。该方法十分灵活,一般情况精度较高,计算相对简单,不需很大计算机内存,其过程如下: (1)根据实际内插要求,解算待定点P的平面坐标( xP , yp )。 (2)为了选取邻近的数据点,以待定点P为圆心,以R为半径作圆(如图1所示) ,凡落在圆内的数据点即被选用。在二次曲面内插时,考虑到计算方便,将坐标原点移至该DEM格网点P ( xP , yp ) 由于二次曲面系数个数为6,要求选用的数据点个数n > 6。当数据点i ( xi , yi )到待定点P ( xP , yp )的距离di = sqr(x i2 - .y i2 ) i < R时,该点即被选用。若选择的点数不够时,则应增大R的数值,直至数据点的个数n 满足要求。(3)选择二次曲面Z =Ax2 +B xy +Cy2 +Dx + Ey +F作为拟合面,则对应点的误差方程为 vi = Ax2 +B xy +Cy2 +Dx + Ey +F - Z i 由n个数据点列出的误差方程为v = MX - Z X = (M T PM) -1M T PZ 由于坐标原点移至该DEM格网点P ( xP , yp ) ,所以系数F就是待定点的内插高程值ZP 二、程序

采用平台：https://www.doczj.com/doc/bd5410669.html,,access数据库（存储已知点）数据：10个点程序代码： Imports System.Math Public Class Form1 Dim conn As OleDb.OleDbConnection = New OleDb.OleDbConnection Dim cmd As OleDb.OleDbCommand = New OleDb.OleDbCommand Dim adapter As OleDb.OleDbDataAdapter = New OleDb.OleDbDataAdapter Dim datareader As OleDb.OleDbDataReader '数据库连接 Dim xe As Double, ye As Double, d(9) As Double, X(9) As Double, Y(9) As Double, Z(9, 0) As Double Dim M(9, 5) As Double, P(9, 9) As Double Dim Xz(5, 0) As Double Private Sub Button1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button1.Click

GPS控制网高程拟合

GPS控制网高程拟合【摘要】通过对沁河防汛工程D级GPS网的高程拟合精度分析，探讨GPS高程拟合成果的精度与起算点分布、起算成果精度、高程拟合数学模型、GPS数据处理软件的关系。 GPS network of Qinhe flood control projects D elevation fitting accuracy, explore the accuracy of the GPS elevation fitting the results with the starting point of distribution, the date the results of precision, the elevation fitting a mathematical model, the relationship of the GPS data processing software. 【关键词】GPS 高程异常值中误差曲面拟合EGM96大地水准面模型前言全球定位系统（Global Positioning System-GPS）是美国从本世纪70年代开始研制，于1994年全面建成，具有在海、陆、空进行全方位实时三维导航与定位能力的新一代卫星导航与定位系统。经近10年我国测绘等部门的使用表明，GPS以全天候、高精度、高效率等显著特点，赢得广大测绘工作者的信赖，并成功地应用于大地测量、工程测量、航空摄影测量、运载工具导航和管制、地壳运动监测、工程变形监测、资源勘察、地球动力学等测绘学科，给测绘领域带来一场深刻的技术革命。目前，大多数的城市首级控制网均采用GPS测量，而其中的高程控制主要采用传统的几何水准测量方法建立高精度的水准网。GPS高程测量却常常被忽视，认为其精度不太可靠。因此，为探讨GPS测量高程拟合成果的精度与起算点分布、起算成果精度、高程拟合数学模型、GPS数据处理软件的关系，我局结合沁河防汛工程D级测量GPS高程拟合的工作，对GPS拟合高程的精度进行了探讨，对于平坦地区以供测量GPS用户参考。 1 GPS网高程拟合的技术要求 1.1 GPS高程拟合成果外部检核 1.1.1 首先对D级GPS网中的所有点联测四等水准或三角高程，选用其中部分点作为GPS高程拟合的

GPS高程拟合方法的比较分析

GPS 高程拟合法的比较分析 (机械工业勘察设计研究院测量公司) 摘要：工程中需要把GPS 高程测量的大地高转换为正常高。通常的做法是采用拟合法建立研究区域的似大地水准面。本文介绍了两种不同的拟合方法：二次曲面拟合法、多面函数拟合法。并结合某区域一定数量已知GPS 高程异常点来内插和外推研究区域内的任一点的高程异常。通过比较发现多面函数拟合法拟合的精度要比二次曲面拟合的精度高。关键词：高程转换；二次曲面拟合法；多面函数拟合法 The elevation of GPS fitting to the comparison and analysis （Machinery industry survey and design institute of measuring company ） Abstract: GPS height measurement of the earth should be converted to normal high in engineering. It is usually to establish the quasi-geoid of the research area by the fitting method. This article introduces two different fitting methods: quadratic surface fitting and multiple-surface function fitting. Combined with a certain number of a region known GPS elevation anomaly points to the interpolation and extrapolation of the height anomaly at any point within the study area. By comparison, the multiple-surface function fitting to the precision is higher than the quadratic surface fitting. Key words ：Elevation conversion; Quadratic surface fitting; Multiple-surface function fitting 1.引言传统的几何水准测量虽然精度高，但耗时长、耗费多、工作效率低。GPS 由于自身测量精度高、速度快、工作效率高等优点被广泛应用于高程测量。GPS 测量的高程坐标是在WGS-84坐标系下的大地高[1]，大地高是地面一点沿参考椭球面的法线到参考椭球面的距离，用符号H 表示。实际应用中需要把GPS 测得的大地高转换为正常高，正常高是地面点到通过该点的铅垂线与似大地水准面的交点的距离，用符号r H 表示。似大地水准面到参考椭球面之间的距离称为高程异常，用符号ζ表示。因此大地高与正常高之间的关系为： r H H ζ=- （1）由于我国采用的高程系统是相对于似大地水准面的正常高，因此如何进行GPS 高程转换成为当前研究的热点问题。拟合法是GPS 高程转换中比较常用的方法，主要的拟合模型

最小二乘曲面拟合

%可用样条曲面拟合，最好对原数据整理一下，拟合的代码如下：x0=2.2:0.1:7;y0=10:5:30; z0 =[ 0.0121 0.0118 0.0129 0.1098 0.0103 0.0116 0.0116 0.0124 0.1007 0.0111 0.0110 0.0113 0.0120 0.0914 0.0119 0.0105 0.0111 0.0116 0.0820 0.0128 0.0099 0.0109 0.0112 0.0726 0.0136 0.0094 0.0107 0.0108 0.0635 0.0144 0.0090 0.0105 0.0105 0.0547 0.0151 0.0085 0.0104 0.0101 0.0465 0.0158 0.0081 0.0102 0.0098 0.0391 0.0164 0.0078 0.0101 0.0096 0.0325 0.0170 0.0075 0.0100 0.0094 0.0270 0.0174 0.0073 0.0099 0.0092 0.0228 0.0177 0.0072 0.0099 0.0091 0.0200 0.0179 0.0071 0.0098 0.0091 0.0187 0.0180 0.0071 0.0098 0.0091 0.0183 0.0180 0.0071 0.0098 0.0091 0.0179 0.0180 0.0071 0.0098 0.0091 0.0176 0.0180 0.0072 0.0099 0.0091 0.0172 0.0180 0.0072 0.0099 0.0091 0.0169 0.0180 0.0072 0.0099 0.0091 0.0165 0.0180 0.0072 0.0099 0.0091 0.0162 0.0180 0.0072 0.0099 0.0091 0.0159 0.0180 0.0073 0.0099 0.0091 0.0156 0.0179 0.0073 0.0100 0.0092 0.0154 0.0179 0.0074 0.0100 0.0092 0.0151 0.0178 0.0075 0.0101 0.0093 0.0149 0.0178 0.0076 0.0101 0.0093 0.0147 0.0177 0.0077 0.0102 0.0094 0.0144 0.0177 0.0078 0.0102 0.0095 0.0142 0.0176 0.0079 0.0103 0.0095 0.0140 0.0175 0.0081 0.0104 0.0096 0.0139 0.0174 0.0082 0.0105 0.0097 0.0137 0.0173 0.0084 0.0106 0.0099 0.0135 0.0171 0.0086 0.0107 0.0100 0.0134 0.0170 0.0089 0.0108 0.0101 0.0133 0.0168 0.0091 0.0109 0.0103 0.0131 0.0166 0.0094 0.0111 0.0105 0.0130 0.0164 0.0097 0.0112 0.0107 0.0129 0.0162 0.0100 0.0114 0.0109 0.0128 0.0160 0.0104 0.0115 0.0111 0.0128 0.0157 0.0108 0.0117 0.0114 0.0127 0.0155 0.0112 0.0119 0.0116 0.0126 0.0152