楊秋偉 陳 華 李翠紅

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

0 引言

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

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

1 模型修正方程組

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

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

u0=K-1l

(1)

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

(2)

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

(3)

方程(2)代入(3)，利用Neumann級數(shù)展開并忽略高階項可得：

(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可以通過初始有限元模型由方程(6)計算獲得，向量y可以由實驗測試值及初始有限元模型由方程(8)計算得到，修正系數(shù)向量x為未知的向量，也是要求解的列向量.因此，模型修正問題最終可以歸結(jié)為線性方程(5)的求解問題.最小二乘估計是求解線性方程組最常用的方法，它滿足最優(yōu)線性無偏性，因此又稱為無偏估計.最小二乘估計主要計算公式[8，9]為：

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

z=B·x

(9)

B=ATA

(10)

z=ATy

(11)

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

xlse=B-1·z

(12)

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

2 一種新的廣義嶺估計方法

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

xre=(B+kI)-1·z

(13)

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

xgre=(B+Λ)-1·z

(14)

(15)

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

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

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

(16)

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

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

(17)

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

(18)

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

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

3 數(shù)值算例

3.1算例1

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

(19)

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

eΔx=‖xestimate-xtrue‖2

(20)

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

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

3.2算例2

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

ρ=7.8×103kg/m3，各桿件長度均為L=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 本文所提廣義嶺估計方法計算結(jié)果

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

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

(21)

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

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

4 結(jié)論

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

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

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