精华内容
下载资源
问答
  • 数据插值可以根据有限个点的取值状况,合理估算出附近其他点的取值,从而节约大量的实验和测试资源,节省大量的人力、物力和财力。 ②数据插值能够根据已知数据推算未知数据,这使得人们解决问题的能力得到了拓展...

    在这里插入图片描述

    ①数据插值可以根据有限个点的取值状况,合理估算出附近其他点的取值,从而节约大量的实验和测试资源,节省大量的人力、物力和财力。
    ②数据插值能够根据已知数据推算未知数据,这使得人们解决问题的能力得到了拓展和延伸。

    1、引例-零件加工问题

    例1、在飞机制造中,机翼的加工是一项关键技术。由于机翼尺寸很大,通常在图纸中只能标出一些关键点的数据。下表给出了某型飞机机翼的下缘轮廓线数据,求x每改变0.1时y的值。
    在这里插入图片描述在这里插入图片描述
    它前段采样点稀疏,后段采样点密集,说明这段曲线前面可能比较规律、平滑,后段比较复杂

    x = [0 3 5 7 9 11 12 13 14 15];
    y = [0 1.2 1.7 2 2.1 2 1.8 1.2 1 1.6];
    plot(x,y)    %为插值前的曲线
    hold on;
    x1 = 0:0.1:15;    %x每改变0.1要插一个值
    y1 = interp1(x,y,x1,'spline');    %插值函数interp1:计算出这些插值点在y方向上的值存入y1中,y1很长
    plot(x1,y1)   %完成插值后的曲线,曲线变得光滑很多
    legend('y原曲线','y1插值后的曲线')
    title('零件加工问题','color','r','fontsize',18)
    

    在这里插入图片描述
    2、数据插值的计算机制
    从数学上来说,数据插值是一种函数逼近的方法。
    在这里插入图片描述(1)interp1( ):一维插值函数。调用格式:
    Y=interp1(X,Y,X1,method)
    根据X、Y的值,计算函数在×1处的值。其中,X、Y是两个等长的已知向量,分别表示采样点和采样值。X是一个向量或标量,表示要插值的点。

    3、数据插值的实现方法
    (1)一维插值method用于指定插值方法,常用的取值有以下四种:
    ①linear:线性插值,默认方法。
    将与插值点靠近的两个数据点用直线连接,然后在直线上选取对应插值点的数据。
    ②nearest:最近点插值。
    选择最近样本点的值作为插值数据。如果是中间点,则取后一个数据点的值。
    ③pchip:分段3次埃尔米特插值。
    采用分段三次多项式,除满足插值条件,还需满足在若干节点处相邻段插值函数的一阶导数相等,使得曲线光滑的同时,还具有保形性。
    ④ spline :3次样条插值。
    每个分段内构造一个三次多项式,使其插值函数除满足插值条件外,还要求在各节点处具有连续的一阶和二阶导数。(进—步提高了曲线的光滑性)

    例2:以例1为例,进行四种数据插值的实现

    x = [0 3 5 7 9 11 12 13 14 15];
    y = [0 1.2 1.7 2 2.1 2 1.8 1.2 1 1.6];
    x1 = 0:0.1:15;    %x每改变0.1要插一个值
    % plot(x,y)
    % hold on
    subplot(2,2,1)
    y1 = interp1(x,y,x1,'linear');    %插值函数interp1:计算出这些插值点在y方向上的值存入y1中,y1很长
    plot(x,y,x1,y1)   %完成插值后的曲线
    xlabel('interp1(x,y,x1,''linear'')','color','b')
    legend('y','y1','location','northwest')
    title('method = ''linear''','color','r','fontsize',14)
    
    subplot(2,2,2)
    y1 = interp1(x,y,x1,'nearest');    %插值函数interp1:计算出这些插值点在y方向上的值存入y1中,y1很长
    plot(x,y,x1,y1)   %完成插值后的曲线
    xlabel('interp1(x,y,x1,''nearest'')','color','b')
    legend('y','y1','location','northwest')
    title('method = ''nearest''','color','r','fontsize',14)
    
    subplot(2,2,3)
    y1 = interp1(x,y,x1,'pchip');    %插值函数interp1:计算出这些插值点在y方向上的值存入y1中,y1很长
    plot(x,y,x1,y1)   %完成插值后的曲线
    xlabel('interp1(x,y,x1,''pchip'')','color','b')
    legend('y','y1','location','northwest')
    title('method = ''pchip''','color','r','fontsize',14)
    
    subplot(2,2,4)
    y1 = interp1(x,y,x1,'spline');    %插值函数interp1:计算出这些插值点在y方向上的值存入y1中,y1很长
    plot(x,y,x1,y1)   %完成插值后的曲线
    xlabel('interp1(x,y,x1,''spline'')','color','b')
    legend('y','y1','location','northwest')
    title('method = ''spline''','color','r','fontsize',14)
    

    在这里插入图片描述
    四种方法的比较:
    (1)线性插值和最近点插值方法比较简单。其中线性插值方法的计算量与样本点n无关。n越大,误差越小。
    (2)3次埃尔米特插值和3次样条插值都能保证曲线的光滑性。相比较而言,3次埃尔米特插值具有保形性;而3次样条插值要求其二阶导数也连续,所以插值函数的性态更好。

    (2)interp2():二维插值函数。调用格式:
    Z1=interp2(X,Y,Z,X1,Y1,method)

    其中,X、Y是两个向量,表示两个参数的采样点,Z是采样点对应的函数值。X1、YI是两个标量或向量,表示要插值的点。method与一维插值方法相同,但不支持pchip方法。

    为什么这两种插值方法都用3次多项式而不用更高次的?
    答:多项式次数并非越高越好。次数越高,越容易产生震荡而偏离原函数,这种现象称为龙格(Runge)现象。

    在这里插入图片描述

    展开全文
  • matlab代码 [x,y]=meshgrid(-3:1:3); z=peaks(x,y) [xi,yi]=meshgrid(-3:0.25:3); figure(1) surfc(x,y,z) title('原始数据') zi1=interp2(x,y,z,xi,yi,'spline') zi2=interp2(x,y,z,xi,yi,'linear') zi3=interp2(x,...

    matlab代码

    [x,y]=meshgrid(-3:1:3);
    z=peaks(x,y)
    [xi,yi]=meshgrid(-3:0.25:3);
    figure(1)
    surfc(x,y,z)
    title('原始数据')
    zi1=interp2(x,y,z,xi,yi,'spline')
    zi2=interp2(x,y,z,xi,yi,'linear')
    zi3=interp2(x,y,z,xi,yi,'nearest')
    zi4=interp2(x,y,z,xi,yi,'cubic')
    figure(2)
     subplot(2,2,1);
     surf(xi,yi,zi1);
     subplot(2,2,2);
     surf(xi,yi,zi2);
     subplot(2,2,3);
     surf(xi,yi,zi3);
     subplot(2,2,4);
     surf(xi,yi,zi4);
    

    在这里插入图片描述
    有棱有角
    在这里插入图片描述
    分别为:spline linear nearest cubic
    cubic:在这里插入图片描述

    在这里插入图片描述

    图像处理

    比较不同的差值方法(图像处理方面)
    最邻近元法The nearest计算量较小,但可能会造成插值生成的图像灰度上的不连续,在灰度变化的地方可能出现明显的锯齿状。
    双线性内插法Bilinear Interpolation的计算比最邻近点法复杂,计算量较大,但没有灰度不连续的缺点,结果基本令人满意。它具有低通滤波性质,使高频分量受损,图像轮廓可能会有一点模糊。
    三次内插法Cubic interpolation:该方法利用三次多项式S(x)求逼近理论上最佳插值函数sin(x)/x, 其数学表达式为:
    在这里插入图片描述
    待求像素(x, y)的灰度值由其周围16个灰度值加权内插得到
    “Inverse Distance to a Power(反距离加权插值法)”、
    “Kriging(克里金插值法)”、
    “Minimum Curvature(最小曲率)”、
    “Modified Shepard’s Method(改进谢别德法)”、
    “Natural Neighbor(自然邻点插值法)”、
    “Nearest Neighbor(最近邻点插值法)”、
    “Polynomial Regression(多元回归法)”、
    “Radial Basis Function(径向基函数法)”、
    “Triangulation with Linear Interpolation(线性插值三角网法)”、
    “Moving Average(移动平均法)”、
    “Local Polynomial(局部多项式法)”

    空间分析

    空间分析中的插值
    在这里插入图片描述

    优点比较
    1)OK 插值法:以空间统计学作为理论基础,可以克服内插中误差难以分析的问题,能够对误差作为这逐点的理论估计,不会产生回归分析的边界效应,插值精度较高,唯一性很强,外推能力较强。
    2)TS 插值法:不需要对空间结构进行预先估计和作统计假设,只进行局部区块的拟合,用以补充或修改局部区块的空间变量分布曲面,而不用处理不涉及局部区块修改的其他部分的数据,当表面很平滑时,也不牺牲精度。
    3)IDW 插值法:计算开销少,具有普适性,不需要根据数据的特点对方法加以调整,当样本数据的密度足够大时,几何方法一般能达到满意的精度。
    缺点比较
    1)OK 插值法:作为一种统计学方法,复杂,计算量大,运算速度慢,变异函数需要根据经验人为选定。
    2)TS 插值法:作为一种函数方法,难以满足对于利用有限的观测数据进行缺值预测和内插格网的精度要求,也难以对误差进行估计,样本点稀疏时插值效果不好。
    3)IDW 插值法作为一种几何方法,插值结果受 r 值的影响很大,根据不同 r 值估算的同一未知点的值会有很大的差别。当任何一个 d r i = 0 ( 1 ≤ i ≤ n , 1 ≤ r ≤ n ) dri=0(1≤i≤n,1≤r≤n) dri=01in1rn时,该点权值为无穷大,导致该点的输出数据不连续,计算时会得到其实际测量值,在进行外插值时,IDW会不恰当地将这些估计值回归为观测数据的平均值;当选用的插值距离的次数为偶数时,dir 为非负数,插值结果Z总是满足: Z m i n < Z < Z m a x [ 5 ] Zmin<Z<Zmax[5] Zmin<Z<Zmax[5]。因此,在山峰、山谷及数据点以外的区域,IDW 的插值常常产生一些与我们直觉相矛盾的结果。IDW 权值的大小直接影响了标准差的大小,决定了插值的整体精度情况,对插值结果造成重大的影响。

    3 种方法的适用范围比较
    1)OK 插值法适用性较广,普遍用于局部的、区域较小的范围,需满足内蕴假设,其区域化变量的平均值是未知的常数[6]。
    2)TS 插值法适合于非常平滑的表面,一般要求有连续的一阶和二阶导数,适用于密度较大的点内插求等值线,对整体的、点的数量较多的数据有很好的插值效果。
    3)IDW 插值法适用于整体的样本点的密度较大且样本点的分布比较均匀的数据。

    matlab代码的解释:

    meshgrid
    [X,Y] = meshgrid(x,y) 基于向量 x 和 y 中包含的坐标返回二维网格坐标。X 是一个矩阵,每一行是 x 的一个副本;Y 也是一个矩阵,每一列是 y 的一个副本。坐标 X 和 Y 表示的网格有 length(y) 个行和 length(x) 个列。
    示例
    [X,Y] = meshgrid(x) 与 [X,Y] = meshgrid(x,x) 相同,并返回网格大小为 length(x)×length(x) 的方形网格坐标。
    示例
    [X,Y,Z] = meshgrid(x,y,z) 返回由向量 x、y 和 z 定义的三维网格坐标。X、Y 和 Z 表示的网格的大小为 length(y)×length(x)×length(z)。
    示例
    [X,Y,Z] = meshgrid(x) 与 [X,Y,Z] = meshgrid(x,x,x) 相同,并返回网格大小为 length(x)×length(x)×length(x) 的三维网格坐标。
    在这里插入图片描述

    surfc

    surfc(X,Y,Z) 使用 Z 来代表颜色数据和曲面高度。X 和 Y 是用于定义曲面的 x 和 y 分量的向量或矩阵。如果 X 和 Y 为向量,length(X) = n 且 length(Y) = m,其中 [m,n] = size(Z)。在这种情况下,曲面顶点是 (X(j), Y(i), Z(i,j)) 三元组。要创建任意域的 X 和 Y 矩阵,请使用 meshgrid 函数。

    [X,Y,Z] = peaks(30);
    figure
    surfc(X,Y,Z)
    

    在这里插入图片描述

    peaks
    peaks 是从高斯分布转换和缩放得来的包含两个变量的函数,在演示 mesh、surf、pcolor、contour 等函数中很有用。
    Z = peaks; 返回一个 49×49 矩阵。
    Z = peaks(n); 返回一个 n×n 矩阵。
    Z = peaks(V); 返回一个 n×n 矩阵,其中 n = length(V)。
    Z = peaks(X,Y); 在给定的 X 和 Y(必须大小相同)处计算 peaks 并返回大小相同的矩阵。
    peaks(…)(无输出参数)使用 surf 绘制 peaks 函数。使用先前语法中的任意输入参数组合。
    [X,Y,Z] = peaks(…); 返回另外两个矩阵 X 和 Y 用于参数绘图,例如 surf(X,Y,Z,del2(Z))。如未作为输入参数给出,基础矩阵 X 和 Y 是
    [X,Y] = meshgrid(V,V)
    其中 V 是给定向量,或者 V 是长度为 n 的向量,其元素从 -3 到 3 均匀间隔。如果未给出输入参数,默认的 n 是 49。
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    n增大,peak(n)提高来采样率
    interp2
    Interpolation — increase sample rate by integer factor
    提高采样率的整数倍
    Syntax
    y = interp(x,r)
    y = interp(x,r,n,alpha)
    [y,b] = interp(x,r,n,alpha)
    Description
    y = interp(x,r) increases the sample rate of x, the input signal, by a factor of r.
    y = interp(x,r,n,alpha) specifies two additional values:
    n is half the number of original sample values used to interpolate the expanded signal. n是用于插值扩展信号的原始样本值的一半。
    alpha is the normalized cutoff frequency of the input signal, specified as a fraction of the Nyquist frequency. alpha是输入信号的归一化截止频率,指定为奈奎斯特频率的一部分。
    [y,b] = interp(x,r,n,alpha) also returns a vector, b, with the filter coefficients used for the interpolation.返回一个向量b,带有用于插值的滤波系数。

    z是怎么来的:

    z=peaks(x,y)
    

    7*7的矩阵

    z =
    
        0.0001    0.0034   -0.0299   -0.2450   -0.1100   -0.0043   -0.0000
        0.0007    0.0468   -0.5921   -4.7596   -2.1024   -0.0616    0.0004
       -0.0088   -0.1301    1.8559   -0.7239   -0.2729    0.4996    0.0130
       -0.0365   -1.3327   -1.6523    0.9810    2.9369    1.4122    0.0331
       -0.0137   -0.4808    0.2289    3.6886    2.4338    0.5805    0.0125
        0.0000    0.0797    2.0967    5.8591    2.2099    0.1328    0.0013
        0.0000    0.0053    0.1099    0.2999    0.1107    0.0057    0.0000
    

    引用自 刘光孟 空间分析中几种插值方法的比较

    展开全文
  • 数据预处理——插值算法matlab实现

    千次阅读 2020-11-24 14:23:44
    常常需要根据已知的函数点进行数据、模型的处理和分析,而有时候现有的数据是极少的,不足以支撑分析的进行,这时就需要使用一些数学的方法,“模拟产生”一些新的但又比较靠谱的值来满足需求,这就是插值的作用。...

    笔记来自《清风数学建模》系列课程课程链接

    目录

    一、插值算法的用途

    二、一维插值问题

    三、分段线性插值

    1.简单分段线性插值

     2. 拉格朗日插值与牛顿插值

    3. 埃尔米特 (Hermite)插值

    (1)原理

    (2)分段三次埃尔米特插值编程实现(建模比赛中最常用)

    4. 三次样条插值

    (1)原理

    (2)三次样条插值编程实现(建模中最常用)

    5. n维数据的插值(了解)

     四、插值算法用于短期预测

     五、实例——数据预处理(补全数据)

    1. 三次埃米尔特插值

    2. 三次样条插值

    3. 三次样条插值和三次埃尔米特插值的对比



    一、插值算法的用途

    数模比赛中,常常需要根据已知的函数点进行数据、模型的处理和分析,而有时候现有的数据是极少的,不足以支撑分析的进行,这时就需要使用一些数学的方法,“模拟产生”一些新的但又比较靠谱的值来满足需求,这就是插值的作用。

    二、一维插值问题

     

     

    限制多项式的次数,插值的函数是唯一的。 但是插值多项式次数越高误差越小吗?

    红色的为标准正确的函数图像。后面是不同插值次数的图像。可以发现,高次插值会产生龙格现象,即在两端处波动极大,产生明显的震荡。在不熟悉曲线运动趋势的前提下,不要轻易使用高次插值。

    由于插值多项式的次数高但是精度未必高,相反,次数越高可能产生的误差越大。为了减少这种误差,我们可以使用分段线性插值。

    三、分段线性插值

    1.简单分段线性插值

    如何提高插值的精度呢?我们可以采用分段线性插值。比如要得到C点的y值,可以得到A、B两点构成的线段的函数关系式,通过函数关系式求得C点的y值。线性是值找离它最近的两个点求得线性函数进行求解。

    也可以选择离C点最近的3个点计算得到抛物线的函数关系式从而进行求解。这种就叫分段二次插值(实际建模时用的比较多)。

     2. 拉格朗日插值与牛顿插值

     上述两种插值仅仅要求插值多项式在插值节点处与被插函数有相等的函数值,而这种插值多项式却不能全面反映被插值函数的性态。然而在许多实际问题中,不仅要求插值函数与被插值函数在所有节点处有相同的函数值,它也需要在一个或全部节点上插值多项式与被插函数有相同的低阶甚至高阶的导数值。对于这些情况,拉格朗日插值和牛顿插值都不能满足。

    3. 埃尔米特 (Hermite)插值

    (1)原理

    不但要求在节点上的函数值相等,而且还要求对应的导数值也相等,甚至要求高阶导数也相等,满足这种要求的插值多项式就是埃尔米特插值多项式。(保证函数和导数相等)

     

    (2)分段三次埃尔米特插值编程实现(建模比赛中最常用)

    直接使用Hermite插值得到的多项式次数较高,也存在着龙格现象,因此在实际应用中,往往使用分段三次 Hermite 插值多项式 (PCHIP)。

    Matlab有内置的函数(实现过程已经帮我们封装好了,会调用就行了):

    p = pchip(x,y, new_x)

    x是已知的样本点的横坐标;y是已知的样本点的纵坐标;new_x是要插入处对应的横坐标

    % 分段三次埃尔米特插值
    x = -pi:pi;%题目中给定的x值步长为1
    y = sin(x); 
    new_x = -pi:0.1:pi;%需要插值的x值步长为0.1
    p = pchip(x,y,new_x);
    plot(x, y, 'o', new_x, p, 'r-')

    plot函数用法:
    plot(x1,y1,x2,y2)
    线方式: ‐ 实线 :点线 ‐. 虚点线 ‐ ‐ 波折线
    点方式: . 圆点 +加号 * 星号 x x形 o 小圆
    颜色: y黄; r红; g绿; b蓝; w白; k黑; m紫; c青

    4. 三次样条插值

    (1)原理

    (2)三次样条插值编程实现(建模中最常用)

    Matlab有内置的函数:

    p = spline (x,y, new_x)

    x是已知的样本点的横坐标;y是已知的样本点的纵坐标;new_x是要插入处对应的横坐标

    % 三次样条插值和分段三次埃尔米特插值的对比
    x = -pi:pi; 
    y = sin(x); 
    new_x = -pi:0.1:pi;
    p1 = pchip(x,y,new_x);   %分段三次埃尔米特插值
    p2 = spline(x,y,new_x);  %三次样条插值
    figure(2);
    plot(x,y,'o',new_x,p1,'r-',new_x,p2,'b-')
    legend('样本点','三次埃尔米特插值','三次样条插值','Location','SouthEast')   %标注显示在东南方向

    【分析】通过图发现,在这个例子中,三次样条插值的图像更为光滑,精确度要更高一些,更加接近于sin(x)这个原始函数。

    可以看出,三次样条生成的曲线更加光滑。在实际建模中,由于我们不知道数据的生成过程,因此这两种插值都可以使用。

    说明:

    legend(string1,string2,string3, …)

    % 分别将字符串1、字符串2、字符串3……标注到图中,每个字符串对应的图标为画图时的图标。
    % ‘Location’用来指定标注显示的位置

    figure(1),figure(2)是为了防止多个图像存在时图像被覆盖。

    5. n维数据的插值(了解)

    p = interpn(x1,x2,...,xn, y, new_x1,newx_2,...,new_xn, method)

    x1,x2,...,xn是已知的样本点的横坐标
    y是已知的样本点的纵坐标坐标
    new_x1,newx_2,...,new_xn是要插入点的横坐标
    method是要插值的方法
    ‘linear’:线性插值(默认算法);
    ‘cubic’:三次插值;
    ‘ spline’ :三次样条插值法; ( 最为精准 )
    ‘nearest’:最邻近插值算法。

    % n维数据的插值
    x = -pi:pi; y = sin(x); 
    new_x = -pi:0.1:pi;
    p = interpn (x, y, new_x, 'spline');
    % 等价于 p = spline(x, y, new_x);
    figure(3);
    plot(x, y, 'o', new_x, p, 'r-')

     

     四、插值算法用于短期预测

    根据过去10年的中国人口数据,预测接下来三年的人口数据。过去10年的数据:

    年份人口(万)
    2009133126
    2010133770
    2011134413
    2012135069
    2013135738
    2014136427
    2015137122
    2016137866
    2017138639
    2018139538
    % 人口预测(注意:一般我们很少使用插值算法来预测数据,随着课程的深入,后面的章节会有更适合预测的算法供大家选择,例如灰色预测、拟合预测等)
    population=[133126,133770,134413,135069,135738,136427,137122,137866,138639, 139538];
    year = 2009:2018;
    p1 = pchip(year, population, 2019:2021)  %分段三次埃尔米特插值预测
    p2 = spline(year, population, 2019:2021) %三次样条插值预测
    figure(4);
    plot(year, population,'o',2019:2021,p1,'r*-',2019:2021,p2,'bx-')
    legend('样本点','三次埃尔米特插值预测','三次样条插值预测','Location','SouthEast')

     五、实例——数据预处理(补全数据)

    原始数据:需要插值得到偶数周轮虫、溶氧等的值。

    '周数'13579111315
    '轮虫'19131945192022052260230223852420
    '溶氧'5.123.26.723.362.44.146.434.6
    'COD'21.92026.827.7323.422.7525.3626.03
    '水温'24.825.726.82830.43027.630.8
    'PH值'9.319.149.149.299.229.339.169.26
    '盐度'1.82.31.92.12.11.11.51.5
    '透明度'2824262222201923
    '总碱度'425.11457.99492.48492.08501.93598.48604.44623.89
    '氯离子'628.1639.2648.87640.33616.43614.72507.14580
    '透明度'2824262222201923
    '生物量'30.5836.1949.7560.5856.5860.0667.9967.74

    Z.mat的数据

    13579111315
    19131945192022052260230223852420
    5.123.26.723.362.44.146.434.6
    21.92026.827.7323.422.7525.3626.03
    24.825.726.82830.43027.630.8
    9.319.149.149.299.229.339.169.26
    1.82.31.92.12.11.11.51.5
    2824262222201923
    425.11457.99492.48492.08501.93598.48604.44623.89
    628.1639.2648.87640.33616.43614.72507.14580
    2824262222201923
    30.5836.1949.7560.5856.5860.0667.9967.74

    【分析】在进行插值时,以周数为自变量,轮虫量、溶氧量等分别为因变量y1,y2……,使用循环进行插值。

    比如原始的自变量:x:1 3 5 7 9 11 13 15

    轮虫量y1:1913    1945    1920    2205    2260    2302    2385    2420

    需要插值的自变量:new_x,为偶数周,或者为了代码编写方便,可以写成1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

    插值后的轮虫量p1=spline(x,y,new_x)

    使用方法:subplot(m,n,p)或者subplot(m n p)。

    subplot是将多个图画到一个平面上的工具。其中,m表示是图排成m行,n表示图排成n列,也就是整个figure中有n个图是排成一行的,一共m行,如果m=2就是表示2行图。p表示图所在的位置,p=1表示从左到右从上到下的第一个位置。在第p块创建坐标系,并返回它的句柄。

    1. 三次埃米尔特插值

    %插值预测中间周的水体评价指标
    load Z.mat
    x=Z(1,:); %Z的第一行是星期Z: 1     3     5     7     9    11    13    15
    [n,m]=size(Z);%n为Z的行数,m为Z的列数
    % 注意Matlab的数组中不能保存字符串,如果要生成字符串数组,就需要使用元胞数组,其用大括号{}定义和引用
    ylab={'周数','轮虫','溶氧','COD','水温','PH值','盐度','透明度','总碱度','氯离子','透明度','生物量'};  % 等会要画的图形的标签
    disp(['共有' num2str(n-1) '个指标要进行插值。'])
    disp('正在对一号池三次埃尔米特插值,请等待')%一号池共有十一组要插值的数据,算上星期所在的第一行,共十二行
    P=zeros(11,15);%对要储存数据的矩阵P赋予初值
    for i=2:n%从第二行开始都是要进行插值的指标
        y=Z(i,:);%将每一行依次赋值给y
        new_x=1:15;%要进行插值的x
        p1=pchip(x,y,new_x);%调用三次埃尔米特插值函数
        subplot(4,3,i-1);%将所有图依次变现在4*3的一幅大图上
        plot(x,y,'ro',new_x,p1,'g-');%画出每次循环处理后的图像
        axis([0 15,-inf,inf])  %设置坐标轴的范围,这里设置横坐标轴0-15,纵坐标不变化
        %  xlabel('星期')%x轴标题
        ylabel(ylab{i})%y轴标题  这里是直接引用元胞数组中的字符串哦
        P(i-1,:)=p1;%将每次插值之后的结果保存在P矩阵中       
    end
    legend('原始数据','三次埃尔米特插值数据','Location','SouthEast')%加上标注,注意要手动在图中拖动标注到图片右下角哦
    P = [1:15; P]  %把P的第一行加上周数

    插值后的结果放在p矩阵中,记得保留2位小数就可以了。

    123456789101112131415
    1913.001936.561945.001932.501920.002050.972205.002238.072260.002279.982302.002344.322385.002407.282420.00
    5.123.583.204.966.725.233.362.692.403.024.145.536.436.004.60
    21.9020.2420.0023.2026.8027.4727.7325.7123.4022.9322.7523.9225.3625.8326.03
    24.8025.2325.7026.2326.8027.3428.0029.4030.4030.2930.0028.7127.6028.4530.80
    9.319.199.149.149.149.229.299.269.229.289.339.259.169.189.26
    1.802.172.302.101.902.002.102.102.101.601.101.301.501.501.50
    28.0025.1324.0025.0026.0024.0022.0022.0022.0021.1720.0019.3319.0020.1923.00
    425.11441.35457.99479.44492.48492.28492.08494.77501.93551.04598.48601.72604.44612.03623.89
    628.10633.83639.20645.33648.87646.17640.33627.21616.43615.60614.72560.51507.14523.19580.00
    28.0025.1324.0025.0026.0024.0022.0022.0022.0021.1720.0019.3319.0020.1923.00
    30.5832.6036.1942.4649.7556.6760.5858.5856.5857.7260.0664.6367.9967.9667.74

    2. 三次样条插值

    %插值预测中间周的水体评价指标
    load Z.mat
    x=Z(1,:); %Z的第一行是星期Z: 1     3     5     7     9    11    13    15
    [n,m]=size(Z);%n为Z的行数,m为Z的列数
    % 注意Matlab的数组中不能保存字符串,如果要生成字符串数组,就需要使用元胞数组,其用大括号{}定义和引用
    ylab={'周数','轮虫','溶氧','COD','水温','PH值','盐度','透明度','总碱度','氯离子','透明度','生物量'};  % 等会要画的图形的标签
    disp(['共有' num2str(n-1) '个指标要进行插值。'])
    disp('正在对一号池三次样条插值,请等待')%一号池共有十一组要插值的数据,算上星期所在的第一行,共十二行
    P=zeros(11,15);%对要储存数据的矩阵P赋予初值
    for i=2:n%从第二行开始都是要进行插值的指标
        y=Z(i,:);%将每一行依次赋值给y
        new_x=1:15;%要进行插值的x
        p1=spline(x,y,new_x);%调用三次样条插值函数
        subplot(4,3,i-1);%将所有图依次变现在4*3的一幅大图上
        plot(x,y,'ro',new_x,p1,'m-')%画出每次循环处理后的图像
        axis([0 15,-inf,inf])  %设置坐标轴的范围,这里设置横坐标轴0-15,纵坐标不变化
        %  xlabel('星期')%x轴标题
        ylabel(ylab{i})%y轴标题  这里是直接引用元胞数组中的字符串哦
        P(i-1,:)=p1;%将每次插值之后的结果保存在P矩阵中       
    end
    legend('原始数据','三次样条插值','Location','SouthEast')%加上标注,注意要手动在图中拖动标注到图片右下角哦
    P = [1:15; P]  %把P的第一行加上周数

    %插值预测中间周的水体评价指标
    load Z.mat
    x=Z(1,:); %Z的第一行是星期Z: 1     3     5     7     9    11    13    15
    [n,m]=size(Z);%n为Z的行数,m为Z的列数
    % 注意Matlab的数组中不能保存字符串,如果要生成字符串数组,就需要使用元胞数组,其用大括号{}定义和引用
    ylab={'周数','轮虫','溶氧','COD','水温','PH值','盐度','透明度','总碱度','氯离子','透明度','生物量'};  % 等会要画的图形的标签
    disp(['共有' num2str(n-1) '个指标要进行插值。'])
    disp('正在对一号池三次样条插值,请等待')%一号池共有十一组要插值的数据,算上星期所在的第一行,共十二行
    %P=zeros(11,15);%对要储存数据的矩阵P赋予初值
    for i=2:n%从第二行开始都是要进行插值的指标
        y=Z(i,:);%将每一行依次赋值给y
        new_x=2:2:15%要进行插值的x,只插值偶数周
        p1=spline(x,y,new_x);%调用三次样条插值函数
        subplot(3,4,i-1);%将所有图依次变现在3*4的一幅大图上
        plot(x,y,'ro',new_x,p1,'m*')%画出每次循环处理后的图像
        axis([0 15,-inf,inf])  %设置坐标轴的范围,这里设置横坐标轴0-15,纵坐标不变化
        %  xlabel('星期')%x轴标题
        ylabel(ylab{i})%y轴标题  这里是直接引用元胞数组中的字符串哦
        %P(i-1,:)=p1;%将每次插值之后的结果保存在P矩阵中       
    end
    legend('原始数据','三次样条插值','Location','SouthEast')%加上标注,注意要手动在图中拖动标注到图片右下角哦
    %P = [1:15; P]  %把P的第一行加上周数

     

    3. 三次样条插值和三次埃尔米特插值的对比

    %插值预测中间周的水体评价指标
    load Z.mat
    x=Z(1,:); %Z的第一行是星期Z: 1     3     5     7     9    11    13    15
    [n,m]=size(Z);%n为Z的行数,m为Z的列数
    % 注意Matlab的数组中不能保存字符串,如果要生成字符串数组,就需要使用元胞数组,其用大括号{}定义和引用
    ylab={'周数','轮虫','溶氧','COD','水温','PH值','盐度','透明度','总碱度','氯离子','透明度','生物量'};  % 等会要画的图形的标签
    disp(['共有' num2str(n-1) '个指标要进行插值。'])
    P=zeros(11,15);%对要储存数据的矩阵P赋予初值
    for i=2:n%从第二行开始都是要进行插值的指标
        y=Z(i,:);%将每一行依次赋值给y
        new_x=1:15;%要进行插值的x
        p1=pchip(x,y,new_x);%调用三次埃尔米特插值函数
        p2=spline(x,y,new_x);%调用三次样条插值函数
        subplot(4,3,i-1);%将所有图依次变现在4*3的一幅大图上
        plot(x,y,'ro',new_x,p1,'g-',new_x,p2,'m-');%画出每次循环处理后的图像
        axis([0 15,-inf,inf])  %设置坐标轴的范围,这里设置横坐标轴0-15,纵坐标不变化
        %  xlabel('星期')%x轴标题
        ylabel(ylab{i})%y轴标题  这里是直接引用元胞数组中的字符串哦
        %P(i-1,:)=p1;%将每次插值之后的结果保存在P矩阵中       
    end
    legend('原始数据','三次埃尔米特插值函数','三次样条插值','Location','SouthEast')%加上标注,注意要手动在图中拖动标注到图片右下角哦
    %P = [1:15; P]  %把P的第一行加上周数

    最后插值的结果可以求两种插值方法计算得到的平均值。 

     

    展开全文
  • matlab克里格插值

    千次阅读 多人点赞 2020-04-18 22:45:59
    为履行前期承若,现在公开matlab的克里格插值的代码。 原理介绍 普通克里格法:假定Z(x)Z(x)Z(x)是满足本征假设的一个随机过程,该随机过程有nnn个观测值z(xi)z(xi)z(xi),要预测未采样点x0处的值,则线性预测值Z∗...

    前言

    为履行前期承诺,现在公开matlab的克里格插值的代码。

    原理介绍

    普通克里格法:假定 Z ( x ) Z(x) Z(x)是满足本征假设的一个随机过程,该随机过程有 n n n个观测值 z ( x i ) z(xi) z(xi),要预测未采样点x0处的值,则线性预测值 Z ∗ ( x 0 ) Z^{*}(x_0) Z(x0)可以表示如下:
    Z ∗ ( x 0 ) = ∑ i = 1 n λ i z ( x i ) Z^{*}(x_0)=\sum_{i=1}^n\lambda_iz(x_i) Z(x0)=i=1nλiz(xi)
    Kriging是在使预测无偏并有最小方差的基础上,去确定最优的权重值,满足以下两个条件:
    (1)无偏性条件: E [ Z ∗ ( x 0 ) − Z ( x ) ) ] = 0 E[Z^{*}(x_0)-Z(x_))]=0 E[Z(x0)Z(x))]=0
    (2)最优条件: v a r [ Z ∗ ( x 0 ) − Z ( x 0 ) ] = m i n var[Z^*(x_0)-Z(x_0)]=min var[Z(x0)Z(x0)]=min
    有如下
    E [ Z ∗ ( x 0 ) − Z ( x 0 ) ] = E [ ∑ i = 1 n λ i z ( x i ) − Z ( x 0 ) ] = ∑ i = 1 n λ i E [ Z ( x i ) ] − E [ Z ( x 0 ) ] E\left[Z^{*}\left(x_{0}\right)-Z\left(x_{0}\right)\right]=E\left[\sum_{i=1}^{n} \lambda_{i} z\left(x_{i}\right)-Z\left(x_{0}\right)\right]=\sum_{i=1}^{n} \lambda_{i} E\left[Z\left(x_{i}\right)\right]-E\left[Z\left(x_{0}\right)\right] E[Z(x0)Z(x0)]=E[i=1nλiz(xi)Z(x0)]=i=1nλiE[Z(xi)]E[Z(x0)]
    根据本征假设 E [ Z ( x i ) ] = E [ Z ( x 0 ) ] = m E[Z(x_i)]=E[Z(x_0)]=m E[Z(xi)]=E[Z(x0)]=m
    上式进一步表示:
    ∑ i = 1 n λ i E [ Z ( x i ) ] − E [ Z ( x 0 ) ] = ∑ i = 1 n λ i m − m = m ( ∑ i = 1 n λ i − 1 ) = 0 \sum_{i=1}^{n} \lambda_{i} E\left[Z\left(x_{i}\right)\right]-E\left[Z\left(x_{0}\right)\right]=\sum_{i=1}^{n} \lambda_{i} m-m=m\left(\sum_{i=1}^{n} \lambda_{i}-1\right)=0 i=1nλiE[Z(xi)]E[Z(x0)]=i=1nλimm=m(i=1nλi1)=0
    因此要满足如下条件:
    ∑ i = 1 n λ i = 1 \sum_{i=1}^n\lambda_i=1 i=1nλi=1
    普通克里格插值的条件如下:
    1.在研究区域内,区域化变量 Z ( x ) Z(x) Z(x)的增量的数学期望对任意 x x x h h h存在且等于0
    E [ Z ( x ) − Z ( x + h ) ] = 0 E[Z(x)-Z(x+h)]=0 E[Z(x)Z(x+h)]=0
    2.在研究区域内,区域化变量的增量[ Z ( x ) − Z ( x + h ) ] Z(x)-Z(x+h)] Z(x)Z(x+h)]的方差对任意x和h存在且平稳
    Var ⁡ [ Z ( x ) − Z ( x + h ) ] = E [ { Z ( x ) − Z ( x + h ) } 2 ] = 2 r ( h ) \operatorname{Var}[Z(x)-Z(x+h)]=E\left[\{Z(x)-Z(x+h)\}^{2}\right]=2 r(h) Var[Z(x)Z(x+h)]=E[{Z(x)Z(x+h)}2]=2r(h)

    在本证假设条件下
    var ⁡ [ Z ∗ ( x 0 ) − Z ( x 0 ) ] = E [ { Z ∗ ( x 0 ) − Z ( x 0 ) } 2 ] = E [ { ∑ i = 1 n λ i z ( x i ) − Z ( x 0 ) } 2 ] = 2 ∑ i = 0 n λ i γ ( x i , x 0 ) − ∑ i = 1 n ∑ j = 1 n λ i λ j γ ( x i , x j ) \begin{array}{l} \operatorname{var}\left[Z^{*}\left(x_{0}\right)-Z\left(x_{0}\right)\right]=E\left[\left\{Z *\left(x_{0}\right)-Z\left(x_{0}\right)\right\}^{2}\right]=E\left[\left\{\sum_{i=1}^{n} \lambda_{i} z\left(x_{i}\right)-Z\left(x_{0}\right)\right\}^{2}\right] \\ =2 \sum_{i=0}^{n} \lambda_{i} \gamma\left(x_{i}, x_{0}\right)-\sum_{i=1}^{n} \sum_{j=1}^{n} \lambda_{i} \lambda_{j} \gamma\left(x_{i}, x_{j}\right) \end{array} var[Z(x0)Z(x0)]=E[{Z(x0)Z(x0)}2]=E[{i=1nλiz(xi)Z(x0)}2]=2i=0nλiγ(xi,x0)i=1nj=1nλiλjγ(xi,xj)
    根据方差最小原则,借助拉格朗日乘子,普通克里格的预测方程组为
    { ∑ i = 1 n λ i = 1 ∑ i = 1 n λ i γ ( x i , x j ) + φ ( x 0 ) = γ ( x j , x 0 ) \left\{\begin{array}{l} \sum_{i=1}^{n} \lambda_{i}=1 \\ \sum_{i=1}^{n} \lambda_{i} \gamma\left(x_{i}, x_{j}\right)+\varphi\left(x_{0}\right)=\gamma\left(x_{j}, x_{0}\right) \end{array}\right. {i=1nλi=1i=1nλiγ(xi,xj)+φ(x0)=γ(xj,x0)
    预测方差为
    σ 2 ( x 0 ) = ∑ i = 1 n λ i γ ( x i , x 0 ) + ϕ \sigma^2(x_0)=\sum_{i=1}^n\lambda_i\gamma(x_i,x_0)+\phi σ2(x0)=i=1nλiγ(xi,x0)+ϕ

    克里格公式也可以用矩阵的形式表示,对点状克里格,有:
    A ∙ [ λ φ ] = B A \bullet\left[\begin{array}{l} \lambda \\ \varphi \end{array}\right]=B A[λφ]=B
    对应的矩阵如下:
    [ γ ( s 1 , s 1 ) γ ( s 1 , s 2 ) ⋯ γ ( s 1 , s n ) 1 γ ( s 2 , s 1 ) γ ( s 2 , s 2 ) ⋯ γ ( s 2 , s n ) 1 ⋮ ⋮ ⋯ ⋮ ⋮ γ ( s n , s 1 ) γ ( s n , p 2 ) ⋯ γ ( p n , s n ) 1 1 1 ⋯ 1 0 ] × [ λ 1 λ 2 ⋮ λ n μ ] = [ γ ( s 1 , s 0 ) γ ( s 2 , s 0 ) ⋮ γ ( s n , s 0 ) 1 ] \left[\begin{array}{ccccc} \gamma\left(s_{1}, s_{1}\right) & \gamma\left(s_{1}, s_{2}\right) & \cdots & \gamma\left(s_{1}, s_{n}\right) & 1 \\ \gamma\left(s_{2}, s_{1}\right) & \gamma\left(s_{2}, s_{2}\right) & \cdots & \gamma\left(s_{2}, s_{n}\right) & 1 \\ \vdots & \vdots & \cdots & \vdots & \vdots \\ \gamma\left(s_{n}, s_{1}\right) & \gamma\left(s_{n}, p_{2}\right) & \cdots & \gamma\left(p_{n}, s_{n}\right) & 1 \\ 1 & 1 & \cdots & 1 & 0 \end{array}\right] \times\left[\begin{array}{c} \lambda_{1} \\ \lambda_{2} \\ \vdots \\ \lambda_{n} \\ \mu \end{array}\right]=\left[\begin{array}{c} \gamma\left(s_{1}, s_{0}\right) \\ \gamma\left(s_{2}, s_{0}\right) \\ \vdots \\ \gamma\left(s_{n}, s_{0}\right) \\ 1 \end{array}\right] γ(s1,s1)γ(s2,s1)γ(sn,s1)1γ(s1,s2)γ(s2,s2)γ(sn,p2)1γ(s1,sn)γ(s2,sn)γ(pn,sn)11110×λ1λ2λnμ=γ(s1,s0)γ(s2,s0)γ(sn,s0)1

    计算实例

    以实验数据为例,对于不规则样点来进行克里格插值,首先要在预定的空间范围下设置网格,确定出左下角的xy起始坐标,将其进行坐标系平移,可以减少计算复杂度,将每个格网的中心点设置为待插值样点,以每个样点为待插值样点,搜索最近的 n n n个点,计算第 i i i行第 j j j列的格网的中心坐标是
    ( x m i n + g r i d s i z e / 2 + ( j − 1 ) ∗ g r i d s i z e , y m i n + g r i d s i z e / 2 + ( m − i ) ∗ g r i d s i z e ) (xmin+gridsize/2+(j-1)*gridsize,ymin+gridsize/2+(m-i)*gridsize) (xmin+gridsize/2+(j1)gridsize,ymin+gridsize/2+(mi)gridsize)
    基于我前面缩写mink函数搜索得到最近的 n n n个点,计算他们之间的距离代入到协方差函数之中
    得到每个点的 A A A B B B矩阵,基于矩阵左除的方法对齐进行求解 λ 、 μ \lambda、\mu λμ
    结果展示
    在这里插入图片描述
    最终插值效果如上图所示

    总结

    本次总结主要在代码的简介度较高,以较快的速度来实现比较简单的克里格插值的这一算法。本次代码均在附录中展示,供大家学习,从代码层次对算法进一步了解。

    附录

    其中covfun_automatch函数需要自己写,根据自己变异函数模型来写对应的协方差函数

    %输入插值格网大小girdsize,邻近点个数num_neighbor,插值的空间范围即沙洋县shp图层
    %输入模型参数result_coef,其中1,2,3,4分别为c0,c1,a,m,块金值,基态值,变程,模型参数(1为球状,2为指数,3为高斯)
    clear;
    clc;
    digits(20);   %设定默认的精度
    tic;
    result_coef=[0.0377,0.0814,42541,2];num_neighbor=4;gridsize=100;%输入参数
    data=xlsread('F:\geostatistical_test\test\data_and_materials\delete_repeat.xlsx');
    %data(:,1),data(:,2)对应于xy坐标,data(:,3)为对应的属性
    Map=shaperead('F:\geostatistical_test\test\data_and_materials\沙洋边界\sy_xzbj.shp'); %读取沙洋县边获取空间范围
    xyrange=[];
    for i=1:size(Map,1)
        xyrange=[xyrange;Map(i).BoundingBox];
    end
    maxxy=max(xyrange);minxy=min(xyrange);xylength=maxxy-minxy;%xylength是最小外接矩形的长x和宽y
    minx=minxy(1);miny=minxy(2);%取出最小xy坐标
    data(:,1)=data(:,1)-minx;data(:,2)=data(:,2)-miny;%这一步是让左下角的xy坐标为0,这样计算距离就比较快
    row_col=ceil(xylength./gridsize);%生成网格的大小,row_col(1)为所需列数,row_col(2)为所需行数
    xc=zeros(row_col(2),row_col(1));yc=zeros(size(xc));%xc,yc分别存放格网中心的坐标xy
    colmat=repmat(1:row_col(1),row_col(2),1);rowmat=repmat((1:row_col(2))',1,row_col(1));
    xc=gridsize/2+(colmat-1).*gridsize;yc=gridsize/2+(row_col(2)-rowmat).*gridsize;
    %这一行是将最小外接矩形的左下角xmin,ymin,与格网大小建立关系,在第i行第j列的格网的中心坐标是
    % xmin+gridsize/2+(j-1)*gridsize,ymin+gridsize/2+(m-i)*gridsize,因为matlab索引从左上角开始,因此,是m-i
    result=zeros(size(xc));sigma2=zeros(size(xc));
    for i=1:row_col(2)
        for j=1:row_col(1)
            xy=[xc(i,j),yc(i,j)];%对第i,j个点寻找最近距离的点
            %计算样点离xy距离平方
            D2=(data(:,1)-xy(1)).^2+(data(:,2)-xy(2)).^2;
            [D,tag]=mink(D2,num_neighbor);%tag表示最近的4个点D表示对应的距离平方
            Dc_neighbor=sqrt(D);
            tag_neighbor=tag;
            %取出在邻域内的样点序号(tag_neighbor)以及对应插值点离样点的距离Dc_neighbor
            D_neighbor=squareform(pdist(data(tag,1:2)));
            %这一步利用的是matlab索引均为矢量的优势,会对里面矢量每个元素穷尽组合来组成矩阵,这里正好
            %可以生成插值点邻域内采样点点对之间的距离矩阵D_neighbor
            covD_neighbor=covfun_automatch(D_neighbor,result_coef);%生成采样点间对应的协方差矩阵
            covDc_neighbor=covfun_automatch(Dc_neighbor,result_coef);%%生成采样点与插值点之间对应的协方差矩阵
            A=[covD_neighbor,ones(size(covD_neighbor,1),1)];
            A=[A;ones(1,size(A,2))];
            A(end,end)=0;
            B=[covDc_neighbor;1];
            lambdau=A\B;%lambdau矩阵得到的最终的各个点的权重与u值
            lambda=lambdau(1:end-1);
            u=lambdau(end);
            result(i,j)=sum(lambda.*data(tag_neighbor,3));
            %这里soil_properties(tag_neighbor-1)减1的原因是前面距离矩阵上是加上了插值的点,因此要获取原始的
            %在data数据里的索引只需要将每个tag_neighbor-1即可
            sigma2(i,j)=mean(diag(covD_neighbor))-sum(lambda.*covDc_neighbor)+u;
            %克里格方差计算公式
        end
    end
    t=toc%运行时间
    data_validij=[row_col(2)-round((data(:,2)-gridsize/2)./gridsize),round((data(:,1)-gridsize/2)./gridsize)+1];
    %求出验证点对应的在格网中的行列号i,j索引,第一列为i,第二列为j
    valid_predict_result=diag(result(data_validij(:,1),data_validij(:,2)));%这里仍然利用的是matlab行列索引均为矢量的优势,这样对角线上的元素就是我们想要找到一一对应的元素
    %因此要用到diag函数得到了验证点对应的插值结果
    valid_actual_value=data(:,3);%验证点的真实值
    [R,RMSE,erroravg]=test_kriging(valid_actual_value,valid_predict_result);
    %返回相关系数R,均方根误差RMSE,平均误差erroravg
    figure;
    imagesc([minx,minx+row_col(2)*gridsize],[miny,miny+row_col(1)*gridsize],result);colorbar;%定义xy范围画出工具条
    title('Ordinary Kriging interpolation','FontSize',20);xlabel('x(m)','FontSize',14);ylabel('y(m)','FontSize',14);
    set(gca,'Fontname','times new Roman');
    figure;
    imagesc([minx,minx+row_col(2)*gridsize],[miny,miny+row_col(1)*gridsize],sigma2);c1=colorbar;%定义xy范围画出工具条
    title('Kriging Estimated Variance','FontSize',20);xlabel('x(m)','FontSize',14);ylabel('y(m)','FontSize',14);
    set(gca,'Fontname','times new Roman');
    
    展开全文
  • 数据插值,拟合,方程求解的MATLAB实现一.简介二.数据插值1.实验数据的一维插值2.实验数据的二维插值三.数据拟合1.实验数据的拟合四.方程和方程组1.方程和方程组的求解 一.简介 在生产和科学实验中,自变量x与因...
  • 克里金插值matlab程序

    热门讨论 2012-09-10 22:08:17
    克里金插值matlab程序 克里金(Kriging)插值法又称空间自协方差最佳插值法,它是以南非矿业工程师D.G.Krige的名字命名的一种最优内插法。克里金法广泛地应用于地下水模拟、土壤制图等领域,是一种很有用的地质统计...
  • MATLAB插值问题

    2021-04-21 20:40:13
    一、一元函数插值已知函数y=f(x)在区间[a,b]上的n+1个不同点的函数值为,若存在一个简单函数F(x), 使,称F(x)为f(x)在区间[a,b]上的插值函数,称(xi, yi)为插值节点。若F(x)为多项式,称为多项式插值(或代数插值) ;...
  • matlab自带的插值函数interp1的几种插值方法

    万次阅读 多人点赞 2018-04-27 14:27:51
    插值法 插值法又称“内插法”,是利用函数f (x)在某区间中已知的若干点的函数值,作出适当的特定函数,在区间的其他点上用这特定函数的值作为函数f (x)的近似值,这种方法称为插值法。如果这特定函数是多项式,就称...
  • matlab插值函数

    2021-06-14 14:13:20
    插值 x=0:2*pi; y=sin(x); xx=0:0.5:2*pi; %interp1对sin函数进行分段线性插值,调用interp1的时候,默认的是分段线性插值 y1=interp1(x,y,xx); figure plot(x,y,'o',xx,y1,'r') title('分段线性插值') %临近...
  • 1. 拉格朗日多项式插值了解概念插值多项式插值节点范德蒙特(Vandermonde)行列式截断误差、插值余项特点函数实现function y=lagrange(x0,y0,x)n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=...
  • 系列文章目录 提示:这里可以添加系列文章的所有文章的...读入数据总结 一、pandas是什么? 实验1.2 误差传播与算法稳定性 问题提出: 考虑一个简单的由积分定义的序列 ∫01xnex−1dx,n=1,2,……\int_{0}^{1}x^ne
  • Matlab实现插值与拟合

    千次阅读 2018-01-25 23:10:01
    1、拉格朗日插值 ...%n 个节点数据以数组 x0, y0 输入(注意 Matlat 的数组下标从1开始), %m 个插值点以数组 x 输入,输出数组 y 为 m 个插值 n=length(x0);m=length(x); for i=1:m z=x(i); s
  • 一维插值 插值不同于拟合。插值函数经过样本点,拟合函数一般基于最小二乘法尽量靠近所有样本点穿过。常见插值方法有拉格朗日插值法、分段插值法、样条插值法。 拉格朗日插值多项式:当节点数n较大时,拉格朗日插值...
  • Matlab插值方法
  • 1、拉格朗日插值新建如下函数:function y=lagrange(x0,y0,x)%拉格朗日插值函数%n 个节点数据以数组 x0, y0 输入(注意 Matlat 的数组下标从1开始),%m 个插值点以数组 x 输入,输出数组 y 为 m 个插值n=length(x0);...
  • 实验平台与数据本算法使用Matlab语言实现,实验平台为Windows 8 32位操作系统、4GB内存(可用为2.31GB)、Matlab2013b。数据1: 大小为:256*256 的lena灰度图像,将使用实现的算法对其进行2倍放大操作,如下图1所示...
  • MATLAB回归、插值、逼近、拟合总结

    万次阅读 多人点赞 2017-06-29 11:11:58
    2、多项式插值:用一个多项式来近似代替数据列表函数,并要求多项式通过列表函数中给定的数据点。(插值曲线要经过型值点。) 3、多项式逼近:为复杂函数寻找近似替代多项式函数,其误差在某种度量意义下最小。...
  • 原标题:简单调研多维插值方法■2020-11-07 11:36:50以前的时候用过二维的插值, 见 二维三次卷积插值算法及Fortran代码 [1], 也用过matlab自带的插值方法, 见 Matlab插值曲面及其曲率的计算 [2]. 最近因为需要自己...
  • 利用matlab的interp1()对矩阵进行插值

    万次阅读 2018-10-21 12:35:54
    利用matlab的interp1()对矩阵进行插值 对矩阵来说,所谓插值其实就是对矩阵...一维数据插值函数interp1()的用法: interp1(X,Y,X1,method) 参数说明 其中: X为原数组x坐标,对矩阵来说,x可以使用索引。 ...
  • 一,直接献上克里金插值MATLAB工具箱 链接:https://pan.baidu.com/s/1p0lt2G8KQ-els7Dyw5jPFg 提取码:wcss 下载后将该程序添加到MATLAB安装文件夹中的/toolbox下,然后在设置路径中添加该程序文件为默认路径即可...
  • 克里金插值(arcgis克里金插值步骤)

    千次阅读 2021-04-24 16:15:36
    1. 克里格方法概述 克里格方法(Kriging)又称空间局部插值法,是以变异函数理论和结构分析为基础, 在有限区域内对区域化变量进行无偏最优估计的一种方法,是地.克里金差值最后的出来的克里金误差有什么意义?是说...
  • 今天我们来谈谈数据插值以及它在matlab中的应用。 实际问题中观察到的数据点,譬如不同时间段的温度,湿度,在力的作用下物体的变形程度等,常常是客观的,真实发生的数据,具有客观世界的准确性,然而一般情况下...
  • variogram 计算各个维度的各向同性和各向异性实验变异函数。... 虽然这对于小数据集(n<2000)来说很快,但由于内存限制,它可能会在大数据集上失败。 因此,对于大型数据集,您可能需要使用 'subsample' 选项。
  • 较大剂量会增加INR,较剂量会降低INR。患者由护士定期监测,当他们的INR超出目标范围时,他们的剂量和测试频率会发生变化。 该文件INR.mat包含在五年内对患者进行的INR测量。该文件包括一个datetime数组,其中...
  • 克里金(Kriging)插值法又称空间自协方差最佳插值法,它是以南非矿业工程师D.G.Krige的名字命名的一种最优内插法。克里金法广泛地应用于地下水模拟、土壤制图等领域,是一种很有用的地质统计格网化方法。它首先考虑...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 2,838
精华内容 1,135
关键字:

matlab数据插值变小

matlab 订阅