基于广义逆矩阵理论的分段线性拟合方法

在采用最小二乘法解决某些拟合问题时，有时候通过简单的线性拟合并不能取得较为满意的结果，此时则可能需要采用分段线性拟合来解决该问题。以疲劳领域中的裂纹扩展问题为例，图2给出了7075-T7351铝合金板的裂纹扩展速率da/dN与裂纹尖端的有效应力强度因子幅值ΔK_eff之间的关系。

图1 7075-T7351铝合金板裂纹扩展速率曲线

可以看到，裂纹扩展速率与裂纹尖端的有效应力强度因子幅值并没有严格满足对数线性关系。但如果将其划分为多个区段，通过合理选择相邻区段之间的断点，则在每一个区段里面可以认为其近似满足对数线性关系。因此，可以采用分段线性拟合方法来拟合图1所示的裂纹扩展速率数据。

虽然我们可以将数据划分若干个区段，然后针对每个区段的数据分别进行线性拟合，但这样的方法并不能保证相邻拟合直线在断点处的连续性。为此，本文基于广义逆矩阵理论，提出了一种分段线性拟合方法，仅需指定待拟合数据点以及相邻区段之间的断点，即可自动计算分段拟合直线的参数，并保证拟合直线之间的连续性。

考虑如图2所示的分段线性拟合问题，图中黑色的实心圆点为待拟合的数据点，黑色直线为位于不同分段的拟合直线，红色的实心圆点为不同拟合分段之间的断点。

图2 分段线性拟合示意图

设第i个待拟合数据点的坐标为(x_i ,y_i) (i=1,…,N)，其中N为待拟合数据点的总数。并且，设第j个拟合区段（j=1,…,M; s_j-1≤x<s_j）的拟合直线的表达式为：

式中：α_j和β_j为第j个区段的拟合直线的系数，M为拟合区段的总数。s_j为第j个拟合区段的断点，为避免出现数据外插，这里取s₀=min{x_i; i=1,…,N}，s_M=max{x_i; i=1,…,N}。

假定所有的拟合数据点均位于拟合直线上，则必然满足：

注意上式中下标j的取值取决于待拟合数据点x_i所在的拟合区段。

将所有分段的拟合直线的系数表示为如下所示的系数向量：

则有：

式中：[y]为待拟合数据点的y坐标构成的向量，[A]为表征等价关系的系数矩阵，可分别表示为：

需要注意的是，{Coeff}中的系数并不是独立的，这是因为相邻区段的拟合直线在断点处必然满足连续性条件，这要求：

对于第一个连续性条件（j=1），我们有：

将其整理后可得：

同理，根据其他连续性条件我们有：

式中：

这里，我们将[C_j]称为第j个连续性条件矩阵。

通过第一个连续性条件矩阵，我们可以消去第一个区段的非独立系数α₁，即有：

式中：

式中：[O₁]为行数为1，列数为{2M-4}的零矩阵；[I₁]为行数和列数均为{2M-1}的单位矩阵。

同理，根据第2个区段的连续性条件，我们有：

式中：

式中：[O₂]是行数为1，列数为{2M-6}的零矩阵；[I₂]为行数和列数均为{2M-3}的单位矩阵。

再次根据第3个区段的连续性条件，我们有：

式中：

式中：[I_3,U]是一个行数和列数均为2的单位矩阵；[O₃]是一个行数为1，列数为{2M-8}的零矩阵；[I_3,L]是一个行数和列数均为{2M-5}的单位矩阵。

不难发现，如果继续使用其余分段的连续性条件，并重复上述过程，最终可以得到：

总结前面的规律可知，上式中[T_j]_j≥2的定义为：

式中：[I_j,U]为行数和列数均为{j-1}的单位矩阵；[O_j]为行数为1，列数为{2(M-j-1)}的零矩阵；[I_j,L]为行数和列数均为{2M-2j+1}的单位矩阵。

令

则有

可以看到，矩阵[Π_T]定义了独立变量{CoeffI}与非独立变量{Coeff}之间的变换关系。将上式代入到原始的线性方程组中，可得：

如果令

则有

由于拟合区段数量M在一般情况下远小于待拟合数据点数量N，因此上式是一个矛盾线性方程组。但根据广义逆矩阵理论，可以得到上述方程组的最小二乘解为：

式中：[E]⁺为[E]的广义逆矩阵，在MATLAB中可以通过pinv()函数进行求解。

有关广义逆矩阵理论的详细介绍可参见《基于广义逆矩阵理论求解线性最小二乘问题》。

则得到各个分段的拟合直线的独立系数后，则各个分段的拟合直线的所有系数可表示为：

基于上述方法，本文编制了相应的MATLAB计算程序，为验证计算方法的正确性，本文构造一个包含数值噪声的正弦函数：

并且设置分段线性直线的断点为π/2和3π/2，通过上述方法计算得到待拟合数据点和相应的分段线性拟合结果如图3所示。

 

图3 含数值噪声的正弦函数的分段线性拟合结果

从图中可以看到，基于本文提出的方法可以很好地对正弦函数进行分段线性拟合，并且相邻拟合直线在断点处能够保持连续性。事实上，只需对本文所述方法进行稍加改进，很容易即可将其扩展到分段样条拟合等更高次的拟合函数中。

但需要指出的是，本文提出的分段线性拟合方法需要提前指定断点的数量和位置，而分段断点的数量和布置情况同样对拟合效果有很大影响，本文由于已知待拟合的数据点近似呈正弦函数的分布，因此可以将极值作为分段断点。但在实际情况中，如何合理布置分段断点也是在进行分段线性拟合时可能面临的问题。例如，对于图4中的待拟合数据点，采用不同的断点布置均能获得较为满意的拟合效果，但显然其拟合误差并不相同。当断点数量较少时，我们可以通过穷举的方法来寻找各个断点的最佳位置；而如果断点数量较多，则可能需要用到一些优化算法，通过合理选择断点位置来实现对拟合数据的最佳逼近。

 

 

 

图4 采用不同断点布置的分段线性拟合结果

clear,clc
%***************************待拟合数据及参数输入***************************
Data_x=linspace(0,2*pi,101)';
Data_y=sin(Data_x)+(0.4*rand(101,1)-0.2);       %待拟合数据点
Data_num=size(Data_x,1);                        %待拟合数据点总数
Data_s=[pi/2;3*pi/2];                           %分段直线断点布置
Interval_num=size(Data_s,1)+1;                  %区段数量
Data_s=[Data_s;max(Data_x)];
%**************************************************************************

%********************确定待拟合数据点所在区段的索引值**********************
Interval_index=zeros(Data_num,1);
Start_index=1;
for i=1:Data_num
  for j=Start_index:Interval_num
    %区段索引为1的特殊情况
    if j==1
      if Data_x(i)<=Data_s(j)
        Interval_index(i)=j;
        break;
      end
      continue;
    end
    %区段索引不为1的一般情况
    if (Data_x(i)>Data_s(j-1)&&Data_x(i)<=Data_s(j))
      Interval_index(i)=j;
      Start_index=j;
      break;
    end
  end
end
%**************************************************************************

%***************************确定等价关系矩阵A******************************
%初始化矩阵A
Matrix_A=zeros(Data_num,2*Interval_num);
%计算矩阵A的分量
for i=1:Data_num                   %行循环
  Matrix_A(i,2*(Interval_index(i)-1)+1)=Data_x(i);
  Matrix_A(i,2*(Interval_index(i)-1)+2)=1;
end
%**************************************************************************

%***************************生成连续性条件矩阵*****************************
%初始化存放不同区段的连续性条件矩阵的元胞数组
Matrix_C=cell(Interval_num-1,1);
%循环创建不同区段的连续性条件矩阵
for j=1:Interval_num-1
  Matrix_C{j,1}=[-1/Data_s(j) 1 1/Data_s(j)];
end
%**************************************************************************

%**************************组装转换矩阵分矩阵******************************
%初始化存放不同区段的转换矩阵分矩阵的元胞数组
Matrix_T=cell(Interval_num-1,1);
%考虑i=1的特殊情况[T1]
Matrix_T{1}=zeros(2*Interval_num,2*Interval_num-1);
Matrix_T{1}(1,1:3)=Matrix_C{1};               %填充连续性条件矩阵[C]
Matrix_T{1}(2:end,:)=eye(2*Interval_num-1);   %填充单位矩阵
%考虑i>1的一般情况[Ti]
for j=2:Interval_num-1
  Matrix_T{j}=zeros(2*Interval_num-j+1,2*Interval_num-j);
  Matrix_T{j}(1:j-1,1:j-1)=eye(j-1);          %填充上侧单位矩阵
  Matrix_T{j}(j,j:j+2)=Matrix_C{j};           %填充连续性条件矩阵
  Matrix_T{j}(j+1:end,j:end)=eye(2*Interval_num-2*j+1);   %填充下侧单位矩阵
end
%**************************************************************************

%****************************计算分段拟合直线系数**************************
%计算转换矩阵[Pai_T]
Matrix_Pai_T=Matrix_T{1};
for j=2:Interval_num-1
  Matrix_Pai_T=Matrix_Pai_T*Matrix_T{j};     %通过连乘分矩阵获得转换矩阵
end
%计算独立变量系数计算矩阵[E]
Matrix_E=Matrix_A*Matrix_Pai_T;
%计算独立变量系数计算矩阵[E]的广义逆矩阵
Matrix_E_inv=pinv(Matrix_E);
%计算分段样条曲线的系数矩阵
Matrix_Coeff=Matrix_Pai_T*Matrix_E_inv*Data_y;
%**************************************************************************

%*****************************绘制分段拟合直线*****************************
%初始化用于绘制分段拟合直线的元胞数组
Fitting_x=cell(Interval_num,1);
Fitting_y=cell(Interval_num,1);
for i=1:Interval_num
  if i==1
    Fitting_x{i}=linspace(min(Data_x),Data_s(i),51)';
  else
    Fitting_x{i}=linspace(Data_s(i-1),Data_s(i),51)';
  end
end
%计算分段样条拟合直线
for i=1:Interval_num
  %基于分段样条拟合直线的系数进行计算
  Fitting_y{i}=Matrix_Coeff(2*(i-1)+1)*Fitting_x{i}+...
    Matrix_Coeff(2*(i-1)+2);
end
%绘制分段样条直线
figure;
scatter(Data_x,Data_y,70,"Filled","Black");hold on;
for i=1:Interval_num
  plot(Fitting_x{i},Fitting_y{i},"LineWidth",2,"Color","red");
end
xlabel("x");ylabel("y");title("分段线性拟合结果");hold off;grid on;
%**************************************************************************