楊海星，黎 浩，曹 凈

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

0 引言

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

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

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

(1)

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

(2)

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

(3)

2 泛函極值的求解

定義：給定一個(gè)劃分Δ：a=x0

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

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

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 奇次光順樣條的計(jì)算

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

Δ：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次多項(xiàng)式，則

由此可得到

(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時(shí)

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

(12)

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

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

(13)

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

φ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)重要。ρ選取過(guò)大會(huì)使σρ(x)過(guò)分依賴數(shù)據(jù)y，而須知y是有誤差的。ρ選取過(guò)小又會(huì)產(chǎn)生一個(gè)基本不依賴于數(shù)據(jù)的過(guò)分光順的樣條，以致使變分為0。為了得到合適的ρ值，Reinsoh曾暗示[12]，如果方差σ2已知，那么可選取ρ使成立。

(18)

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

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

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

其中

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

整理可得：

式中A=B+ρ-1E。

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

(19)

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

4.1 算法

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

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

2)計(jì)算bij=Bj(xi)，i，j=1，2，…，N，并組成矩陣B；

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

4)計(jì)算A=B+ρ-1E，形成矩陣A；

4.2 匯編程序

利用Matlab軟件中的文本文件編輯器，創(chuàng)建了5個(gè)M文件(1個(gè)主文件，4個(gè)子文件)，通過(guò)主M文件中的Matlab指令完成對(duì)其他子M文件的調(diào)用，最終完成算法的Matlab實(shí)現(xiàn)。各M文件的功能見(jiàn)表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ù)點(diǎn)的橫坐標(biāo)值；y為向量y= [y1，y2，y3，…，yN]，平滑數(shù)據(jù)點(diǎn)的縱坐標(biāo)值；s為給定的平滑程度控制值，人為給定，一般s=Nσ2；w為行向量w=[w1，w2，w3，..，wN]，wi表示第i個(gè)數(shù)據(jù)點(diǎn)xi的權(quán)重，若不賦值，則取默認(rèn)值wi=1；Y為向量Y=[Y1，Y2，…，YN]，用于存儲(chǔ)經(jīng)平滑處理后的y值；N為數(shù)據(jù)點(diǎn)的個(gè)數(shù)，要求N≥kk+2；rou為光順與逼近得權(quán)函數(shù)，用牛頓迭代法求解；c為列向量c=[c1，c2，…，cN]，為各基函數(shù)的系數(shù)；

%程序部分

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

%計(jì)算平滑數(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，用于存儲(chǔ)各基函數(shù)在各節(jié)點(diǎn)的值；E為矩陣E，用于存儲(chǔ)各基函數(shù)在各節(jié)點(diǎn)的(2q-1)階導(dǎo)數(shù)的跳躍量；kk為奇次自然光順樣條函數(shù)的次數(shù)；x為插值點(diǎn)列(x1，x2，x3，…，xN)；N為插值點(diǎn)列的個(gè)數(shù)，N≥kk+2；q，q=(kk+1)/2

%)程序

%%計(jì)算各基函數(shù)在各節(jié)點(diǎn)的值

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用于存儲(chǔ)各基函數(shù)在X的值；kk為奇次自然光順樣條函數(shù)的次數(shù)；x為插值點(diǎn)列(x1，x2，x3，…，xN)；N為插值點(diǎn)列的個(gè)數(shù)，N≥kk+2；q=(kk+1)/2；X為計(jì)算點(diǎn)的坐標(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ì)算第j個(gè)基第在X處的值；kk為奇次自然光順樣條函數(shù)的次數(shù)；x為插值點(diǎn)列(x1，x2，x3，…，xN)，X為計(jì)算點(diǎn)處的坐標(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為計(jì)算基函數(shù)的(2q-1)階導(dǎo)數(shù)在X處的跳躍量；kk為奇次自然光順樣條函數(shù)的次數(shù)；x為插值點(diǎn)列(x1，x2，x3，…，xN)；X為計(jì)算點(diǎn)處的坐標(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àn)證

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

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

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

6 結(jié)語(yǔ)

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

圖3 平滑曲線曲率

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