趋势面模型分析实例（python）

python 趋势面模型分析实例

已知五个气象站，这五个站围绕着未知的0号站。下图为各个站点的x、y坐标及及其已知值

第一行为0号点x，y坐标 z值未知暂设为0
而后依次为1，2，3，4，5号点xyz值

 对未知点0进行插值。最小二乘法通常用于计算方程：

中的待定系数b0，b1，b2。因此。第一步是建立如下三个正规方程：

以上方程改写成如下矩阵：

将左边第一个逆矩阵与右边的矩阵相乘，计算出系数b
而后利用公式：

从而求出未知值

具体代码如下：

import numpy as np
from numpy.linalg import *
import matplotlib.pyplot as plt
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[2]     http://t.csdn.cn/qgJ8N