DEM移动拟合法逐点内插算法作业（python）

     ps1:用python写的第一个程序，怎么也要纪念一下。

     ps2:写完后不是很想整理，代码可能有些混乱。

     ps3:该作业要求中待求点固定且只有一个，附近点已给出且固定，不包含动态圆半径法等方法搜索临近点过程，所以程序不存在较强泛用性，但是或许可以给完整的逐点内插方法提供部分过程思路。

1.数据录入和预处理

      由于在移动曲面法逐点内插中，需要将坐标系平移至待求点P，所以需要对提供的坐标数据进行处理，其X、Y值改写为以待求点为中心的相对值，Z值保持不变。数据源此处按下不表了，改写后的坐标数值如下表。

点号	X	Y	Z
1	-8	0	15
2	-1	3	18
3	-5	5	19
4	-7	-7	17
5	-2	-5	21
6	-5	-2	15
7	5	-6	20
8	8	-2	15
9	6	3	17
10	3	8	22

     引入处理矩阵的numpy库（这个不是python自带的，没有的小伙伴要搜搜下载教程），录入list格式的x,y坐标，以array格式录入z坐标。

import numpy as np#引入numpy库

x = [0,-8,-1,-5,-7,-2,-5,5,8,6,3] #x和y属于List，只需要元素参与计算，不进行矩阵整体计算
y = [0,0,3,5,-7,-5,-2,-6,-2,3,8]  #x0,y0是起始点
z = np.array([[15],[18],[19],[17],[21],[15],[20],[15],[17],[22]]) #z要进行矩阵运算，所以是待转化为矩阵的数组
Z = np.matrix(z)#将数列z转化为矩阵Z

2.数据处理与计算

    根据x,y坐标，求出方程Z=AX^2+BXY+CY^2+DX+EY+F中由X,Y组成的矩阵m，使用numpy中的函数np.transpose()对其进行转置处理。

    计算每个点与待求点之间的距离d，根据公式p=1/d^2求出矩阵p。

    检查法方程中所需要素是否均为矩阵格式，若不是，使用np.matrix方法将其转化为矩阵。ps:数据是list,arry还是matrix的影响挺大，比如matrix的运算只跟matrix玩，matrix不好提取单个数据但是array可以，这里一定要注意。

    使用逆矩阵方法对MTPM进行求逆，得到（MTPM）^-1，再将各元素带入法方程W=（MTPM）^-1(MTPZ)中进行计算，得到六个系数值所组成的矩阵W。ps:此处为了求逆矩阵尝试了np.linalg.inv()方法，但是好像不能运行就放弃了。

m = np.array([
             [x[1]**2,x[1]*y[1],y[1]**2,x[1],y[1],1],
             [x[2]**2,x[2]*y[2],y[2]**2,x[2],y[2],1],
             [x[3]**2,x[3]*y[3],y[3]**2,x[3],y[3],1],
             [x[4]**2,x[4]*y[4],y[4]**2,x[4],y[4],1],
             [x[5]**2,x[5]*y[5],y[5]**2,x[5],y[5],1],
             [x[6]**2,x[6]*y[6],y[6]**2,x[6],y[6],1],
             [x[7]**2,x[7]*y[7],y[7]**2,x[7],y[7],1],
             [x[8]**2,x[8]*y[8],y[8]**2,x[8],y[8],1],
             [x[9]**2,x[9]*y[9],y[9]**2,x[9],y[9],1],
             [x[10]**2,x[10]*y[10],y[10]**2,x[10],y[10],1]
             ])#注意0给待求点了，这里从1开始
M = np.matrix(m)#转为矩阵
mt = np.transpose(m)#M的转置矩阵
MT = np.matrix(mt)#转化为矩阵
# Mm=np.linalg.inv(m) #求逆矩阵达咩的方法，还是得转化为矩阵后再用.I方法

p = np.array([#权重矩阵1/d^2
             [1/((x[0]-x[1])**2+(y[0]-y[1])**2),0,0,0,0,0,0,0,0,0],
             [0,1/((x[0]-x[2])**2+(y[0]-y[2])**2),0,0,0,0,0,0,0,0],
             [0,0,1/((x[0]-x[3])**2+(y[0]-y[3])**2),0,0,0,0,0,0,0],
             [0,0,0,1/((x[0]-x[4])**2+(y[0]-y[4])**2),0,0,0,0,0,0],
             [0,0,0,0,1/((x[0]-x[5])**2+(y[0]-y[5])**2),0,0,0,0,0],
             [0,0,0,0,0,1/((x[0]-x[6])**2+(y[0]-y[6])**2),0,0,0,0],
             [0,0,0,0,0,0,1/((x[0]-x[7])**2+(y[0]-y[7])**2),0,0,0],
             [0,0,0,0,0,0,0,1/((x[0]-x[8])**2+(y[0]-y[8])**2),0,0],
             [0,0,0,0,0,0,0,0,1/((x[0]-x[9])**2+(y[0]-y[9])**2),0],
             [0,0,0,0,0,0,0,0,0,1/((x[0]-x[10])**2+(y[0]-y[10])**2)]
            ])
P = np.matrix(p)
Mm = np.matrix(MT@P@M).I#求MT*P*M的逆矩阵,@可以乘三个以上的矩阵
W=Mm@MT@P@Z#求最后的系数矩阵

3.结果提取与输出

    完整的W矩阵可以直接输出，但是为了方便地使用其中系数，将W由矩阵格式转化为array格式的w。

    按照完整流程，应将w中的各个系数带入公式Z=AX2+BXY+CY2+DX+EY+F中对Z值进行解算，但是由于此处已将待求点的X和Y值设置为0，所以F系数的值可以直接代表待求点的高程值。

    最后输出矩阵W和Z0即可，该待求点的高程约为17.575m。

w = np.array(W)#转为数列方便取值
z0=w[5,0]
# z0=w[0,0]*x[0]**2+w[1,0]*x[0]*y[0]+w[2,0]*y[0]**2+w[3,0]*x[0]+w[4,0]*y[0]+w[5,0] 完整版，但是x0，y0为0
# print('矩阵MT：\n',MT)
# print('矩阵Mm：\n',Mm) 以上为试验检查时写的
print('系数矩阵w：\n',W)
print("待定点的高程为：",z0)

    程序运行截图

ps：嘉人们谁懂啊，CSDN你的首行缩进在哪里？！