全面地介绍 将二维数据矩阵插值拟合为二维曲面方程——基于matlab

前言

要利用matlab将二维数据矩阵拟合成一个二维曲面方程，首先，我们得明白什么是插值和拟合。下面，从对洋流数据的处理入手，先以一个基本的流程为例，进行问题的概述，然后介绍两种具体的方法。

原理初析

插值

根据已知的数据点，生成一个连续的函数，可以在任意位置求值。matlab 提供了多种插值方法，如线性插值、$B$-样条插值、多项式插值等。

这些方法的优缺点和适用场合各有不同，一般来说，线性插值比较简单和快速，但可能不够平滑；$B$ -样条插值比较平滑和灵活，但可能不够精确；多项式插值比较精确和光滑，但可能出现奇异或震荡现象。具体选择哪种方法，要根据数据的特点和目的来决定。

matlab中，可以使用 $interp2$ 函数或 $griddedInterpolant$ 类来进行二维插值，它们都可以指定不同的插值方法。例如，如果已知的数据矩阵是$V$，想要在 $Xq$ 和 $Yq$ 的网格上进行插值，可以使用以下代码：

% 线性插值
Vq1 = interp2(V,Xq,Yq,'linear');
% B-样条插值
Vq2 = interp2(V,Xq,Yq,'spline');
% 多项式插值
Vq3 = interp2(V,Xq,Yq,'pchip');

拟合

根据插值得到的函数，寻找一个合适的数学模型，可以用一个简洁的表达式来近似描述数据的规律。

matlab提供了多种拟合方法，如线性回归、非线性回归、最小二乘法等。这些方法的原理和应用也各有不同，一般来说，线性回归比较简单和稳定，但可能不够灵活；非线性回归比较灵活和通用，但可能不够稳定；最小二乘法比较优化和精确，但可能不够简洁。具体选择哪种方法，要根据数据的分布和形状来决定。matlab中，可以使用 $fit$ 函数或 $fittype$ 类来进行二维拟合，它们都可以指定不同的拟合方法。例如，如果插值得到的函数是 $Vq$，想要用一个二次多项式来拟合，可以使用以下代码：

% 线性回归
ft1 = fittype('poly22'); % 二次多项式
f1 = fit([Xq(:),Yq(:)],Vq1(:),ft1); % 拟合
% 非线性回归
ft2 = fittype('a*exp(b*x)*sin(c*y)'); % 非线性函数
f2 = fit([Xq(:),Yq(:)],Vq2(:),ft2); % 拟合
% 最小二乘法
ft3 = fittype('a*x^2+b*y^2+c*x*y+d*x+e*y+f','problem','f'); % 二次多项式，f是固定参数
f3 = fit([Xq(:),Yq(:)],Vq3(:),ft3,'problem',0); % 拟合

拟合完成后，可以用 f1，f2，f3等变量来表示拟合得到的二维曲面方程，也可以用 $coeffvalues$ 函数来提取方程的系数。

例如，如果想要得到 f1的表达式，可以使用以下代码：

% 提取系数
c1 = coeffvalues(f1);
% 构造表达式
s1 = sprintf('V = %.4f*X^2 + %.4f*Y^2 + %.4f*X*Y + %.4f*X + %.4f*Y + %.4f',c1(1),c1(2),c1(3),c1(4),c1(5),c1(6));

拟合优度的评估

对于方程的拟合，如果要做可行性分析的话，那必然要对拟合优度进行评估。

拟合优度是一种衡量拟合结果与实际数据之间的差异的指标。拟合优度越高，说明拟合结果越能反映数据的真实情况。拟合优度的评估方法有很多种，常见的有残差分析、决定系数、均方根误差等。不同的评估方法有不同的含义和应用场景，比如残差分析可以检验拟合函数的合理性，决定系数可以衡量拟合函数对数据变化的解释程度，均方根误差可以衡量拟合函数的预测精度等。

具体实例的展开

下面，将以两种方法进行具体问题的求解。

B-样条插值+最小二乘法

B-样条插值法

$B$-样条（B-spline）是一种基于多项式的插值方法，它可以构造出平滑且连续的曲线或曲面。b-样条的优点是可以灵活地控制插值的阶数、节点、边界条件等，从而得到不同的插值效果。b-样条的缺点是计算量较大，且难以表示复杂的形状。

$B$ -样条的数学定义如下：

设有 $n+1$ 个数据点 $\left(x_i,y_i\right)$，其中 $i=0,1,\dots ,n$ ，且 $x_0<x_1<\cdots <x_n$ 。$B$-样条插值法的目标是找到一个函数 $f\left(x\right)$，使得 $f\left(x_i\right)=y_i$，对所有 $i$ 成立，且 $f\left(x\right)$ 在每个区间 $\left[x_i,x_{i+1}\right]$ 上是一个 $k$ 阶多项式，其中 $k$ 是给定的插值阶数，通常为 $3$ 或 $4$。

为了构造这样的函数，我们需要引入一个节点序列 $t_0,t_1,\dots ,t_{n+k}$ ，其中 $t_i$ 是一个实数，且满足 $t_0\le t_1\le \cdots \le t_{n+k}$。节点序列的作用是划分数据点的区间，并确定每个区间上的多项式的系数。根据节点序列，我们可以定义一组基函数 $B_{i,k}(x)$ ，其中 $i=0,1,\dots ,n，k=1,2,\dots$，满足以下递推公式：

$$B_{i,1}(x)=\{\begin{array}{ll}1,& t_i\le x<t_{i+1}\\ 0,& \text{otherwise}\end{array}$$
$$B_{i,k}(x)=\frac{x-t_i}{t_{i+k-1}-t_i}{B}_{i,k-1}(x)+\frac{t_{i+k}-x}{t_{i+k}-t_{i+1}}{B}_{i+1,k-1}(x)$$

基函数 $B_{i,k}(x)$ 的意义是在节点序列 $t_i,t_{i+1},\dots ,$ 上的 $k$ 阶 $B$ -样条，它是一个非负的、分段的、最高为 $k-1$ 次的多项式，且在 $[t_i,t_{i+1}]$上为 $1$，在 $[t_{i+k-1},t_{i+k}]$ 上为 $0$，在其他区间上线性变化。

基于基函数，我们可以将插值函数 $f(x)$ 表示为：

$$f(x)=\sum _{i=0}^n{a}_i{B}_{i,k}(x)$$
其中 $a_i$ 是待求的系数，它们可以通过解线性方程组 $$\sum _{i=0}^n{a}_i{B}_{i,k}(x_j)=y_j$$

对 $j=0,1,\dots ,n$ 得到。

这样，我们就得到了一个 $B$ -样条插值函数 $f(x)$ ，它可以在任意给定的 $x$ 值处计算出对应的 $y$ 值。

如果我们要得到一个二维曲面方程 $f(x,y)$ ，我们可以先沿着 $x$ 方向对每个 $y$ 值进行 $B$ -样条插值，得到一组函数

$$f_i(y)=\sum _{j=0}^n{a}_{ij}{B}_{j,k}(x)$$

其中 $i=0,1,\dots ,m，$然后再沿着 $y$ 方向对每个 $x$ 值进行 $B$-样条插值，得到最终的函数

$$f(x,y)=\sum _{i=0}^m{f}_i(y)B_{i,k}(y)$$

这样，我们就得到了一个双向的 $B$ -样条插值曲面，它可以在任意给定的 $(x,y)$ 值处计算出对应的 $z$ 值。

lsqcurvefit 拟合法

接下来，介绍一种主要的最小二乘法的拟合算法在matlab里的应用——$lsqcurvefit$ 拟合法。

$lsqcurvefit$ 是 MATLAB 中的一个函数，它可以用来对数据进行非线性曲线拟合。$lsqcurvefit$ 的优点是可以灵活地指定拟合函数的形式，且可以考虑数据的误差和噪声。$lsqcurvefit$ 的缺点是需要提供拟合函数的初始参数，且可能存在局部最优解的问题。

设有 $n$ 个数据点 $(x_i,y_i)$ ，其中 $i=1,2,\dots ,n$ ，且 $x_i$ 和 $y_i$ 可以是标量或向量。$lsqcurvefit$ 的目标是找到一个函数 $f(x,a)$ ，使得

$$\sum _{i=1}^n\Vert y_i-f(x_i,a){\Vert }^2$$

最小，其中 $a$ 是待求的参数向量，它可以有一些约束条件，比如上下界、线性或非线性不等式或等式。

为了求解这样的问题，$lsqcurvefit$ 使用了一种基于信赖域的迭代算法。

它的基本思想是在每一步，以当前的参数值为中心，构造一个信赖域，然后在信赖域内用一个二次模型来近似目标函数，求解该二次模型的最优解，作为下一步的参数值。这个过程重复进行，直到满足收敛条件或达到最大迭代次数。信赖域的大小会根据每一步的优化效果进行调整，以保证算法的稳定性和效率。

$lsqcurvefit$ 在具体编程的时候需要提供以下输入：

• fun：一个函数句柄，表示拟合函数 $f(x,a)$ ，它接受两个输入参数，分别是自变量 $x$ 和参数向量 $a$ ，并返回一个输出值，表示因变量 $y$ 。
• x0：一个向量，表示参数向量 $a$的初始值，它决定了算法的起始点，对于非线性问题，不同的初始值可能导致不同的最优解，所以应该尽量选择一个合理的初始值。
• xdata：一个向量或矩阵，表示自变量 $x$ 的数据，它的大小应该与因变量 $y$ 的数据相匹配。
• ydata：一个向量或矩阵，表示因变量 $y$ 的数据，它的大小应该与自变量 $x$ 的数据相匹配。
• lb：一个向量，表示参数向量 $a$ 的下界，它的大小应该与参数向量 $a$ 相同，如果没有下界，可以用 Inf 表示。
• ub：一个向量，表示参数向量 $a$ 的上界，它的大小应该与参数向量 $a$ 相同，如果没有上界，可以用 Inf 表示。

$lsqcurvefit$ 的输出是一个向量，表示参数向量 $a$ 的最优值，它使得目标函数的平方和达到最小。

如果我们要得到一个二维曲面方程 $f(x,y,a)$ ，我们可以先将自变量 $x$ 和 $y$ 组合成一个矩阵，然后将其作为 xdata 输入 $lsqcurvefit$ ，得到参数向量 $a$ 的最优值，然后将其代入拟合函数 $f(x,y,a)$ ，得到最终的函数 $f(x,y)$ 。这样，我们就得到了一个非线性的曲面拟合方程，它可以在任意给定的 $(x,y)$ 值处计算出对应的 $z$ 值。

数学实现

为了演示插值和拟合的数学实现，使用一个示例数据矩阵，它包含了 $x$ ，$y$ 位置方向的，$u$ ，$v$方向速度的洋流数据，大小为 $50\times 50$ 。

这里，假设你的洋流数据是从官方网站下载的，所以我的这个算法就假设你直接 load 这个洋流数据文件。

读取数据矩阵，并提取 $x，y，u，v$ 四个变量，分别表示位置和速度。这里我假设数据矩阵已经保存在一个名为 data.mat 的文件中，各位小伙伴们可以根据自己的实际情况修改文件名和路径哒。

load data.mat
x = data(:,1);
y = data(:,2);
u = data(:,3);
v = data(:,4);

将 $x，y，u，v$ 四个变量转换为矩阵形式，以便进行插值和拟合。这里我假设 $x，y$ 是均匀分布的，当然也可以根据自己的实际情况修改矩阵的大小和范围。

X = reshape(x,50,50);
Y = reshape(y,50,50);
U = reshape(u,50,50);
V = reshape(v,50,50);

计算 $z$ 变量，表示速度的大小，它是 $u，v$ 的平方和的开方。这个变量将作为插值和拟合的因变量，表示洋流的强度。

Z = sqrt(U.^2+V.^2);

使用 $B$ -样条插值法，对 $z$ 变量进行插值，得到一个二维曲面方程 $f_1(x,y)$ 。

这里，我使用 MATLAB 的 $CurveFittingToolbox$ 中的 $fit$ 函数，指定插值方法为 $‘biharmonicinterp’$，表示双向的三次 $B$ -样条插值。也可以使用其他的插值方法，比如 $‘thinplateinterp’$ ，表示双向的二次 $B$ -样条插值，或者 $‘linearinterp’$，表示双向的线性插值。

$fit$ 函数的输出是一个 $cfit$ 对象，它包含了插值函数的信息，也可以使用 $feval$ 函数来计算插值函数在任意 $(x,y)$ 值处的 $z$ 值，或者使用 $formula$ 函数来查看插值函数的具体表达式。

注意，由于 b-样条插值函数是分段的，它的表达式可能很长，且难以理解，所以我不建议各位友友直接查看它的表达式，而是使用 $feval$ 函数来计算它的值。

f1 = fit([X(:),Y(:)],Z(:),'biharmonicinterp');
Z1 = feval(f1,X,Y);

接下来，使用 $lsqcurvefit$ 拟合法，对 $z$ 变量进行拟合，得到一个二维曲面方程 $f_2(x,y,a)$。这里我假设拟合函数的形式为 $$f_2(x,y,a)=a_1\sin (a_{2}x+a_{3}y)+a_4\cos (a_{5}x+a_{6}y)+a_7$$

，其中 $a=(a_1,a_2,a_3,a_4,a_5,a_6,a_7)$ 是待求的参数向量。

当然，你也可以根据自己的实际情况修改拟合函数的形式，但请注意，拟合函数应该是非线性的，且能够较好地拟合数据的特征。

$lsqcurvefit$ 函数的输入包括拟合函数的函数句柄，参数向量的初始值，自变量和因变量的数据，参数向量的上下界等。$lsqcurvefit$ 函数的输出是参数向量的最优值，它使得目标函数的平方和达到最小。

根据参数向量的最优值，我们可以得到拟合函数的具体表达式，也可以使用 $feval$ 函数来计算拟合函数在任意 $(x,y)$ 值处的 $z$ 值。

fun = @(a,xdata) a(1)*sin(a(2)*xdata(:,1)+a(3)*xdata(:,2))+a(4)*cos(a(5)*xdata(:,1)+a(6)*xdata(:,2))+a(7);
x0 = [1,1,1,1,1,1,1];
xdata = [X(:),Y(:)];
ydata = Z(:);
lb = [-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf];
ub = [Inf,Inf,Inf,Inf,Inf,Inf,Inf];
a = lsqcurvefit(fun,x0,xdata,ydata,lb,ub);
Z2 = feval(fun,a,xdata);
Z2 = reshape(Z2,50,50);
f2 = @(x,y) a(1)*sin(a(2)*x+a(3)*y)+a(4)*cos(a(5)*x+a(6)*y)+a(7);

代码

接下来，是一个整个的代码实现。

% Load data matrix and extract variables
load data.mat
x = data(:,1);
y = data(:,2);
u = data(:,3);
v = data(:,4);

% Convert variables to matrix form
X = reshape(x,50,50);
Y = reshape(y,50,50);
U = reshape(u,50,50);
V = reshape(v,50,50);



% Calculate z variable as the magnitude of velocity 
Z = sqrt(U^2+V^2);
% Use b-spline interpolation method to interpolate z variable and get a surface equation f1(x,y) 
f1 = fit([X(:),Y(:)],Z(:),'biharmonicinterp'); Z1 = feval(f1,X,Y);
% Use lsqcurvefit fitting method to fit z variable and get a surface equation f2(x,y,a) 
fun = @(a,xdata) a(1)*sin(a(2)*xdata(:,1)+a(3)*xdata(:,2))+a(4)*cos(a(5)*xdata(:,1)+a(6)*xdata(:,2))+a(7); 
x0 = [1,1,1,1,1,1,1]; 
xdata = [X(:),Y(:)]; 
ydata = Z(:); 
lb = [-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf]; 
ub = [Inf,Inf,Inf,Inf,Inf,Inf,Inf]; 
a = lsqcurvefit(fun,x0,xdata,ydata,lb,ub); 
Z2 = feval(fun,a,xdata); 
Z2 = reshape(Z2,50,50); 
f2 = @(x,y) a(1)*sin(a(2)*x+a(3)*y)+a(4)*cos(a(5)*x+a(6)*y)+a(7);
% Plot the original data, the interpolated surface, and the fitted surface 
figure 
subplot(1,3,1) 
surf(X,Y,Z) 
title('Original Data') 
xlabel('x'),ylabel('y'),zlabel('z') 
subplot(1,3,2) 
surf(X,Y,Z1) 
title('Interpolated Surface') 
xlabel('x'), ylabel('y'),zlabel('z') 
subplot(1,3,3) 
surf(X,Y,Z2) 
title('Fitted Surface') 
xlabel('x'),ylabel('y'),zlabel('z')
% Compare the interpolation and fitting results with the original data % Calculate the mean squared error (MSE) and the coefficient of determination (R-squared) for each method 
mse1 = mean((Z(:)-Z1(:)).^2); 
mse2 = mean((Z(:)-Z2(:)).^2); 
r1 = 1 - sum((Z(:)-Z1(:)).^2)/sum((Z(:)-mean(Z(:))).^2); 
r2 = 1 - sum((Z(:)-Z2(:))^.2)/sum((Z(:)-mean(Z(:)))^.2); 
fprintf('The MSE and R-squared for b-spline interpolation are: %.4f and %.4f\n',mse1,r1); 
fprintf('The MSE and R-squared for lsqcurvefit fitting are: %.4f and %.4f\n',mse2,r2);
% Display the interpolation and fitting functions 
fprintf('The interpolation function is: f1(x,y) = %s\n',formula(f1));
fprintf('The fitting function is: f2(x,y) = %s\n',formula(f2));

B-样条插值+非线性回归法

这个方法就不着重讲原理分析了，讲一下拟合评估怎么做，还附带讲一讲洋流的可视化。

数据的预处理

同样，假设你已经有了一个 $50\times 50$ 的二维数据矩阵，比如洋流数据矩阵，其中每个元素表示某个位置的洋流速度。为了进行插值和拟合，你需要将这个矩阵转换为三个向量，分别表示 $x$ 坐标、$y$ 坐标和 $z$ 值。可以使用 matlab 的 $reshape$ 函数来实现这一转换:

% 假设 Z 是 50 x 50 的洋流数据矩阵
n = 50; % 矩阵的大小
x = repmat(1:n,1,n); % 重复 1 到 n，共 n 次，得到 x 坐标向量
y = reshape(repmat(1:n,n,1),1,n^2); % 重复 1 到 n，每个 n 次，然后转换为行向量，得到 y 坐标向量
z = reshape(Z,1,n^2); % 将 Z 转换为行向量，得到 z 值向量

% 假设 x, y, z 是已经转换好的三个向量，表示数据点的坐标和值
m = 100; % 网格的大小
xi = linspace(min(x),max(x),m); % 生成 x 方向的网格坐标
yi = linspace(min(y),max(y),m); % 生成 y 方向的网格坐标
[XI,YI] = meshgrid(xi,yi); % 生成网格矩阵
ZI = griddata(x,y,z,XI,YI,'cubic'); % 使用样条插值得到网格上的值

非线性回归的实现

为了使用 $fit$ 函数进行非线性回归，你需要指定拟合函数的形式和参数的初始值或范围。拟合函数的形式可以是一个匿名函数，也可以是一个字符串，表示函数的表达式。参数的初始值或范围可以是一个向量，也可以是一个矩阵，表示参数的下限和上限：

% 假设 x, y, z 是已经转换好的三个向量，表示数据点的坐标和值
f = @(a,b,c,d,x,y) a*x.^2 + b*y.^2 + c*x.*y + d; % 定义拟合函数的形式，这里是一个二次曲面
s = [1 1 1 1]; % 定义参数的初始值
f = fit([x,y],z,f,'StartPoint',s); % 使用非线性回归得到拟合结果

拟合优度的评估

前面已提及了拟合优度评估的原理，这里讲讲代码的实现。为了使用 $fit$ 函数的输出来评估拟合优度，你可以查看拟合结果的属性，其中包含了一些拟合优度的指标：

% 假设 f 是拟合结果
R2 = f.rsquare; % 得到决定系数
RMSE = f.rmse; % 得到均方根误差
disp(['R-squared = ',num2str(R2)]); % 显示决定系数
disp(['RMSE = ',num2str(RMSE)]); % 显示均方根误差

曲面方程的求解

为了使用 $fit$ 函数的输出来求解曲面方程，你可以使用 $f$ 函数来计算任意位置的 $z$ 值，或者使用 $feval$ 函数来计算网格上的 $z$ 值:

% 假设 f 是拟合结果，XI, YI 是网格矩阵
z1 = f(10,20); % 计算 x = 10, y = 20 时的 z 值
z2 = feval(f,XI,YI); % 计算网格上的 z 值
disp(['z(10,20) = ',num2str(z1)]); % 显示 z(10,20) 的值

矢量场的绘制

对于洋流的可视化，因为洋流由 $x,y,u,v$ 这四组向量组成，其中前两个定位洋流的每个坐标点的位置，后两个代表每个坐标点水平方向和竖直方向的速度，所以可以用矢量场来可视化，为了绘制矢量场，可以使用 $quiver$ 函数来绘制箭头，表示矢量的方向和大小。

在绘制的时候，需要指定箭头的起点和终点的坐标，以及箭头的颜色、线宽、比例等属性：

% 假设 XI, YI, ZI 是网格矩阵，表示插值后的数据
[U,V,W] = gradient(ZI); % 计算网格上的梯度，表示矢量的方向和大小
quiver3(XI,YI,ZI,U,V,W,'r','LineWidth',2,'AutoScaleFactor',0.5); % 绘制矢量场，使用红色的箭头，线宽为 2 ，比例因子为 0.5
xlabel('x'); % 标注 x 轴
ylabel('y'); % 标注 y 轴
zlabel('z'); % 标注 z 轴
title('Vector field of the fitted surface'); % 标注标题

尾声

本文由一个通用方法进行问题的介绍，然后，使用了两种组合方法进行具体问题求解的展开。

将二维数据矩阵转化为二维曲面方程有很多方法，可能本文的解答有局限，请多多见谅。

最后，祝君功不唐捐，玉汝于成！

参考

1. MathWorks. (2021). Interpolate scattered data.
2. MathWorks. (2021). Representing a matrix as a surface.
3. Wikipedia contributors. (2024, March 2). Bilinear interpolation. In Wikipedia, The Free Encyclopedia. Retrieved 02:41, March 4, 2024.
4. Kiusalaas, J. (2013). Interpolation in Two and Three Dimensions. In Numerical Methods in Engineering with Python 3 (pp. 97-122). Cambridge: Cambridge University Press. doi:10.1017/CBO9781139108278.006
5. de Boor, C. (1978). Biharmonic spline interpolation. In A Practical Guide to Splines (pp. 211-247). New York: Springer-Verlag.
6. MathWorks. (n.d.). Two-Dimensional Interpolation. Retrieved March 4, 2024.
7. XlXtrFun. (n.d.). Excel 2D interpolation function. Retrieved March 4, 2024.