基于樣條函數(shù)的監(jiān)測數(shù)據(jù)平滑算法與Matlab實現(xiàn)

楊海星，黎 浩，曹 凈

(1.中國建筑西南勘察設(shè)計研究院有限公司，成都 610051；2.長春工程學(xué)院，長春 130021； 3.昆明理工大學(xué)，昆明 650000)

0 引言

由于巖土體性質(zhì)的復(fù)雜多變性，以及各種計算模型的局限性，僅依靠理論分析和經(jīng)驗估計很難預(yù)測工程結(jié)構(gòu)和土體在施工過程中的變化。為了保證工程安全順利地進(jìn)行，在施工過程中開展嚴(yán)密的監(jiān)測已經(jīng)成為工程建設(shè)必不可少的重要環(huán)節(jié)[1]。通過對監(jiān)測數(shù)據(jù)的分析，可充分挖掘其蘊含的信息，實時掌握工程的動態(tài)，必要時可對原有方案進(jìn)行調(diào)整以避免發(fā)生工程事故或降低工程造價。然而監(jiān)測數(shù)據(jù)往往是借助測量儀器直接或間接獲得的，這一過程中總會不可避免地產(chǎn)生測量誤差，直接對該監(jiān)測數(shù)據(jù)進(jìn)行分析只可大概地了解工程的動態(tài)變化。若想從該監(jiān)測數(shù)據(jù)中挖掘更多關(guān)于工程的信息，則需對監(jiān)測數(shù)據(jù)進(jìn)行更進(jìn)一步的處理，如求導(dǎo)等。但因監(jiān)測數(shù)據(jù)中含有誤差，這些微小的誤差在分析過程中將會被放大，甚至直接影響分析結(jié)果。為此，在利用監(jiān)測數(shù)據(jù)進(jìn)行進(jìn)一步分析前需對其進(jìn)行平滑處理，盡可能地將誤差降低至最小。目前，數(shù)據(jù)平滑處理的方法有很多，如能量法[2-4]、小波分解法[5-6]、最小二乘法[7-8]、選點修改法[9]和節(jié)點棄除法[10]等，這些方法雖計算速度快，但很難保證曲線的曲率均勻變化，平滑處理效果差，且有時會顧此失彼，即平滑后的曲線過于平滑而偏離原始數(shù)據(jù)，逼近效果不好。

基于上述考慮，本文首先依據(jù)一定的光順準(zhǔn)則和逼近準(zhǔn)則建立泛函，然后基于B—樣條函數(shù)，構(gòu)造奇次光順樣條函數(shù)，建立方程組求解泛函極小值，所求泛函極小值即為光順樣條函數(shù)。

1 問題的數(shù)學(xué)描述

(1)

上式稱為光順準(zhǔn)則。在點xi，函數(shù)f(x)要逼近數(shù)組yi，(i=1，2，…，N)，逼近程度可用E(f)來估計：

(2)

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

(3)

2 泛函極值的求解

定義：給定一個劃分Δ：a=x0

1)s(x)在每個子區(qū)間[xi，xi+1](i=1，2，…，N-1)上是2q-1次代數(shù)多項式；

2)s(x)在每個子區(qū)間[a，x1)與(b，x1]上是q-1次代數(shù)多項式；

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

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

(4)

則稱σ(x)為以ρ>0為權(quán)系數(shù)的2q-1次光順樣條函數(shù)，其中

σΔ=σ(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 奇次光順樣條的計算

給定數(shù)據(jù)點(xi，yi)，i=1，2，…，N，并給定以xi為內(nèi)節(jié)點的分劃

Δ：a=x0

(5)

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

(6)

(7)

其中cj(j=1，2，…，N)為待定常數(shù)。又σ(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的確定

當(dāng)j=1，2，…，q時

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

(12)

當(dāng)j=q+1，q+2，…，N-q時

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

(13)

當(dāng)j=N-q+1，N-q+2，…，N時

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

(14)

其中

(15)

(16)

(17)

3.2 權(quán)因子ρ的確定

參數(shù)ρ可以調(diào)節(jié)逼近函數(shù)σρ(x)與數(shù)據(jù)之間的接近程度和逼近函數(shù)σρ(x)本身的“光滑”程度，所以它的選取至關(guān)重要。ρ選取過大會使σρ(x)過分依賴數(shù)據(jù)y，而須知y是有誤差的。ρ選取過小又會產(chǎn)生一個基本不依賴于數(shù)據(jù)的過分光順的樣條，以致使變分為0。為了得到合適的ρ值，Reinsoh曾暗示[12]，如果方差σ2已知，那么可選取ρ使成立。

(18)

參數(shù)ρ的取值可通過牛頓迭代法獲得，具體如下：首先將式(18)改寫為

令式(18)左端為F(ρ)，將式(7)代入可得

對F(ρ)進(jìn)行求導(dǎo)有：

其中

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

整理可得：

式中A=B+ρ-1E。

利用牛頓迭代法可求出權(quán)因子ρ，其迭代關(guān)系式如下：

(19)

4 光順樣條的算法及匯編程序

4.1 算法

對給定數(shù)據(jù)點，求光順樣條函數(shù)的步驟如下：

1)給定y=(y1，y2，…，yN)T、標(biāo)準(zhǔn)差σ以及初始值ρ；

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軟件中的文本文件編輯器，創(chuàng)建了5個M文件(1個主文件，4個子文件)，通過主M文件中的Matlab指令完成對其他子M文件的調(diào)用，最終完成算法的Matlab實現(xiàn)。各M文件的功能見表1。

表1 各M文件的功能

各M文件的代碼如下：

1)function [Y]=smoothdata(kk，x，y，s，deta)。% [Y]=smoothdata(kk，x，y，s，XX，deta)，各參數(shù)含義如下：kk為平滑樣條次數(shù)，kk=2q-1，q為正整數(shù)；x為向量x=[x1，x2，x3，…，xN]，平滑數(shù)據(jù)點的橫坐標(biāo)值；y為向量y= [y1，y2，y3，…，yN]，平滑數(shù)據(jù)點的縱坐標(biāo)值；s為給定的平滑程度控制值，人為給定，一般s=Nσ2；w為行向量w=[w1，w2，w3，..，wN]，wi表示第i個數(shù)據(jù)點xi的權(quán)重，若不賦值，則取默認(rèn)值wi=1；Y為向量Y=[Y1，Y2，…，YN]，用于存儲經(jīng)平滑處理后的y值；N為數(shù)據(jù)點的個數(shù)，要求N≥kk+2；rou為光順與逼近得權(quán)函數(shù)，用牛頓迭代法求解；c為列向量c=[c1，c2，…，cN]，為各基函數(shù)的系數(shù)；

%程序部分

if nargin<5；%給定各個數(shù)據(jù)點的權(quán)重，如果沒有賦值，則取默認(rèn)值為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

%計算平滑數(shù)據(jù)

BB=RB(kk，x，x)；

Y=BB*c；

Y=Y'；

2)function [B，E]=BE(kk，x，w)。%[B，E]=BE(kk，x，w)，各參數(shù)含義如下：B為矩陣B，用于存儲各基函數(shù)在各節(jié)點的值；E為矩陣E，用于存儲各基函數(shù)在各節(jié)點的(2q-1)階導(dǎo)數(shù)的跳躍量；kk為奇次自然光順樣條函數(shù)的次數(shù)；x為插值點列(x1，x2，x3，…，xN)；N為插值點列的個數(shù)，N≥kk+2；q，q=(kk+1)/2

%)程序

%%計算各基函數(shù)在各節(jié)點的值

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)；各參數(shù)含義如下：B用于存儲各基函數(shù)在X的值；kk為奇次自然光順樣條函數(shù)的次數(shù)；x為插值點列(x1，x2，x3，…，xN)；N為插值點列的個數(shù)，N≥kk+2；q=(kk+1)/2；X為計算點的坐標(biāo)

%程序

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)，各參數(shù)含義如下：value為計算第j個基第在X處的值；kk為奇次自然光順樣條函數(shù)的次數(shù)；x為插值點列(x1，x2，x3，…，xN)，X為計算點處的坐標(biāo)。

%程序

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)，各參數(shù)含義如下：beta為計算基函數(shù)的(2q-1)階導(dǎo)數(shù)在X處的跳躍量；kk為奇次自然光順樣條函數(shù)的次數(shù)；x為插值點列(x1，x2，x3，…，xN)；X為計算點處的坐標(biāo)。

%程序部分

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 實例驗證

已知函數(shù)為y=10sin(πx/50)，并在該已知函數(shù)上加上一組隨機(jī)數(shù)，該隨機(jī)數(shù)服從均勻分布U～(0，1)，相當(dāng)于在原函數(shù)上加最大函數(shù)值的10%的隨機(jī)擾動。然后采用本文所提數(shù)據(jù)平滑法(五次樣條函數(shù)平滑法)對該含噪音數(shù)據(jù)進(jìn)行平滑處理，將平滑結(jié)果與原數(shù)據(jù)進(jìn)行比較。原始數(shù)據(jù)值與噪音值如圖1，含噪音數(shù)據(jù)與平滑后的數(shù)據(jù)如圖2，平滑數(shù)據(jù)的一階導(dǎo)數(shù)(曲率)如圖3。通過圖2～3可發(fā)現(xiàn)，經(jīng)過平滑處理后的數(shù)據(jù)能很大程度上消除噪音的影響，可很好地逼近原始數(shù)據(jù)，所得曲線具有較好的光滑性，其曲率變化均勻。通過該算例證明了本文所提數(shù)據(jù)平滑算法的可行性。

圖1 原始數(shù)據(jù)值與噪音值

圖2 含噪音數(shù)據(jù)與平滑后的數(shù)據(jù)

6 結(jié)語

1)工程監(jiān)測數(shù)據(jù)中不可避免地含有誤差，在進(jìn)行數(shù)據(jù)分析前需對原始數(shù)據(jù)進(jìn)行平滑處理，將誤差對數(shù)據(jù)分析結(jié)果的影響降至最低。本文首先依據(jù)一定的光順準(zhǔn)則和逼近準(zhǔn)則建立泛函，將數(shù)據(jù)平滑問題轉(zhuǎn)化為泛函求極值的問題，然后基于B—樣條函數(shù)，構(gòu)造奇次光順樣條函數(shù)，建立方程組求解泛函極小值，所求泛函極小值即為光順樣條函數(shù)，該函數(shù)既有一定的光順性，又具有較好的逼近性能。

圖3 平滑曲線曲率

2)依據(jù)光順樣條函數(shù)的求解過程，給出數(shù)據(jù)平滑算法，并利用Matlab匯編語言將算法程序化，最后通過一實例驗證了該數(shù)據(jù)平滑算法的可行性。