似水无痕
一、二次曲面移动拟合法内插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
二、程序
采用平台：,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
Dim i As Integer
xe = 110
ye = 110 '待定点坐标
DataSet1.Tables.Clear()
DataGridView1.Columns.Clear()
conn.ConnectionString = "Provider=Microsoft.Jet.OLEDB.4.0;Data Source=" & Application.StartupPath & "\DEM.mdb"
mandText = "select 点号,X,Y,Z from DEM "
cmd.Connection = conn
conn.Open()
OleDbDataAdapter1.SelectCommand = cmd
OleDbDataAdapter1.Fill(DataSet1, "DEM")
DataGridView1.DataSource = DataSet1.Tables.Item("DEM")
conn.Close()
'连接数据库，把数据读到Datagridview
For i = 0 To 9
X(i) = DataGridView1.Rows.Item(i).Cells.Item("X").Value
Y(i) = DataGridView1.Rows.Item(i).Cells.Item("Y").Value
Z(i, 0) = DataGridView1.Rows.Item(i).Cells.Item("Z").Value
M(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) - xe
M(i, 4) = Y(i) - xe
M(i, 5) = 1
P(i, i) = 1 / Sqrt(Abs((X(i) - xe) * (X(i) - xe) + (Y(i) - ye) * (Y(i) - ye)))
Next
'对矩阵赋初值
End Sub
Private Sub Button2_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button2.Click
TextBox1.Text = ""
Dim i As Integer
Xz = JZCF(JZCF(JZCF(NJZ(JZCF(JZCF(DZ(M), P), M)), DZ(M)), P), Z)
