楊秋偉 陳 華 李翠紅

(1.紹興文理學(xué)院 土木工程學(xué)院，浙江 紹興 312000； 2.寧波工程學(xué)院 建筑與交通工程學(xué)院，浙江 寧波 315211；3.浙江省土木工程工業(yè)化建造工程技術(shù)研究中心(寧波工程學(xué)院)，浙江 寧波 315211)

0 引言

有限單元法已廣泛用于土木、機(jī)械、航空航天、汽車(chē)、船舶等眾多工程設(shè)計(jì)領(lǐng)域中，已經(jīng)成為解決復(fù)雜工程分析計(jì)算問(wèn)題的有效途徑.利用實(shí)際工程結(jié)構(gòu)的有限元模型，可開(kāi)展力學(xué)計(jì)算、響應(yīng)預(yù)測(cè)、優(yōu)化設(shè)計(jì)和損傷識(shí)別等研究.因此，一個(gè)準(zhǔn)確的有限元模型是上述工作能夠成功開(kāi)展的基礎(chǔ).然而，由于工程結(jié)構(gòu)的復(fù)雜性，用有限元軟件所建立的結(jié)構(gòu)模型與真實(shí)結(jié)構(gòu)總是存在著一定的差異，具體體現(xiàn)在：由有限元模型計(jì)算所得的結(jié)構(gòu)響應(yīng)參數(shù)(如位移、頻率和振型)總是與儀器測(cè)試所得的結(jié)構(gòu)響應(yīng)參數(shù)存在著差異.當(dāng)這種差異較大時(shí)，則必須對(duì)結(jié)構(gòu)的有限元模型進(jìn)行修正，即通過(guò)修改模型中部分單元的物理參數(shù)，來(lái)使得計(jì)算所得的響應(yīng)參數(shù)與測(cè)試所得的響應(yīng)參數(shù)盡可能接近，修正后的有限元模型方可用于力學(xué)分析與計(jì)算中，這種方法稱為模型修正法[1-7].

數(shù)據(jù)誤差是困擾模型修正法成功實(shí)施的一個(gè)關(guān)鍵問(wèn)題.一方面，由于測(cè)試設(shè)備和條件的限制，以及環(huán)境因素(如溫度、濕度、偶然荷載)的干擾，測(cè)試所得的結(jié)構(gòu)響應(yīng)參數(shù)必然存在著或大或小的誤差；另一方面，結(jié)構(gòu)建模時(shí)所采用的物理參數(shù)如彈性模量等的理論值和實(shí)際值也總是存在一定的偏差，當(dāng)采用這些含有誤差的數(shù)據(jù)進(jìn)行模型修正計(jì)算時(shí)，往往會(huì)出現(xiàn)病態(tài)最小二乘問(wèn)題[8-10].為了解決病態(tài)最小二乘問(wèn)題，許多學(xué)者開(kāi)展了嶺估計(jì)和廣義嶺估計(jì)方面的研究[11-23]，這類方法的主要難點(diǎn)在于嶺參數(shù)的選擇，常見(jiàn)的L曲線法或嶺跡法均需要很大的計(jì)算量且所得結(jié)果精度有限.有鑒于此，本文提出一種新的廣義嶺估計(jì)方法，所提方法充分考慮系數(shù)矩陣對(duì)角元素之間的差異，提出一個(gè)底數(shù)可調(diào)節(jié)的正則系數(shù)計(jì)算公式，用于獲得廣義嶺估計(jì)中的正則對(duì)角矩陣，可以盡量減少?gòu)V義嶺估計(jì)中由于引入正則矩陣而帶來(lái)的額外誤差，更好地兼顧解的穩(wěn)定性和準(zhǔn)確性，計(jì)算快捷且精度高.以一個(gè)病態(tài)方程組和一個(gè)桁架結(jié)構(gòu)為例驗(yàn)證了所提方法，結(jié)果表明，所提方法能有效克服數(shù)據(jù)噪聲的不利影響，改善方程組的病態(tài)性，獲得比現(xiàn)有方法精度更高的計(jì)算結(jié)果.

1 模型修正方程組

結(jié)構(gòu)模型修正可以利用靜力測(cè)試數(shù)據(jù)或動(dòng)力測(cè)試數(shù)據(jù)來(lái)進(jìn)行，不失一般性，本文利用靜力位移來(lái)建立結(jié)構(gòu)模型修正方程組，簡(jiǎn)述如下：

設(shè)結(jié)構(gòu)初始有限元模型的剛度矩陣為K，在已知的外載荷l的作用下，結(jié)構(gòu)將產(chǎn)生靜力位移.由初始有限元模型可以計(jì)算得到位移的數(shù)值解u0為[7]：

u0=K-1l

(1)

另一方面，通過(guò)實(shí)際靜力測(cè)試，可以獲得結(jié)構(gòu)在外載荷l作用下的位移實(shí)驗(yàn)值ue，由于建模誤差以及測(cè)量誤差，u0和ue總是存在偏差，需要修正剛度矩陣K以使得計(jì)算所得位移值與測(cè)試所得位移值更接近.設(shè)有限元模型中各單元的修正系數(shù)為ci，則修正后的剛度矩陣Km可表示為[7]

(2)

其中Ki為有限元模型中第i個(gè)單元?jiǎng)偠染仃?，N為單元的總數(shù)目.理論上，由修正后的剛度矩陣Km計(jì)算所得的位移值應(yīng)等于實(shí)測(cè)值ue，即

(3)

方程(2)代入(3)，利用Neumann級(jí)數(shù)展開(kāi)并忽略高階項(xiàng)可得：

(4)

方程(1)代入(4)并整理可得：

A·x=y

(5)

A=[η1，…，ηN]，ηi=K-1KiK-1

(6)

x=(c1，…，cN)T

(7)

y=ue-u0

(8)

上述方程中，系數(shù)矩陣A可以通過(guò)初始有限元模型由方程(6)計(jì)算獲得，向量y可以由實(shí)驗(yàn)測(cè)試值及初始有限元模型由方程(8)計(jì)算得到，修正系數(shù)向量x為未知的向量，也是要求解的列向量.因此，模型修正問(wèn)題最終可以歸結(jié)為線性方程(5)的求解問(wèn)題.最小二乘估計(jì)是求解線性方程組最常用的方法，它滿足最優(yōu)線性無(wú)偏性，因此又稱為無(wú)偏估計(jì).最小二乘估計(jì)主要計(jì)算公式[8，9]為：

首先，方程(5)兩邊乘以AT可得法方程為：

z=B·x

(9)

B=ATA

(10)

z=ATy

(11)

其中方陣B稱為法方程的系數(shù)矩陣.由方程(9)可得x的最小二乘估計(jì)為：

xlse=B-1·z

(12)

然而，當(dāng)系數(shù)矩陣B呈現(xiàn)病態(tài)或嚴(yán)重病態(tài)時(shí)，即矩陣B的條件數(shù)很大或矩陣B近似為奇異矩陣，此時(shí)，由方程(12)計(jì)算所得的最小二乘估計(jì)xlse將嚴(yán)重偏離x的真值，必須尋求x的更好估計(jì).

2 一種新的廣義嶺估計(jì)方法

為了解決病態(tài)最小二乘問(wèn)題，許多學(xué)者開(kāi)展了嶺估計(jì)和廣義嶺估計(jì)方面的研究.從矩陣方程的角度來(lái)看，嶺估計(jì)的基本思想是：在系數(shù)矩陣B中增加一個(gè)正則對(duì)角矩陣kI，以降低系數(shù)矩陣的條件數(shù)，從而獲得更加穩(wěn)定估值，其計(jì)算公式為：

xre=(B+kI)-1·z

(13)

其中xre為x的嶺估計(jì)，I為和系數(shù)矩陣B同維的單位矩陣，系數(shù)k稱為嶺參數(shù).如果把方程(13)中的正則矩陣kI更換為一般的對(duì)角矩陣Λ，則可得廣義嶺估計(jì)xgre的計(jì)算公式為：

xgre=(B+Λ)-1·z

(14)

(15)

顯然，對(duì)比方程(12)、(13)和(14)可知，當(dāng)k1=k2=…=kn=k時(shí)，廣義嶺估計(jì)退化為嶺估計(jì)，當(dāng)k1=k2=…=kn=0時(shí)，嶺估計(jì)退化為最小二乘估計(jì).目前，嶺估計(jì)方法的主要問(wèn)題在于選擇嶺參數(shù)往往需要很大的計(jì)算量，且精度不高.為了進(jìn)一步提高計(jì)算精度，并避免嶺參數(shù)選擇所需要的復(fù)雜運(yùn)算，本文提出一種新的廣義嶺估計(jì)方法如下.

首先，分析廣義嶺估計(jì)的計(jì)算公式(14)可知，正則矩陣Λ的作用是兩方面的：其正面作用在于，引入正則矩陣Λ可以降低系數(shù)矩陣的條件數(shù)，即cond(B+Λ)

設(shè)系數(shù)矩陣B中的元素為bij(i=1～n，j=1～n)，即

(16)

找出方程(16)中B對(duì)角元素的最大值，記為bmax，即

bmax=max(bii)，i=1～n

(17)

則廣義嶺估計(jì)中正則矩陣Λ(方程(15))中對(duì)角元素由以下公式來(lái)計(jì)算：

(18)

其中，τ是一個(gè)可以調(diào)節(jié)的底數(shù)，只需滿足τ≥1即可.顯然，如果τ取為1，則方程(18)退化為k1=k2=…=kn=0.05×bmax，相當(dāng)于廣義嶺估計(jì)退化為常規(guī)的嶺估計(jì).實(shí)際應(yīng)用時(shí)一般可以取τ為2～10之間的某個(gè)整數(shù).

最后，將所提的廣義嶺估計(jì)新方法操作步驟總結(jié)如下：首先，由方程(17)找出系數(shù)矩陣B中最大的對(duì)角元素；然后，由方程(18)計(jì)算正則矩陣Λ各對(duì)角元素，得到正則矩陣Λ；最后，由方程(14)計(jì)算未知量x的新廣義嶺估計(jì)解xgre.

3 數(shù)值算例

3.1算例1

首先，以文獻(xiàn)[24]中所模擬的病態(tài)方程組為例，來(lái)驗(yàn)證所提的廣義嶺估計(jì)方法，并將計(jì)算結(jié)果與現(xiàn)有的嶺估計(jì)方法及其他幾種方法進(jìn)行比較，以說(shuō)明本文所提方法的可行性與優(yōu)勢(shì).該病態(tài)方程組的系數(shù)矩陣和觀測(cè)向量具體為：

(19)

該模型法矩陣的條件數(shù)為2.0837×104，病態(tài)性嚴(yán)重.文獻(xiàn)[24]中給出了各種方法的解算結(jié)果，如表1所示，其中TLS表示整體最小二乘法，LSE表示最小二乘法，RE表示嶺估計(jì)方法，VOM表示虛擬觀測(cè)解法.為了衡量表1中各種方法計(jì)算結(jié)果的精度，定義eΔx為x的計(jì)算值與真值之間的偏差，其計(jì)算公式如下：

eΔx=‖xestimate-xtrue‖2

(20)

其中下標(biāo)“2”表示向量的2-范數(shù).

采用本文所提廣義嶺估計(jì)新方法，τ從1取到6時(shí)的計(jì)算結(jié)果如表2所示.由表2可見(jiàn)，τ=4或τ=5時(shí)計(jì)算結(jié)果精度都比較好，進(jìn)一步對(duì)比表1和表2可知，采用本文所提廣義嶺估計(jì)方法，τ取2到6時(shí)各個(gè)解的精度均要好于表1中的現(xiàn)有方法，進(jìn)一步說(shuō)明了本文方法是合理可行的.

3.2算例2

接下來(lái)以圖1所示的桁架結(jié)構(gòu)模型為例，來(lái)驗(yàn)證所提方法在結(jié)構(gòu)模型修正中的應(yīng)用.該結(jié)構(gòu)中桿件材料的彈性模量為E=200 GPa， 密度為

ρ=7.8×103kg/m3，各桿件長(zhǎng)度均為L(zhǎng)=1 m，桿件橫截面面積A=1.759×10-4m2.靜力荷載為F1=F2=F3=10 kN.不失一般性，假設(shè)各桿件修正系數(shù)ci的真值為c2=c7=0.2，c15=0.1，c10=c20=0.15，

表1 現(xiàn)有方法的解算結(jié)果(算例1)

表2 本文所提廣義嶺估計(jì)方法計(jì)算結(jié)果

其他桿件的修正系數(shù)均假設(shè)為0.在由有限元模型計(jì)算所得的位移差向量中添加隨機(jī)噪聲來(lái)模擬測(cè)量誤差，模擬誤差的公式為：

Δu=Δu×(1+δ×unifrnd(-1，1))

(21)

其中δ表示誤差水平，本例中取δ=0.1；unifrnd(-1，1)表示一個(gè)位于[-1，1]區(qū)間里的隨機(jī)數(shù).分別運(yùn)用最小二乘法、嶺估計(jì)法和本文所提廣義嶺估計(jì)新方法，求解模型修正線性方程(5)，計(jì)算所得結(jié)果列于圖2中.

由圖2可見(jiàn)，和修正系數(shù)的真值(即假設(shè)值)相比，本文所提方法的計(jì)算結(jié)果與之最接近.另外，也可利用公式(20)來(lái)衡量各計(jì)算結(jié)果與真值之間的偏差，三種方法計(jì)算結(jié)果的偏差分別為：0.548 8，0.256 6和0.178 0，也說(shuō)明了本文方法計(jì)算精度最好.因此，本文所提方法可以應(yīng)用于結(jié)構(gòu)模型修正問(wèn)題中，能夠一定程度上克服數(shù)據(jù)噪聲的不利影響，改善方程組的病態(tài)性，提高求解精度.

4 結(jié)論

針對(duì)病態(tài)最小二乘問(wèn)題，本文提出一種新的廣義嶺估計(jì)方法，并將其應(yīng)用于結(jié)構(gòu)有限元模型修正問(wèn)題中.所提方法充分考慮系數(shù)矩陣對(duì)角元素之間的差異，提出一個(gè)底數(shù)可調(diào)節(jié)的正則系數(shù)計(jì)算公式，用于獲得廣義嶺估計(jì)中的正則對(duì)角矩陣，從而盡量減少?gòu)V義嶺估計(jì)中由于引入正則矩陣而帶來(lái)的額外誤差，更好地兼顧了解的穩(wěn)定性和準(zhǔn)確性. 以一個(gè)病態(tài)方程組和一個(gè)桁架結(jié)構(gòu)為例驗(yàn)證了所提方法，結(jié)果表明，所提方法能有效克服數(shù)據(jù)噪聲的不利影響，改善方程組的病態(tài)性，獲得比現(xiàn)有方法精度更高的計(jì)算結(jié)果.

圖1 桁架結(jié)構(gòu)及其靜力加載

圖2 各種方法計(jì)算結(jié)果對(duì)比(算例2)