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

《移动二次曲面拟合内插DEM程序》由会员分享，可在线阅读，更多相关《移动二次曲面拟合内插DEM程序（13页珍藏版）》请在装配图网上搜索。

1、似水无痕 一、二次曲面移动拟合法内插 DEM 的原理DEM 内插就是根据参考点上的高程求出其他待定点上的高程,在数学上属于插值问题。任意一种内插方法都是基于邻近数据点之间存在很大的相关性,从而由邻近点的数据内插出待定点上的数据。移动曲面拟合法内插,是以每一待定点为中心 ,定义一个局部函数去拟合周围的数据点。 该方法十分灵活 ,一般情况精度较高 , 计算相对简单 ,不需很大计算机内存 ,其过程如下 :(1)根据实际内插要求 ,解算待定点 P 的平面坐标 ( xP , yp ) 。为了选取邻近的数据点，以待定点P为圆心，以R为半径作圆(如图1所示)，凡落 在圆内的数据点即被选用。在二次曲面内插时

2、,考虑到计算方便 ,将坐标原点移至 该 DEM 格网点 P ( xP , yp )由于二次曲面系数个数为 6,要求选用的数据点个数 n 6。当数据点 i ( xi , yi ) 到待定点 P ( xP , yp )的距离 di = sqr(xi2 - .yi2 ) i R 时, 该点即被选用。若选择 的点数不够时 ,则应增大 R 的数值,直至数据点的个数 n满足要求。 (3)选择二次曲面 Z =Ax2 +B xy +Cy2 +Dx + Ey +F 作为拟合面 ,则对应 点的误差方程为22vi = Ax +B xy +Cy +Dx + Ey +F - Zi由 n 个数据点列出的误差方程为 v =

3、 MX - ZT - 1 TX = (MTPM) - 1MTPZ由于坐标原点移至该DEM格网点P ( xP , yp ),所以系数F就是待定点的内插高 程值 ZP二、程序采用平台：vs2008-vb. net,acces数据库(存储已知点)数据： 10 个点程序代码：Imports System.MathPublic Class Form1Dim conn As OleDb.OleDbConnection = New OleDb.OleDbConnectionDim cmd As OleDb.OleDbCommand = New OleDb.OleDbCommandDim adapter As

4、 OleDb.OleDbDataAdapter = New OleDb.OleDbDataAdapterDim datareader As OleDb.OleDbDataReader数据库连接Dim xe As Double, ye As Double, d(9) As Double, X(9) As Double, Y(9) AsDouble, Z(9, 0) As DoubleDim M(9, 5) As Double, P(9, 9) As DoubleDim Xz(5, 0) As DoublePrivate Sub Button1_Click(ByVal sender As Syst

5、em.Object, ByVal e AsSystem.EventArgs) Handles Button1.ClickDim i As Integerxe = 110ye = 110 待定点坐标DataSet1.Tables.Clear()DataGridView1.Columns.Clear()conn.ConnectionString = Provider=Microsoft.Jet.OLEDB.4.0;Data Source= &Application.StartupPath & DEM.mdbcmd.CommandText = select 点号 ,X,Y,Z from DEM cm

6、d.Connection = conn conn.Open()OleDbDataAdapter1.SelectCommand = cmdOleDbDataAdapter1.Fill(DataSet1, DEM)DataGridView1.DataSource = DataSet1.Tables.Item(DEM) conn.Close()连接数据库，把数据读到 DatagridviewFor i = 0 To 9X(i) = DataGridView1.Rows.Item(i).Cells.Item(X).ValueY(i) = DataGridView1.Rows.Item(i).Cells

7、.Item(Y).ValueZ(i, 0) = DataGridView1.Rows.Item(i).Cells.Item(Z).ValueM(i, 0)= (X(i) - xe) * (X(i) - xe)M(i, 1)= (X(i) - xe) * (Y(i) - ye)M(i, 2)= (Y(i) - ye) * (Y(i) - ye)M(i, 3)= X(i) - xeM(i, 4)= Y(i) - xeM(i, 5)= 1P(i, i) =1 / Sqrt(Abs(X(i) - xe) * (X(i) - xe) + (Y(i) - ye) * (Y(i) - ye)Next对矩阵赋

8、初值End SubPrivate Sub Button2_Click(ByVal sender As System.Object, ByVal e AsSystem.EventArgs) Handles Button2.ClickTextBox1.Text = Dim i As IntegerXz = JZCF(JZCF(JZCF(NJZ(JZCF(JZCF(DZ(M), P), M), DZ(M), P), Z)TextBox1.Text = TextBox1.Text & 待定点的坐标： X=110,Y=110 & vbCrLf & vbCrLfTextBox1.Text = TextBo

9、x1.Text & 误差方程为：V=AXiA2+BXiYi+CYP2+DXi+EYi+F-Zi & vbCrLf & vbCrLfTextBox1.Text = TextBox1.Text & 由 X=(M(T)PM)M(T)PZ & vbCrLfTextBox1.Text = TextBox1.Text & 其中 X=(A B C D E F )(T) Z=(Z1 Z2 Z3 、Z n)(T) & vbCrLfTextBoxl.Text = TextBoxl.Text & M=(XiA2 XiYi YiA2 Xi Yi 1) &vbCrLf & vbCrLfTextBox1.Text = T

10、extBox1.Text & 求得 X=(For i = 0 To 5TextBox1.Text = TextBox1.Text & Format(Xz(i, 0), 0.000) & NextTextBox1.Text = TextBox1.Text & ) & vbCrLf & vbCrLfTextBox1.Text = TextBox1.Text & 由于坐标原点移至待定点上，因此待定 点的高程 Z=F & vbCrLfTextBox1.Text = TextBox1.Text & 待定点高程 Z= & Format(Xz(5, 0), 0.000)End SubPublic Funct

11、ion NJZ(ByVal q(,) As Double) As Double(,)Dim BS(,) As Double, NJ(,) As Double, Qz As Double, h As Integer, l As Integer, n As Integern = q.GetUpperBound(0)ReDim BS(n, n)ReDim NJ(n, n)BS = BSJZ(q) 伴随矩阵Qz = HLei(q)If Qz 0 ThenFor h = 0 To nFor l = 0 To nNJ(h, l) = (1 / Qz) * BS(h, l) 求解矩阵的逆矩阵NextNext

12、NJZ = NJElseMessageBox.Show(矩阵不可逆”)NJZ = NothingEnd IfEnd FunctionPublic Function BSJZ(ByVal q(,) As Double) As Double(,) 求解矩阵的伴随矩阵Dim i As Integer, j As Integer, k As Integer, n As Integer, BS(,) As Doublei = 0j = 0k = 0n = q.GetUpperBound(0)ReDim BS(n, n)For i = 0 To nFor j = 0 To nBS(j, i) = YZS(

13、q, i, j)NextNextBSJZ = BSEnd Function求解矩阵的余子式的值Public Function YZS(ByVal q(,) As Double, ByVal h As Integer, ByVal l As Integer) As Double 求解矩阵的余子式的值Dim i As Integer, j As Integer, x As Integer, y As Integer, n As Integer, YZ(,) As Doublei = 0j = 0n = q.GetUpperBound(0)ReDim YZ(n - 1, n - 1)For i =

14、0 To nIf i h ThenFor j = 0 To nIf j l Then 求得 h,l 位置的余子式YZ(x, y) = q(i, j)y += 1End IfNext x += 1End IfNextYZS = (-1) A (h + l + 2) * HLei(YZ)余子式的值End Function求解 n 阶行列式的值Public Function HLei(ByVal q(,) As Double) As Double 求解 n 阶行列式的值Dim i As Integer, j As Integer, k As Integer, n As IntegerIf q.Get

15、UpperBound(0) 1 Then 判断是否为 2*2 矩阵，若不是则继续求其余子式n = q.GetUpperBound(0)Dim p(n - 1, n - 1) As Double 余子式方阵For k = 0 To nDim h As IntegerDim l As Integer 余子式的行列号h = 0For i = 1 To nl = 0For j = 0 To nIf j k Then 按第一列展开求余子式 p(h, l) = q(i, j)l += 1End IfNext h += 1NextHLei = HLei + (-1) A (k + 2) * q(0, k)

16、* HLei(p)递归方法求行列式的值NextElseHLei = q(0, 0) * q(1, 1) - q(0, 1) * q(1, 0) 二阶行列式求值End IfEnd FunctionPublic Function DZ(ByVal q(,) As Double) As Double(,) 求解矩阵的转置Dim i As Integer, j As Integer, n As Integer, m As Integer, DZJZ(,) As Doublei = 0j = 0n = q.GetUpperBound(0)m = q.GetUpperBound(1)ReDim DZJZ(

17、m, n)For i = 0 To nFor j = 0 To mDZJZ(j, i) = q(i, j)NextNextDZ = DZJZEnd FunctionPublic Function JZCF(ByVal a(,) As Double, ByVal b(,) As Double) As Double(,) 矩阵乘法Dim ha As Integer, la As Integer, hb As Integer, lb As Integer, c(,) As DoubleDim i As Integer, j As Integer, k As Integerha = a.GetUpperBound(0)la = a.GetUpperBound(1)hb = b.GetUpperBound(0)lb = b.GetUpperBound(1)If la = hb ThenReDim c(ha, lb)For i = 0 To haFor j = 0 To lbFor k = 0 To lac(i, j) = c(i, j) + a(i, k) * b(k, j)NextNextNextJZCF = cElseMessageBox.Show(这两个矩阵不能相乘！”)JZCF = NothingEnd IfEnd FunctionEnd Class