python 趋势面模型分析实例
已知五个气象站,这五个站围绕着未知的0号站。下图为各个站点的x、y坐标及及其已知值
第一行为0号点x,y坐标 z值未知暂设为0
而后依次为1,2,3,4,5号点xyz值
对未知点0进行插值。最小二乘法通常用于计算方程:
中的待定系数b0,b1,b2。因此。第一步是建立如下三个正规方程:
以上方程改写成如下矩阵:
将左边第一个逆矩阵与右边的矩阵相乘,计算出系数b
而后利用公式:
从而求出未知值
代码如下:
print('\n###山人青冥的趋势面分析小程序###\n')
import numpy as np
from numpy.linalg import *
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']
data=np.loadtxt(open('站点坐标.csv'),delimiter=",",skiprows=0) #读取坐标
print('\n站点坐标:\n',data)
x=data[1:,0] #提取x第一列
y=data[1:,1] #提取y第二列
z=data[1:,2] #提取z第二列
sx=sum(x);sy=sum(y);sz=sum(z)
sxy=sum(x*y);sxx=sum(x*x);syy=sum(y*y)
sxz=sum(x*z);syz=sum(y*z)
p1=[len(x),sx,sy,sx,sxx,sxy,sy,sxy,syy]
p2=[sz,sxz,syz]
p1=np.array(p1).reshape(3,3)
p2=np.array(p2).reshape(3,1)
b=np.dot(inv(p1),p2)
z0=b[0]+b[1]*data[0,0]+b[2]*data[0,1]
print('\n0号点未知值计算结果:\n',z0)
data[0,2]=z0
#站点三维坐标图 0号点高程未知,暂设为0
x=data[:,0] #提取x第一列
y=data[:,1] #提取y第二列
z=data[:,2] #提取z第二列
ax= plt.figure().add_subplot(111,projection='3d')
ax.scatter(x,y,z,c='b',marker='D')
ax.set_title('0号站点及周围站点三维坐标图',fontsize=12,color='teal')
for i in range(0,len(x)): #添加点标签
ax.text(x[i],y[i],z[i],i)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()
#画图结束
运行结果:
上述示例为简单模型,大多数自然现象的分布通常比由一阶趋势面生成的倾斜更复杂。因而,需要更高阶的趋势面模型来拟合更复杂的表面。
参考文献:
[1] (美)张康聪 著;陈健飞等译. 地理信息系统导论(第三版). 北京:清华大学出版社, 2009.04.