基于样条函数的监测数据平滑算法与Matlab实现

杨海星，黎 浩，曹 净

(1.中国建筑西南勘察设计研究院有限公司，成都 610051；2.长春工程学院，长春 130021； 3.昆明理工大学，昆明 650000)

0 引言

由于岩土体性质的复杂多变性，以及各种计算模型的局限性，仅依靠理论分析和经验估计很难预测工程结构和土体在施工过程中的变化。为了保证工程安全顺利地进行，在施工过程中开展严密的监测已经成为工程建设必不可少的重要环节[1]。通过对监测数据的分析，可充分挖掘其蕴含的信息，实时掌握工程的动态，必要时可对原有方案进行调整以避免发生工程事故或降低工程造价。然而监测数据往往是借助测量仪器直接或间接获得的，这一过程中总会不可避免地产生测量误差，直接对该监测数据进行分析只可大概地了解工程的动态变化。若想从该监测数据中挖掘更多关于工程的信息，则需对监测数据进行更进一步的处理，如求导等。但因监测数据中含有误差，这些微小的误差在分析过程中将会被放大，甚至直接影响分析结果。为此，在利用监测数据进行进一步分析前需对其进行平滑处理，尽可能地将误差降低至最小。目前，数据平滑处理的方法有很多，如能量法[2-4]、小波分解法[5-6]、最小二乘法[7-8]、选点修改法[9]和节点弃除法[10]等，这些方法虽计算速度快，但很难保证曲线的曲率均匀变化，平滑处理效果差，且有时会顾此失彼，即平滑后的曲线过于平滑而偏离原始数据，逼近效果不好。

基于上述考虑，本文首先依据一定的光顺准则和逼近准则建立泛函，然后基于B—样条函数，构造奇次光顺样条函数，建立方程组求解泛函极小值，所求泛函极小值即为光顺样条函数。

1 问题的数学描述

(1)

上式称为光顺准则。在点xi，函数f(x)要逼近数组yi，(i=1，2，…，N)，逼近程度可用E(f)来估计：

(2)

J(f)=Iq(f)+ρEq(f)。

(3)

2 泛函极值的求解

定义：给定一个划分Δ：a=x0

1)s(x)在每个子区间[xi，xi+1](i=1，2，…，N-1)上是2q-1次代数多项式；

2)s(x)在每个子区间[a，x1)与(b，x1]上是q-1次代数多项式；

3)s(x)∈C(2q-2)[a，b]，1≤q≤N。

σΔ+ρ-1(-1)qd(σ)=y，

(4)

则称σ(x)为以ρ>0为权系数的2q-1次光顺样条函数，其中

σΔ=σ(x1)，σ(x2)，…，σ(xN)T，

d(σ)=d1(σ)，d2(σ)，…，dN(σ)T，

di(σ)=σ(2q-1)(xi+0)-σ(2q-1)(xi-0)，i=1，2，…，N，

y=y1，y2，…，yNT。

3 奇次光顺样条的计算

给定数据点(xi，yi)，i=1，2，…，N，并给定以xi为内节点的分划

Δ：a=x0

(5)

σΔ+ρ-1(-1)qd(σ)=y。

(6)

(7)

其中cj(j=1，2，…，N)为待定常数。又σ(x)还可以表示为

(8)

若把基底Bj(x)也写为(8)的形式

(9)

其中pj(x)是q-1次多项式，则

由此可得到

(10)

则方程组(6)可写为

将其写为矩阵的形式为

(B+ρ-1E)c=y。

(11)

其中

c=(c1，c2，…，cN)T，

y=(y1，y2，…，yN)T。

3.1 βij的确定

当j=1，2，…，q时

φ2q-1x；x1，…，xq+j=

(12)

当j=q+1，q+2，…，N-q时

(xj+q-xj-q)φ2q-1x；xj-q，…，xj+q=

(13)

当j=N-q+1，N-q+2，…，N时

φ2q-1xN-x；xN-xN，…，xN-xj-q=

(14)

其中

(15)

(16)

(17)

3.2 权因子ρ的确定

参数ρ可以调节逼近函数σρ(x)与数据之间的接近程度和逼近函数σρ(x)本身的“光滑”程度，所以它的选取至关重要。ρ选取过大会使σρ(x)过分依赖数据y，而须知y是有误差的。ρ选取过小又会产生一个基本不依赖于数据的过分光顺的样条，以致使变分为0。为了得到合适的ρ值，Reinsoh曾暗示[12]，如果方差σ2已知，那么可选取ρ使成立。

(18)

参数ρ的取值可通过牛顿迭代法获得，具体如下：首先将式(18)改写为

令式(18)左端为F(ρ)，将式(7)代入可得

对F(ρ)进行求导有：

其中

B(B+ρ-1E)-1ρ-2E(B+ρ-1E)-1y。

整理可得：

式中A=B+ρ-1E。

利用牛顿迭代法可求出权因子ρ，其迭代关系式如下：

(19)

4 光顺样条的算法及汇编程序

4.1 算法

对给定数据点，求光顺样条函数的步骤如下：

1)给定y=(y1，y2，…，yN)T、标准差σ以及初始值ρ；

2)计算bij=Bj(xi)，i，j=1，2，…，N，并组成矩阵B；

3)由式(15)～(17)求出βij，并求出eij=(-1)qβij，其中i，j=1，2，…，N，并组成矩阵E；

4)计算A=B+ρ-1E，形成矩阵A；

4.2 汇编程序

利用Matlab软件中的文本文件编辑器，创建了5个M文件(1个主文件，4个子文件)，通过主M文件中的Matlab指令完成对其他子M文件的调用，最终完成算法的Matlab实现。各M文件的功能见表1。

表1 各M文件的功能

各M文件的代码如下：

1)function [Y]=smoothdata(kk，x，y，s，deta)。% [Y]=smoothdata(kk，x，y，s，XX，deta)，各参数含义如下：kk为平滑样条次数，kk=2q-1，q为正整数；x为向量x=[x1，x2，x3，…，xN]，平滑数据点的横坐标值；y为向量y= [y1，y2，y3，…，yN]，平滑数据点的纵坐标值；s为给定的平滑程度控制值，人为给定，一般s=Nσ2；w为行向量w=[w1，w2，w3，..，wN]，wi表示第i个数据点xi的权重，若不赋值，则取默认值wi=1；Y为向量Y=[Y1，Y2，…，YN]，用于存储经平滑处理后的y值；N为数据点的个数，要求N≥kk+2；rou为光顺与逼近得权函数，用牛顿迭代法求解；c为列向量c=[c1，c2，…，cN]，为各基函数的系数；

%程序部分

if nargin<5；%给定各个数据点的权重，如果没有赋值，则取默认值为1

deta=ones(1，length(x))；

end

detamax=max(deta)；

w=detamax./deta；

rou=1*10^(-10)；

F=s+1；

Stepp1=0；

[B，E]=BE(kk，x，w)；

y=y'；

while F>s

A=B+rou^(-1)*E；

c=Ay；

m=B*c-y；

F=m'*m；

if F>s

invA=inv(A)；

dFdr=2*rou^(-2)*m'*B*invA*E*c；

rou=rou-(F-s)/(dFdr)；

else

break

end

end

%计算平滑数据

BB=RB(kk，x，x)；

Y=BB*c；

Y=Y'；

2)function [B，E]=BE(kk，x，w)。%[B，E]=BE(kk，x，w)，各参数含义如下：B为矩阵B，用于存储各基函数在各节点的值；E为矩阵E，用于存储各基函数在各节点的(2q-1)阶导数的跳跃量；kk为奇次自然光顺样条函数的次数；x为插值点列(x1，x2，x3，…，xN)；N为插值点列的个数，N≥kk+2；q，q=(kk+1)/2

%)程序

%%计算各基函数在各节点的值

q=(kk+1)/2；

N=length(x)；

for i=1∶1∶N

for j=1∶1∶q

xx=x(1，1∶q+j)；

B(i，j)=bjxi(kk，xx，x(i))；

if i<=q+j&&i>=1

E(i，j)=w(i)^(-1)*

beta1(kk，xx，x(i))；

end

end

for j=q+1∶1∶N-q

xx=x(1，j-q∶1∶j+q)；

B(i，j)=(x(j+q)-x(j-q))*bjxi(kk，xx，x(i))；

if i<=j+q&&i>=j-q

E(i，j)=w(i)^(-1)*

(x(j+q)-x(j-q))*beta1(kk，xx，x(i))；

end

end

for j=N-q+1∶1∶N

xx=x(N)-x(1，N∶-1∶j-q)；

B(i，j)=bjxi(kk，xx，x(N)-x(i))；

if i<=N&&i>=j-q

E(i，j)=w(i)^(-1)*beta1(kk，

xx，x(N)-x(i))；

end

end

end

E=(-1)^q*E；

3) function [B]=RB(kk，x，X)。%[B]=RB(kk，x，X)；各参数含义如下：B用于存储各基函数在X的值；kk为奇次自然光顺样条函数的次数；x为插值点列(x1，x2，x3，…，xN)；N为插值点列的个数，N≥kk+2；q=(kk+1)/2；X为计算点的坐标

%程序

q=(kk+1)/2；

N=length(x)；

n=length(X)；

for i=1∶1∶n

for j=1∶1∶q

xx=x(1，1∶q+j)；

B(i，j)=bjxi(kk，xx，X(i))；

end

for j=q+1∶1∶N-q

xx=x(1，j-q∶1∶j+q)；

B(i，j)=(x(j+q)-x(j-q))*bjxi(kk，xx，X(i))；

end

for j=N-q+1∶1∶N

xx=x(N)-x(1，N∶-1∶j-q)；

B(i，j)=bjxi(kk，xx，x(N)-X(i))；

end

end

4) function [value]=bjxi(kk，xx，X)。%[value]=bjxi(kk，xx，X)，各参数含义如下：value为计算第j个基第在X处的值；kk为奇次自然光顺样条函数的次数；x为插值点列(x1，x2，x3，…，xN)，X为计算点处的坐标。

%程序

n=length(xx)；

value=zeros(1，length(X))；

for i=1∶1∶length(X)

XX=X(i)；

for k=1∶1∶n

if (xx(k)-XX)>0

m=(xx(k)-XX)^(kk)；

else

m=0；

end

for l=1∶1∶n

if l～=k

m=m/(xx(k)-xx(l))；

end

end

value(i)=value(i)+m；

end

end

5) function[beta]=beta1(kk，xx，X)。%[beta]=beta1(kk，xx，X)，各参数含义如下：beta为计算基函数的(2q-1)阶导数在X处的跳跃量；kk为奇次自然光顺样条函数的次数；x为插值点列(x1，x2，x3，…，xN)；X为计算点处的坐标。

%程序部分

n=length(xx)；

m=factorial(kk)；

for j=1∶1∶n

if X～=xx(j)

m=m/(X-xx(j))；

end

beta=m；

end

5 实例验证

已知函数为y=10sin(πx/50)，并在该已知函数上加上一组随机数，该随机数服从均匀分布U～(0，1)，相当于在原函数上加最大函数值的10%的随机扰动。然后采用本文所提数据平滑法(五次样条函数平滑法)对该含噪音数据进行平滑处理，将平滑结果与原数据进行比较。原始数据值与噪音值如图1，含噪音数据与平滑后的数据如图2，平滑数据的一阶导数(曲率)如图3。通过图2～3可发现，经过平滑处理后的数据能很大程度上消除噪音的影响，可很好地逼近原始数据，所得曲线具有较好的光滑性，其曲率变化均匀。通过该算例证明了本文所提数据平滑算法的可行性。

图1 原始数据值与噪音值

图2 含噪音数据与平滑后的数据

6 结语

1)工程监测数据中不可避免地含有误差，在进行数据分析前需对原始数据进行平滑处理，将误差对数据分析结果的影响降至最低。本文首先依据一定的光顺准则和逼近准则建立泛函，将数据平滑问题转化为泛函求极值的问题，然后基于B—样条函数，构造奇次光顺样条函数，建立方程组求解泛函极小值，所求泛函极小值即为光顺样条函数，该函数既有一定的光顺性，又具有较好的逼近性能。

图3 平滑曲线曲率

2)依据光顺样条函数的求解过程，给出数据平滑算法，并利用Matlab汇编语言将算法程序化，最后通过一实例验证了该数据平滑算法的可行性。