曲面拟合原理与实例
matlab曲线曲面拟合讲解及实例
例如：
x1=[1.0500 1.0520 1.0530 1.0900 1.0990 1.1020 1.1240 1.1420...
1.1490 1.0500 1.0520 1.0530 1.0900 1.0990 1.1020 1.1240 1.1420 1.1490];
x2=[3.8500 1.6500 2.7500 5.5000 7.7000 3.3000 4.9500 8.2500 11.5500...
52.5000 62.0000];
data=[x1;x2]; ห้องสมุดไป่ตู้ %类似于将x1 x2整合成一个2维数组。
a0= [-0.0014,0.07];
option=optimset('MaxFunEvals',5000);
format long;
同理用Symbolic Math Toolbox可以直接执行
>> ezmesh('X.^2 - Y.^2', [-2 2], [-2 2])
surf函数
在函数不能表示成z = f(x, y)时，需要用surf函数。比如x2+y2+z2=1.
先需要用柱面坐标或者球坐标来表示。这里用柱面坐标表示为 r2+z2=1
x = sqrt(1-z2)cosθ, x = sqrt(1-z2)sinθ;
执行matlab指令：
>> [theta, Z] = meshgrid((0:0.1:2)*pi, (-1:0.1:1));
>> X =sqrt(1 - Z.^2).*cos(theta);
>> Y =sqrt(1 - Z.^2).*sin(theta);
思拓力加权平均 平面拟合 曲面拟合 垂直平差
思拓力加权平均平面拟合曲面拟合垂直平差思拓力加权平均，平面拟合，曲面拟合，垂直平差是地球物理勘探中常用的数据处理技术。
本文将详细介绍这些技术的原理、应用和优缺点。
思拓力加权平均是处理GPS数据时经常使用的方法之一。
在实际情况中，由于各种误差的存在，每次测量所得的数据极其不稳定，这就需要采取一定的方法来消除这些误差。
思拓力加权平均是一种最小二乘法，是一种能够消除后验误差的参数估计法。
通过这种方法，可以得到更加精准的GPS观测结果，从而提高勘探的准确度和效率。
平面拟合是另一种常用的数据处理方法。
在地球物理勘探中，经常遇到一些复杂的地形和地貌，这就需要进行平面拟合，以方便对地形和地貌进行分析和研究。
平面拟合可以对大量的数据进行处理，并将这些数据最终汇集到一个平面上，以方便对这些数据进行分析和处理，从而得到准确的结果。
平面拟合可以用于处理地图、图像等各种数据，是一个非常重要的技术。
曲面拟合是一个加强版的平面拟合。
在地球物理勘探中，我们经常需要处理的是比较复杂的地形和地貌，因此平面拟合并不能满足我们的需求。
曲面拟合可以更好地适应各种复杂的地形和地貌，可以得到更加准确的结果。
曲面拟合是一种用于处理空间数据的重要方法，可以处理各种三维数据，例如地形图、地球模型等等。
垂直平差是一种用于高精度测量数据处理的方法。
在地球物理勘探中，垂直平差是处理GPS测量数据的一种重要技术。
它是一种非常精确的计算方法，通过将已知数据和未知数据结合起来，从而得到准确的结果。
垂直平差可以用于处理各种复杂的地形和地貌，可以精确到毫米级别的计算，是地球物理勘探中不可或缺的方法之一。
综上所述，思拓力加权平均，平面拟合，曲面拟合，垂直平差是地球物理勘探中常用的数据处理技术。
它们的优点在于可以准确地处理大量的数据，得到精确的结果。
缺点在于可能需要专业的软件和设备支持，有一定的学习门槛。
对于地球物理勘探工作者而言，熟练运用这些技术能够提高工作效率、提高勘探的准确性，为地球物理勘探事业的发展做出更大的贡献。
基于局部曲面拟合法
基于局部曲面拟合法
基于局部曲面拟合法是一种图形处理技术，它可以对图像中的特征的几何形状进行拟合，以获得最佳的几何特征表达。
局部曲面拟合法主要用于图像处理中的几何特性检测，其可以有效地构建出特征的精确结构表达。
局部曲面拟合法主要是通过拟合出合适的局部曲面，计算得到最佳的拟合曲面，以拟合图像中的特征。
具体来说，首先计算出曲面的三角剖分，然后按照局部曲面拟合的原理拟合出最佳曲面。
相比于其他拟合方法，局部曲面拟合法可以更好地描述几何特征，并可以可靠地定位几何特征。
局部曲面拟合法的精确度取决于输入的点集精度，因此，在多视图三维重构中，这种方法非常有用。
多视图三维重构可以重新构造出物体的三维几何形状，但是由于多个视图的投影误差，重构的三角分割精度不如期望值。
局部曲面拟合可以有效的改善这一问题，它可以起到可靠的三维几何表达，而且是一种高效的拟合方法。
此外，局部曲面拟合还可以应用于物体分割、对象检测和跟踪等应用领域。
物体分割是一种分离输入图像中不同物体的技术，它可以帮助更好地描述场景中的物体，从而有助于更准确地识别出输入图像中的物体。
局部曲面拟合可以有效地提取出更多的几何特征，从而更好地进行物体分割。
此外，在对象检测和跟踪中，局部曲面拟合可以有效提取出目标物体的特征，用于更准确地识别和跟踪物体。
总而言之，局部曲面拟合法可以有效构建出特征的精确结构表达，
这种方法可以用于多视图三维重构、物体分割、对象检测和跟踪等应用领域，在这些应用中都可以得到良好的结果。
拟合曲面函数
拟合曲面函数拟合曲面函数是数据分析中一项重要的技术，用于通过已知的数据点来建立一个连续的曲面函数。
在实际应用中，拟合曲面函数通常用于曲面拟合、数据可视化及预测等方面。
本文将介绍拟合曲面函数的相关知识。
1、什么是曲面拟合？曲面拟合是指用一个函数表示一组数据点所在的曲面。
它是拟合算法中的一种常见形式。
曲面拟合可以用于描绘地形、海洋气象、建筑设计等问题中。
其基本思想是，在误差最小化的约束下，尽可能地逼近数据点所在的曲面。
曲面拟合的理论基础是多项式拟合和最小二乘法。
2、曲面拟合的类型(1) 多项式拟合：多项式拟合是将数据点拟合到一条曲线或曲面上。
它的优点是较为简单，但是拟合的精度不如其他方法，不能适用于复杂的数据情况。
(2) 核函数方法：核函数方法是一种非参数方法，利用核函数对数据点进行进行拟合，适用于异常值较多的复杂数据情况。
(3) 计算机图形学方法：计算机图形学方法主要适用于曲面拟合和模型近似问题，它将数据表达为曲面网格，并对曲面进行分段处理，适用于采样有密度梯度的曲面。
3、曲面拟合方法(1) 插值法：在已知数据点间插值，得到一个连续的曲面。
插值法的优点是可以完全保证数据点被准确地拟合，但是对输入数据要求较高。
(2) 最小二乘法：使用最小二乘法模拟函数的拟合过程，得到一个拟合函数。
最小二乘法的优点是它对数据的要求不高，适用于大多数数据情况。
但是它不能完全保证数据点被准确地拟合，会产生一定的误差。
(3) 最大似然估计法：最大似然估计法是针对样本数据的拟合方法，通过优化统计模型参数，得到一个能够最优描述样本数据的模型。
最大似然估计法的优点是它可以根据数据情况选择不同的分布，适用于不同类型的数据。
二、基本步骤拟合曲面函数的基本步骤如下：1、读入数据：在拟合曲面函数之前，必须要读入一组有序的数据点。
2、选择拟合函数的类型：通过观察数据情况，根据实际需要选择适合的拟合函数类型。
4、确定拟合函数的参数：利用拟合方法，确定拟合函数的参数。
第五章 曲线拟合
)2
j 1
j 1 i0
对ak求偏导数（k＝0,1…m)
ak
nm
2
(
ai
x
i j
j1 i0
y
j
)
x
k j
0
m
m
n
化简得
ai
xik j
y
j
x
k j
i0 j 1
j 1
n
n
记
x
k j
Sk
y
j
x
k j
Tk
j 1
j 1 m
aiSki Tk (k 0,1m)
i0
写成矩阵形式
S0 S1 S2 Sm S1 S2 S3 Sm1 S2 S3 S4 Sm2 Sm Sm1 Sm2 S mm
i1 j 1
i 1
用矩阵形式给出即： AT Ax ATb 法方程组
例
用最小二乘法解下列超定方程组的近似解
2x1 x2 1 8x1 4x2 0 2x1 x2 1 7x1 x2 8 4x1 3
解： A=
2 1 8 4 2 1 7 1 4 0
2 1
AT A
如何衡量接近程度？
最小二乘原理
一、什么是最小二乘原理
是衡量接近程度的一种方法
x x0 x1 xn 已知 y y0 y1 yn
设p(x) a0 a1x an xn
n
n
求a0 , a1,an 使 Ri2 (P(xi ) yi )2 最小。
io
i0
用最小二乘原理进行曲线拟合的方法称为最小二乘法。
这里（m<n)，适当的选取 a0 , a`,am 使得
n
(a0 , a1,am ) [ p(x j ) y j ]2 为最小值 j 1
B样条曲面拟合
谢谢
这样的 节点矢量定义了一般非均匀B样条基.
分段贝齐尔曲线 其节点矢量中二端节点重复度与类型2相同, 为k+1。所不同的是,点数减1必 须等于次数的正整数倍,
二、曲面拟合的条件
待拟合曲面的控制点数 样条基次数 两个方向的节点向量
B样条曲面拟合
一、B样条曲线的种类及主要性质
均匀B样条曲线 均匀B样条曲线节点矢量中节点为沿参数轴 均匀或等距分布,所有节点区间长=Δi=ti+ti常数 >0(i=0,1,…,n+k)。这样的节点矢量定义了均匀 B样条基。 准均匀B样条曲线 其节点矢量中二端节点具有重复度k+1,即t0 =t1=…=tk,tn+1=tn+2=…=tn+k+1,所有内 节点均匀分布,重复度为1。
3、节点向量
节点向量是由数据点的参数化确定的，给定 曲面的数据点赋两个参数值，使位于拟和曲 面上的这些点与平面参数域内的点建立一一 对应的关系。 参数化的手段有：均与参数化、累加弦长参 数化、基面参数化。
得到数据点集的参数值后，就可以利用最小 二乘拟合反算曲面控制点
三、曲面拟合的实例
例1： 图1是一条三次B样条曲线与控制多 边形，符号“○”表示曲线初始控制点， “*”号表示曲线经一次4尺度加细后得到的 控制点，由图易见，加细后的控制点与曲线 更加接近。因此，可通过控制点的多尺度迭 代来生成B样条曲线增加了生成样条曲线的 速度。
例2 图2是利用第2节算法进行4尺度三层加细得 到的双三次B样条拟合曲面，0D为 (1+3)×(1+3)网格，经过三次4尺度加细 后3D为(64+3)×(64+3)网格，重建的B 样条曲面即有很好的光滑形状，又有很高的 逼近精度，同时，在上机实现时算法速度比 以前基于Binary细分多层加细算法的速度稍 快，主要由于减少了耗时的Step 2的次数。
曲面方程拟合算法 python
曲面拟合是一个在各个领域都有广泛应用的主题。
在Python中，你可以使用numpy、scipy 等库来实现。
以下是一个简单的二次曲面拟合的例子：首先，我们需要导入所需的库：pythonimport numpy as npfrom scipy.optimize import curve_fitimport matplotlib.pyplot as plt假设我们有一些数据点(x, y)，我们希望找到一个表达式f(x, y) = ax^2 + by^2 + cxy + dx + ey + f 来拟合这些数据。
其中a, b, c, d, e 和f 是我们需要找到的参数。
我们可以定义一个函数来拟合我们的数据：pythondef func(x, y, a, b, c, d, e, f):return a*x**2 + b*y**2 + c*x*y + d*x + e*y + f然后，我们可以使用curve_fit 函数来找到最佳参数：pythonpopt, pcov = curve_fit(func, (x, y), z)这里，(x, y) 是我们的输入数据点，z 是我们的目标数据。
curve_fit 会返回两个值：popt 是最优参数，pcov 是参数的协方差矩阵。
最后，我们可以使用matplotlib 来可视化我们的拟合结果：pythonplt.scatter(x, y, color='b') # 散点图表示原始数据点plt.plot(x, func(x, y, *popt), color='r') # 红色曲线表示拟合的曲面plt.show()注意，这只是一个非常简单的例子。
在实际应用中，你可能需要考虑更多的因素，例如噪声、异常值等。
此外，如果数据的维度很高，或者你希望拟合的模型很复杂，你可能需要使用更复杂的优化算法，例如Levenberg-Marquardt 算法等。
曲面拟合是啥原理图的应用
曲面拟合是啥原理图的应用1. 曲面拟合的概念曲面拟合是一种数学建模技术，用于将一组离散点数据拟合成平滑的曲面。
它通过寻找最适合给定点集的曲面来实现数据的近似和拟合。
曲面拟合在计算机图形学、CAD/CAM、工程设计和地理信息系统等领域得到了广泛应用。
2. 曲面拟合的原理曲面拟合的原理基于数学最优化方法，旨在找到一个曲面模型，使其最接近给定的离散点数据。
常见的曲面拟合方法包括最小二乘法和样条曲面拟合等。
2.1 最小二乘法最小二乘法是曲面拟合中常用的一种方法。
它通过最小化数据点与曲面之间的距离来确定最佳拟合曲面。
最小二乘法可以分为线性最小二乘法和非线性最小二乘法。
2.1.1 线性最小二乘法线性最小二乘法适用于拟合线性模型的情况。
其基本原理是建立一个与数据点相匹配的线性模型，并通过最小化残差平方和来确定最佳拟合曲面。
线性最小二乘法的数学公式可以表示为：min E = Σ (yi - f(xi))^2其中，E为残差平方和，yi为实际观测值，f(xi)为线性模型的预测值。
2.1.2 非线性最小二乘法非线性最小二乘法适用于拟合非线性模型的情况。
其原理与线性最小二乘法类似，不过在计算残差平方和时，需要通过迭代的方式逼近最佳拟合结果。
非线性最小二乘法的数学公式可以表示为：min E = Σ (yi - f(xi;θ))^2其中，θ为模型参数，f(xi;θ)为非线性模型的预测值。
2.2 样条曲面拟合样条曲面拟合是一种使用控制点和插值方法构造曲面的技术。
它将拟合问题转化为一个插值问题，在给定的控制点上生成一个平滑的曲面。
样条曲面拟合的原理是通过插值方法将数据点与控制点相连，并在控制点上生成一个曲面模型，以实现数据的拟合。
3. 曲面拟合的应用曲面拟合在许多领域都有广泛的应用。
以下是一些常见的应用场景：•计算机图形学：曲面拟合可以用于生成光滑的曲线和曲面，用于渲染和动画效果的生成。
•CAD/CAM：曲面拟合可以用于设计和制造曲面形状的产品，例如汽车、飞机等。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
问题：给定一组坐标(,,)g g g x y z ，1,2,g =…,n ，表示有n 个点。
要求用以下二元多项式函数对所给的坐标进行拟合：,11111,111(,)p qp qi j i j ijij ij i j f x y a xya x y ----=====∑∑∑即211112131212122232111211123111211123(,)q q q q i i i i q i i i iq p p p p q p p p pq a a y a y a y f x y a x a xy a xy a xy a x a x ya x y a x y a x a x y a x y a x y ------------++++=+++++++++++++++L L M L M L设1112121222221211,,q q p p pq p q a a a x y a a a x y a a a x y ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦⎣⎦x y A L LM M O M M M L则函数又可表示为(,)T f =x y x Ay ，拟合的目标就是求出系数矩阵A 。
最小二乘法：构造关于系数ij a 的多元函数：2112111111(,,)[(,)]()pqnni j pq g g g g g ij g g g i j s a a f x y z a x y z ωω--=====-=-∑∑∑∑L点（11a ，…，pq a ）是多元函数11(,,)pq s a a L 的极小点，其中g ω为权函数，默认为1，所以点（11a ，…，pq a ）必须满足方程组0ijsa ∂=∂ 在1g ω=的情况下，有{}21111111111[(,)]2[(,)][(,)]2[(,)]2(,)nggg g ij ijng g g g g g ij ni j g g g g g g n i j i j g gg g g g g g s f x yz a a f x y z f x y a f x y z x y x y f x y x y z ==--=----=∂∂=-∂∂⎧⎫∂⎪⎪=-⎨⎬∂⎪⎪⎩⎭=-⎡⎤=-⎣⎦∑∑∑∑因此可得111111(,)nni j i j ggg g g g g g g xyf x y x y z ----===∑∑1111111111p qnni j i j ggggg g g g g xya x y x y z αβαβαβ------=====∑∑∑∑ ,11111111,11p qnni j i j ggg gg g g g g xya x y x y z αβαβαβ------====∑∑∑ ,1111111,111()p qn n i j i j g g g g g g g g g a x y x y x y z αβαβαβ------===⎡⎤=⎢⎥⎣⎦∑∑∑ 令11111(,)()ni j g g g gg u i j x y x yαβαβ----==∑，111(,)ni j gg g g v i j x y z --==∑ 则,1,1(,)(,)p qa u i j v i j αβαβαβ==∑(,)(11),(,)i j p q =，…, 上式实际共有p q ⨯个等式，可将这p q ⨯个等式写成矩阵的形式有：111111(1,1)(1,1)(1,1)(,)(,)(,)pq pq pq u u a v u p q u p q a v p q ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦L M O M M M L 也就是U*a=V 的形式，其中1111(1,1)(1,1)(,)(,)pq pq u u u p q u p q ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦U L M OM L ，11pq a a ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦a M ，(1,1)(,)v v p q ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦V MU 为pq pq ⨯阶矩阵，实现函数为function A=leftmatrix(x,p,y,q)；V 为长pq 的列向量，实现函数为function B=rightmatrix(x,p,y,q,z)。
这样就可以算出列矩阵a ,然后转化成A 。
例子：某地区有一煤矿，为估计其储量以便于开采，先在该地区进行勘探。
假设该地区是一长方形区域，长为4公里，宽为5公里。
经勘探得到如下数据：请你估计出此地区内（51,42≤≤≤≤y x ）煤的储量，单位用立方米表示，并用电脑画出该煤矿的三维图象。
如果直接画出三维曲面图形：clear;x=1:4;y=1:5;[X,Y]=meshgrid(x,y) Z=[13.72 25.80 8.47 25.27 22.32; 15.47 21.33 14.49 24.83 26.19; 23.28 26.48 29.14 12.04 14.58; 19.95 23.73 15.35 18.01 16.29]' surf(X,Y,Z);X =1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 1 2 3 4 Y =1 1 1 12 2 2 23 3 3 34 4 4 45 5 5 5 Z =13.7200 15.4700 23.2800 19.950025.8000 21.3300 26.4800 23.73008.4700 14.4900 29.1400 15.350025.2700 24.8300 12.0400 18.010022.3200 26.1900 14.5800 16.2900粗略计算体积：底面积乘以平均高度。
p=sum(Z);q=p(:,[2,3,4]);h=sum(q')/15v=2000*4000*hh =20.0773v =1.6062e+008进行线性插值：xi=linspace(1,4,31);yi=linspace(1,5,41);[XI,YI]=meshgrid(xi,yi); ZI=interp2(X,Y,Z,XI,YI,'linear');surf(XI,YI,ZI);进行三次多项式插值：xi=linspace(1,4,31);yi=linspace(1,5,41);[XI,YI]=meshgrid(xi,yi); ZI=interp2(X,Y,Z,XI,YI,'cubic');surf(XI,YI,ZI);进行插值后计算体积：底面积乘以平均高度。
xi=linspace(1,4,61);yi=linspace(1,5,81);[XI,YI]=meshgrid(xi,yi); ZI=interp2(X,Y,Z,XI,YI,'cubic');surf(XI,YI,ZI);H=0;n=0;for j=21:61for i=1:81H=H+ZI(i,j);n=n+1;endendnH=H/nS=2000*4000;V=S*Hn =3321H =20.8222V =1.6658e+008上面是插值的方法解题，下面用拟合的方法解题。
为此编写了几个M函数：leftmatrix.mfunction U=leftmatrix(x,p,y,q)% U*a=V a为系数列矩阵，长度为p*q% U为左边p*q乘p*q矩阵% x,y 为长度一致的列矩阵，给定点的坐标% p,q 为拟合的函数中x,y的幂的最高次数m=length(x);if (nargin~=4) & (m~=length(y))error('error check check!');endU_length=p*q; % U 为p*q阶方阵U=zeros(U_length,U_length); % 赋值0，目的是分配内存for i=1 : p*qfor j= 1 : p*qx_z=quotient(j-1,q)+quotient(i-1,q); % x 的幂的次数，quotient为求商 y_z=mod(j-1,q)+mod(i-1,q); % y 的幂的次数U(i,j)=qiuhe(x,x_z,y,y_z);endendrightmatrix.mfunction V=rightmatrix(x,p,y,q,z)% U*a=V% V 为一个列向量长为p*q% x y z 为点的坐标%p q 分别为x y幂的最高次数if nargin~=5error('error check check! rightmatrix')endV=zeros(p*q,1);for i=1 : p*qx_z=quotient(i-1,q);y_z=mod(i-1,q);V(i,1)=qiuhe(x,x_z,y,y_z,z);endquuotient.mfunction sh=quotient(x,y)% sh 为 x/y 的商sh=(x-mod(x,y))/y;qiuhe.mfunction he=qiuhe(x,p,y,q,z)% he x^p*y^q 从1－>m的和% x,y 向量长度相同% p,q分别为x,y的幂的次数m=length(x);if (nargin<4 )&(m~=length(y)) %输入量至少为四，x,y行向量长度必需一样error('error check check!');endif nargin==4 %没有 z , 默认为元素全部为1的向量z=ones(m,1);endhe=0;for i=1:mhe=he+x(i)^p * y(i)^q*z(i); % 1-->m 求和end下面一段程序先进行拟合，然后验证拟合的效果，具体操作：先输入x=…y=…z=…p=…q= (注意x,y,z是向量)；拟合得到系数a，也就是得到了拟合的函数；根据拟合函数计算给定点(xx, yy)的函数值zz=f(xx, yy)并进行画图检验。
程序保存于M文件fit.m。
fit.mclear;[X,Y]=meshgrid(1:4,1:5);Z=[13.72 25.80 8.47 25.27 22.32;15.47 21.33 14.49 24.83 26.19;23.28 26.48 29.14 12.04 14.58;19.95 23.73 15.35 18.01 16.29]';x=reshape(X,20,1);y=reshape(Y,20,1);z=reshape(Z,20,1);p=4;q=5;U=leftmatrix(x,p,y,q); % U*a_n=VV=rightmatrix(x,p,y,q,z);%a_n=inv(U)*V;a_n=U\V;for i=1 : length(a_n) % 把长为p*q 的列向量a_n转换成p*q的矩阵aa ii=quotient(i-1,q)+1; % quotient求商jj=mod(i-1,q)+1;aa(ii,jj)=a_n(i,1);endaam=31;n=41;%m=4;n=5;[XI,YI]=meshgrid(linspace(1,4,m),linspace(1,5,n));xx=reshape(XI,m*n,1);yy=reshape(YI,m*n,1);zz=zeros(m*n,1);xy=zeros(m*n,1);xt=zeros(m*n,1);yt=zeros(m*n,1);%zz=0; % zz 是 xx,yy 代入所拟合的函数求出的函数值for i=1 : p % 函数为Σaa(i,j)*x^i*y^j，(i=1...p,j=1...q) for j=1 : q % aa 为pxq的系数的矩阵xt=xx.^(i-1);yt=yy.^(j-1);xy=xt.*yt;zz=zz+aa(i,j).*xy;endendZI=reshape(zz,n,m);surf(XI,YI,ZI); %axis([1 4 1 5 0 30])aa =1.0e+003 *0.1465 -0.2678 0.2132 -0.0624 0.0058-0.7287 1.3972 -0.9275 0.2412 -0.02100.4416 -0.8415 0.5487 -0.1407 0.0122-0.0680 0.1295 -0.0839 0.0214 -0.0018注意：权函数在拟合的函数非常重要，不过她只能按你遇到的具体问题来取，我这里为1；当p，q越大时，拟合的函数与原数据的方差越小，但是有可能函数本身抖动非常厉害，可以画图看出来。