# 最小二乘法 模型_利用最小二乘法实现曲线拟合

1. 简介

在工程问题的计算中，我们经常需要处理一些离散数据的拟合问题，而最小二乘法是处理曲线拟合问题的常用方法。目前，许多软件都提供有基于最小二乘法进行曲线拟合的功能，例如在Origin和Excel中均可直接利用离散数据进行曲线拟合。然而，这些软件只能处理一些简单函数的拟合问题，当需要拟合的函数较为复杂时，或者无法用简单的表达式来表述时，则往往无法直接进行拟合。为此，本文将对最小二乘法的基本原理做简单介绍，随后介绍如何利用Matlab的lsqnonlin函数处理复杂函数的拟合问题。

1.1 曲线拟合的最小二乘法原理

利用最小二乘法进行曲线拟合的本质为寻找某个近似函数φ(x)，使得该函数与离散点之间尽可能地逼近。若将偏差定义为近似函数的近似值φ(x_i)与离散点y_i*之间的差值：

显然，偏差的大小可以用来反映近似函数的好坏。若该近似函数可使得偏差的平方和达到最小，即：

则这种求近似函数的方法即为离散数据曲线拟合的最小二乘法，函数φ(x)为最小二乘拟合函数，一般可选为一些简单的函数，例如选为多项式函数，则也称为最小多项式拟合。假定该多项式可表示为如下形式：

其中m为多项式阶次，n为离散数据对的个数。

根据最小二乘法的原理，则该问题实际上为选取参数a_j(j=0,1,…,m)使得如下所示的目标函数达到最小值：

根据多元函数求极值的必要条件，可得如下所示的方程组：

将上式移项并展开可得到如下形式的线性方程组：

求解上述线性方程组即可得到拟合多项式的系数。

2. 利用Matlab处理曲线拟合问题

基于上述计算原理，Matlab提供了polyfit函数用于处理多项式曲线的拟合问题，对于一些较为复杂但仍可通过简单表达式进行表述的函数，也可以利用Matlab的拟合工具箱(Curve fitting)进行拟合。但在某些情况，当拟合函数非常复杂，以致于无法用简单表达式进行表述时(例如分段函数以及涉及到条件语句)，则无法使用拟合工具箱进行拟合。对于此类问题，可以使用Matlab优化工具箱中的lsqnonlin函数进行解决。

2.1 lsqnonlin函数

lsqnonlin函数用于求解以下述形式表示的非线性最小二乘法拟合问题：

在使用该函数进行最小二乘法拟合时，lsqnonlin函数并不需要用户提供min ||f(x)||(平方和)，而是需要用户提供自定义函数fun，用于计算矢量形式表示的f(x)：

lsqnonlin函数常用语法为：

x = lsqnonlin(fun,x0)x = lsqnonlin(fun,x0,lb,ub)

其中fun为用户自定义函数，x0为计算采用的初始值，lsqnonlin函数首先利用x₀通过自定义函数fun计算f_i(x)的取值并计算平方和，随后通过优化算法调整x的取值直至得到平方和的最小值。此外，lb和ub还可以用于定义x的取值范围，使得x满足lb≤x≤ub。

例如，对于节1.1中所述的多项式，根据最小二乘法的定义，则自定义函数f(x)应表示为：

注意此时f(x)中的x为以向量形式表示的多项式P(x)的系数：

在计算时，用户需要指定多项式系数的初始值，则lsqnonlin函数将利用最小二乘法计算多项式系数。

下面，本文将以笔者所在领域常用的NASGRO方程为例，介绍如何利用lsqnonlin函数处理此类复杂函数的曲线拟合问题。

2.2 NASGRO方程简介

在进行基于断裂力学的损伤容限分析时，应力强度因子和裂纹扩展速率模型是最为重要的输入。一般来说，应力强度因子可以通过经验公式或数值方法进行计算，而裂纹扩展速率模型则需要通过裂纹扩展速率试验获得的试验数据拟合得到。例如，大量的试验结果表明，在裂纹扩展的中速率区域，应力强度因子幅值ΔK和裂纹扩展速率da/dN满足良好的对数线性关系，可以通过Paris公式进行描述：

其中C和m为材料常数。

尽管Paris公式已经得到广泛的应用，但是Paris公式仅仅描述了裂纹在中速率区域的扩展行为，没有描述近门槛区域和接近断裂的高速率区域的扩展行为，也没有考虑应力比R和裂纹闭合效应对裂纹扩展速率的影响，因此给出的计算结果将过于保守。另一个常用的裂纹扩展速率模型为Newman提出的NASGRO模型，该模型基于Forman模型改进了裂纹扩展速率模型，同时比Paris和Walker模型更加全面，不仅考虑了应力强度因子门槛值和断裂韧度，还体现了应力比以及裂纹闭合效应对裂纹扩展速率da/dN的影响，如图2.1所示，其表达式如下：

其中R为应力比，ΔK为应力强度因子幅值，ΔK_th为应力强度因子幅值门槛值，K_max为最大应力强度因子，可表示为：

K_c为材料断裂韧度，C, m, p, q为材料常数，f为Newman裂纹张开函数，可表示为：

其中，

其中，α为塑性约束因子，σ_max为最大应力水平，σ_F为流动应力，计算时可以取为单轴屈服应力和单轴强度极限的平均值。

图2.1 NASGRO方程

NASGRO方程中的应力强度因子门槛值ΔKth可采用下面的经验公式进行估算：

其中A₀为裂纹张开函数中的多项式系数，ΔK₁是R=1时的应力强度因子门槛值，C_th是对于正应力(