用matlab进行线性拟合（基于matlab的信号描述及其运算实验原理）

在声发射中可能有这么几个应用场景。

1、在做加载试验的时候，获得一组声发射事件数随载荷的变化，后期当我们只测得声发射信号事件数的时候，能否评估当前的载荷是多少，或者推断当前试验所处的状态？

2、对一个材料做疲劳试验的时候，随着疲劳的加深，声发射事件率也逐渐增加，当我们采集完数据之后能否通过声发射事件率判断材料当前所处的疲劳状态？

3、一个材料可能有多个模态的信号源，每个信号源产生的信号可能有明显的区别，如下图，其中一类数据点随着信号幅度的增加能量快速的增加，另一类信号幅度增加的同时能量增加得就没那么快了。能否通过信号的幅度和能量推测信号产生的模式呢？

以上分析，就可以使用本节介绍的数据拟合功能了。

数据拟合又称曲线拟合，俗称拉曲线，是一种把现有数据透过数学方法来代入一条数式的表示方式。科学和工程问题可以通过诸如采样、实验等方法获得若干离散的数据，根据这些数据，我们往往希望得到一个连续的函数（也就是曲线）或者更加密集的离散方程与已知数据相吻合，这过程就叫做拟合(fitting)。

其实数据拟合就是通过一系列的离散点找出一条曲线比如Y=a*X+c；后期当我们测到一个X值的时候，就可以推测对应的Y值。（曲线可能有不同类型，数据也可以是三维的甚至更多维度的，本处主要以二维为例）。

二、数据拟合的基本算法：最小二乘法

2.1、引言

我们以前做实验，为了提高数据的准确性，经常用的一个方法叫做多次测量取平均值。比如要测量一个东西的长度的时候，我们可能会测量7次，然后在7次的结果里去掉最大和最小值，剩下的数据取平均值。

假设测得的数据如下|：

我们先会将序号2和7的数据排除，然后将其他的数据求平均值。

xr=(10.2+9.8+10+9.8+9.9)/5=9.94

利用matlab将其绘制出来，得到如下图形

clear allclose allx=1:5;y1=[10.2,9.8,10,9.8,9.9];plot(y1,'r*');set(gca,'YLim',[9,11]);

2.2、最小二乘法

我们假设我们要测量的东西的实际长度是X，我们将其画在上图中。

我们可以得到每个测量点和实际长度之间的距离，也可以理解为测量值和真实值之间的误差。因为取绝对值计算起来麻烦，所以用平方来代替误差，5次测量的总的误差可以表示为：

我们假设虽然我们每次测量都会有误差，但是离真实值越近我们得到的可能会越大，比如我们可能会测量出结果是9.9,9.8,9.7，但是我们正常情况下是不会测量出结果是1的情况的。换言之，当Y值是真实值的时候，以上公式计算的结果是最小的。

以上公式是一个一元二次方程，且开口向上，根据微积分的原理我们可以知道当其导数是0的时候，公式取得最小的值。

对其求导：

进而得出

可见对于这个案例，最小二乘法的结果就是算术平均值。

2.3、一次方程

以上是最简单的方程，但是已经说明了最小二乘法的原理。

我们假设声发射断铅实验中信号的幅度和断铅距离传感器的距离做一个案例，假设数据如下：

绘制其曲线如下：

可以看出线性关系变成了f(x)=a*x+b；（上一个案例中是f(x)=a*x）

这个时候每个点和实际点的误差是f(xi)-yi，最小误差和可以表示为：

函数有两个参数a,b，分别对ab求导，两个导数都为0的时候对应的是误差最小的f(x0);

以上是二元一次方程组，解方程即可得到两个参数ab，根据ab即可以画出拟合曲线。

2.4、二次方程

当我们的拟合曲线是二次方程的时候，f(x)=a*x^2+b*x+c。

每个点距离真实点的距离是f(xi)-yi。

所有的的误差和为

根据微积分原理，分别对abc求导，可以得到3元一次方程组，方程组都为0时候得到的值是最小值，以下是a=1,b=0.3,c=2时候的数据和拟合曲线。

clear allclose allx=1:5;y1=x.^2+0.3*x+2;y1=y1+randn(1,5);plot(y1,'r*');set(gca,'XLim',[0,6]);set(gca,'YLim',[0,30]);

2.5、指数等其他方程

除了以上方程还可以用其他方程取拟合数据，比如指数方程f(x)=a*e^b，计算的方式也是如上，分别对未知参数ab求导，导数为0时候误差和最小。然后根据导数方程组求二元方程。

也可以根据实际写一个自定义的方程，比如f(x)=a*exp(b*x)+c*x^2+d*sin(e*x);这个就是包含了指数、次方、三角波的拟合函数。

三、matlab图表中自带的拟合功能介绍

3.1、使用代码生成数据，绘图结果如下所示

x=1:100;y1=x+randn(1,100)*5;y2=0.3*x+randn(1,100)*5;y3=x.*x/100;plot(x,y1,'r*',x,y2,'b+',x,y3,'g.');ylim([-120,120]);

3.2、图形界面选择“工具”->“基本拟合”，弹出设置对话框。

3.3、设置模块分为几个模块。分别是选择数据、选择拟合曲线样式、显示方程、绘制残差图和更多。

以下分别介绍

1、选择数据：我们的数据是两个数据，分别是y1和y2。这个数据可以由聚类分析等方法区分开。

2、选择拟合曲线样式：

如果选择样条插值或者保形插值，会将每个点连接。

如果选择直线，matlab会用一条直线拟合数据。

如果选择二次或者多次方程，会用二次方程式进行拟合

3、勾选显示方程，界面上会出现拟合曲线的方程式

4、勾选绘制残差，会在图形上显示每点和拟合曲线的差

5、点击更多按钮，右侧将出现具体的数值结果

如果只需要简单的分析或者大概寻找规律，基本就以上功能可以满足。

四、matlab中cftool工具箱的使用

在声发射数据中，声发射数据也可能不是标准的次方规律，可能是指数增长，可能是周期变化，也可能是多种变化的组合，那么，下面将介绍这种情况下的拟合。

本节继续介绍matlab中非常实用的数据拟合工具箱cftool。简单实用而强大的函数拟合工具箱，提供丰富的拟合算法，能实现多种线性、非线性的函数拟合。

相较于图表自带的拟合功能，cftool提供了更多的拟合函数选择、甚至可以自己编辑拟合的函数，可以应对更复杂的声发射信号分析。

4.1、生成几个测试用的数据

x=1:1:100;rd=randn(size(x))*3;y1=5*(x+rd)+3;%带随机的线性数据y2=3*x.*(x+rd)+0.5*x+7;%带随机的二次方数据y3=exp(0.08*x+rd*0.05);%指数分布的数据y4=sin(0.3*x)*5+rd*0.1;%三角函数分布的数据y5=sin(0.3*x)+0.1*x;%三角函数和线性叠加的数据

4.2、打开cftool工具箱

在命令行输入“cftool”，或者在菜单中点击“curve fitting”菜单，弹出工具箱操作界面。

4.3、选择数据源

首先给拟合一个命名“fity1”，Xdata选择x，Ydata选择y1；我们主要操作二维的，所以Zdata选择空，如果要做三维拟合可以选择一个Z数据，weights也选择空，暂时每个点的权重都是一样的，我们后面再介绍不同权重的点的拟合。

4.4、拟合数据

工具箱默认的是线性拟合，所以此时图中出现了数据点和拟合的曲线，如果没有自动出现拟合结果，可以手动点击界面上的‘FIT’按钮进行数据拟合。

在拟合曲线的左侧是拟合的结果，分别是拟合模型：线性y=p1*x+p2；拟合曲线的参数：p1p2的值；后面的拟合曲线对数据的解释能力，也就是拟合的准确性评估：

其中R-square越接近1表示拟合越完美。

4.5、选择其他拟合函数

点击拟合函数下拉框，会出现如下拟合函数，可以选择一个，之后会出现拟合的函数的示例方程，如果选择了自动拟合，就会即时出现拟合曲线和拟合结果。

4.6、拟合比较

我们选择y2作为数据，y2是2次方数据加随机数形成，选择线性拟合结果如下1，选择二次函数拟合结果如下2；可以通过拟合结果确定选择拟合函数的适用性。

4.7、其他拟合曲线

同理y3的最佳拟合函数是函数或者，其求出的值也和我们生成数据的值比较接近。

y4的最佳拟合函数是指数函数，其求出的值也和我们生成数据的值比较接近

4.8、自定义函数拟合

y5的数据有有一个线性和正弦信号叠加，所以我们用自定义函数。通过观察数据，可以写出拟合函数是y5=a*sin(b*x)+c*x+d，写入输入框中。采用默认的参数，数据拟合可能不准确，再次观察数据，可以获得正弦周期性变化的周期大概范围，点击fitopetion，给参数b一个范围，获得拟合的结果如下所示。

通过以上的数据就可获得数据的一个变化近似函数，将声发射的过程进行量化。

五、拟合的优化：分段拟合与局部加权拟合

上面介绍的主要都是基于基本的最小二乘法的拟合操作，本节是响应的优化。

5.1分段拟合

在有些实验中，实验的过程是可以按照时间段区分的。比如在负荷加载中，有加压——恒压——加压——恒压——卸压的过程，在数据分析中，可以分析一次加恒压过程中声发射的规律。如下，随着加压，声发射事件快速增加。恒压阶段增速减慢后基本保持稳定。

在数据拟合的时候，可能没有非常合适的函数能够表述上面的这条曲线，那么就可以使用分段拟合的方法。

分段拟合就是将数据分成多个连续的区间段，对每个区间段进行拟合。比如下面这个图，虽然也可以拟合成一条曲线，但是拟合的损失值会比分段拟合大。另外在实验室条件下的实验很多都是加有外部作用力的，拐点会比较明显，所以根据实验的过程进行分段拟合是个不错的选择。

最后可以拟合出如下的方程组

f(x)=x,when(x<10)

f(x)=3*x+5,when(x>10)

5.2局部加权拟合

局部加权回归使一种非参数方法(Non-parametric)。在每次预测新样本时会重新训练临近的数据得到新参数值。

局部加权回归优点：可以有效避免欠拟合，减小了较远数据的干扰，仅与较近的数据有关。缺点是计算量比较大。

如上图，我们要预测x点对应的y值，根据经验可以知道，距离x比较近的点对x的预测影响比较大，而较远的点影响比较小，如上虚线内的影响较大而虚线外的影响比较小。所以在对x进行预测的时候，就根据x的范围重新拟合一曲线，如上斜线，将x带入斜线方程求得对应的y值。

5.3在cftool工具中加权回归的功能

也有加权回归的功能，通常加权值和自变量取同样的数据。

我们在上节的y1数据上，将y1y10y20y30都改成偏离很大的值，如下图中图上的4个黑点。

我们将权重也选为x的时候，拟合相似度明显提升。

六、数据拟合在声发射数据分析中的实际使用

6.1《输气管道阀门内漏声发射多参数检测实验研究》

针对输气管道阀门内漏喷流噪声声场特征复杂特点，应用声发射检测系统进行阀门内漏过程的检测实验研究，提取阀门不同内漏率下声发射信号波形频谱特征和时域、频率参数特征；应用最小二乘法拟合函数曲线，获得声发射信号特征参数与内漏率对应关系。实验结果表明，声发射信号特征参数（均方根、标准差、均值、频域峰值）与内漏率之间都存在明显对应关系。