秦玉靈，孔憲仁，羅文波

（1.哈爾濱工業(yè)大學 衛(wèi)星技術研究所，哈爾濱 150001；2.北京空間飛行器總體設計部，北京 100094）

0 引言

有限元模型修正技術在航空航天工程中應用廣泛，高精度的有限元模型是準確進行結構力學特性分析的基礎。但由于建模誤差及實際結構在使用過程中的損傷等因素的影響，有限元模型計算結果與實測結果間總存在差異。利用結構實測響應修正有限元計算響應，使得修正后有限元模型計算響應值與試驗測量值一致的過程即為有限元模型修正過程[1]。修正過程充分利用結構試驗和有限元分析兩者優(yōu)點，用少量的試驗數據對有限元模型進行修正，修正后模型可以代替實際結構進行多種分析，節(jié)約試驗成本和時間。

用有限元軟件的優(yōu)化功能直接進行模型修正時，由于每次修正量的變化都要重新調用有限元軟件計算，因此計算效率較低。近年來響應面方法開始在模型修正領域被廣為應用，其基本思想是用多項式等構造響應面對原結構進行模擬，以較少的試驗代價得到可接受的結果，在工程中有重要應用價值[2-3]。響應面模型構造過程需要對模型中的待修正參數進行篩選，用正交設計方法設計試驗，并用統(tǒng)計分析方法中的方差分析選擇模型參數中對響應影響程度大的參數作為待修正參數。響應面模型形式簡潔，易于與粒子群算法等優(yōu)化算法融合，修正效率高，但一組響應面只能對應一組分析過程，因此將響應面修正后參數代入原模型，可以進行多種力學特性分析，提高了分析效率。

本文用某雷達衛(wèi)星簡化有限元模型計算基于正交設計的各水平參數下模態(tài)頻率，用方差分析確定待修正參數并構造二次響應面模型，用最小二乘法確定多項式系數，以響應面計算結果與實測結果的差值構造適應度函數并用其引導含混沌搜索機制的改進粒子群算法對待修正參數的偏移量進行尋優(yōu)，修正后參數代入原有限元模型得到修正模型。修正后模型計算所得測試頻段內模態(tài)頻率與基準模型一致，且能以一定精度預測測試頻段外的模態(tài)頻率，從而證實了基于響應面和改進粒子群算法的模型修正方法的有效性。

1 正交試驗設計及響應面構造

正交試驗設計方法是一種研究多因子多水平試驗問題的重要數學方法。該方法根據正交性從全面試驗中挑選出“均勻分散，整齊可比”的樣本點，是一種高效經濟的試驗設計方法[4]。正交試驗完成后可用方差分析評價參數顯著性，從中選擇顯著性高的參數進行后續(xù)計算。

方差分析（analysis of variance，ANOVA）又稱“變異數分析”或“F檢驗”，其基本思想是將數據總離差平方和分解為各因素及誤差的離差平方和，再用F檢驗確定各因素影響的顯著性水平，從而確定待修正參數。模型參數A的F值檢驗按公式(1)進行[5]。

式中：α為置信度；FA為因素A 的顯著性水平；SA和Se為因素A和誤差e的偏差；fA為因素A和誤差e的自由度。

響應面法是用多元多項式或非多項式模型（如人工神經網絡）來描述系統(tǒng)自變量和響應特征的復雜關系，從而替代有限元仿真和其他復雜模型進行更有效設計或計算的一種方法。響應面表達式一方面要盡可能簡單，方便計算；另一方面又要考慮到能靈活反映各種真實曲面形狀，擬合精度和效率是評價響應面建模方法好壞的重要標準[6]。二階多項式響應面形式簡潔，準確性高，在工程模擬中應用廣泛，其形式為

響應面模型質量可用兩個標準檢驗，即均方根誤差（root mean squared error，RMSE）的相對值和決定系數R2。

2 引進混沌搜索機制的粒子群算法

粒子群算法（particle swarm optimization，PSO）通過模擬鳥群覓食過程中的遷徙和聚集行為實現對問題的優(yōu)化。算法中的每一個粒子都代表搜索空間的一個解，先隨機初始化每個粒子的位置和速度，然后粒子在自身及整個種群最優(yōu)位置的引導下飛向最優(yōu)解[7-9]。標準粒子群算法中粒子速度和位置進化公式為

式中：pi,d(t)和pgd(t)分別為第t次迭代時第i個粒子第d維分量的最優(yōu)位置和群體最優(yōu)位置；常量c1＞0，c2＞0（c1+c2≈4），分別代表對自身經歷的最優(yōu)位置和種群最優(yōu)位置的記憶能力；r1、r2用來增強粒子飛行的隨機性以保持種群的多樣性，r1∈rand(0,1)，r2∈rand(0,1)；ω 為慣性因子，ω＞0，當 ω較大時群體以較大步長進行全局搜索，較小時算法以較小步長進行細致局部搜索。

混沌搜索的基本思想是：把混沌變量線性映射到原設計空間，然后利用混沌變量進行有效搜索。由于混沌的遍歷性，混沌優(yōu)化算法易跳出局部最優(yōu)解，是一種很好的搜索機制[10-11]。Logistic映射是一種常用的混沌映射，其表達式為

式中：μ為控制變量，一般在(3.5699456,4]之間取值；cx∈(0,1)。μ= 4時為完全混沌狀態(tài)，此時cx在(0,1)范圍內遍歷。Logistic映射有3個不動點——0.25，0.5，0.75，在參數設定過程中應避免使用這些不動點。

混沌變量cx與普通變量x可以通過式(7)和式(8)進行往復映射。

式中：minx為設計空間中變量x的取值下限；maxx為取值上限。

引入混沌搜索機制的改進粒子群算法流程如圖1所示。

圖1 引入混沌搜索機制的改進粒子群算法流程Fig.1 The improved PSO with chaos-search mechanism

3 衛(wèi)星有限元模型修正

基于響應面方法和改進粒子群算法的有限元模型修正包括響應面模型的建立與粒子群算法的尋優(yōu)兩方面內容。其模型修正流程如圖2所示。

圖2 基于響應面方法的模型修正流程Fig.2 Flow chart of the RSM-based model updating

1）用等效板理論計算構造衛(wèi)星壁板和太陽帆板的蜂窩板的等效結構參數；

2）用正交試驗設計方法進行試驗設計，通過F值檢驗判定各結構參數對頻率響應的影響，從而確定用于構造響應面的待修正參數；

3）用ANSYS中的SHELL63單元和BEAM4單元建立簡化衛(wèi)星模型并計算其在各組參數下的模態(tài)響應；

4）構造二次多項式響應面，參數歸一化并用最小二乘法則確定響應面函數系數向量；

5）響應面模型有效性判定，若有效則進入下一步分析，否則返回第2步；

6）用響應面計算所得頻率?y與有限元計算頻率y的偏差構造適應度函數；

7）用帶混沌搜索機制的粒子群算法修正結構參數；

8）判斷是否滿足終止條件，若滿足則輸出修正量，否則返回第7步繼續(xù)計算；

9）將修正量代入有限元模型，模型計算質量得到改善，進行多種分析計算。

3.1 衛(wèi)星簡化有限元模型建立

某雷達衛(wèi)星簡化有限元模型如圖3所示，衛(wèi)星壁板和太陽帆板為蜂窩板結構，其等效板由ANSYS中的SHELL63殼單元建立，由等效板理論[12]計算所得壁板和帆板的等效參數均為：蜂窩板尺寸1 m×1 m，板厚req=0.06 m，等效彈性模量Eeq=831.4 MPa，等效泊松比μeq=0.3，等效密度ρeq=107 kg/m3。衛(wèi)星箱體與蜂窩板之間由4根梁進行連接，用BEAM4梁單元建立，結構參數為：彈性模量E=70 GPa，密度 ρ=2.7×103kg/m3，泊松比 μ=0.3，截面積a =10-4m2，z向和y向抗彎慣性矩Izz=Tyy=10-8m4。

圖3 衛(wèi)星簡化有限元模型Fig.3 Simplified FEM of the satellite

3.2 正交試驗設計及待修正參數確定

對有限元模型中的4個等效結構參數req、Eeq、ρeq、μeq進行篩選，同時考慮結構損傷及阻尼材料的影響，將連接梁的3個材料參數E、ρ、a也作為待修正參數，每個參數取2個水平，即七因素二水平，水平1為7個參數原值，水平2為7個參數攝動后值（各減小20%）。正交設計方法設計試驗時，首先用正交表從所有可能搭配中挑出若干必需的試驗，然后再用統(tǒng)計分析方法對試驗結果進行綜合處理，可以減少試驗次數，提高分析效率。本例選取正交表L8(27)設計8組試驗數據如表1所示。

表1 不確定參數正交設計表Table 1 Orthogonal design for uncertain parameters

用有限元分別計算8組參數前6階模態(tài)頻率值，對計算所得的48個樣本點進行方差分析。常用F(f,f )值檢驗水平有4個，即α =[0.25,0.10,αAe 0.05,0.01]，各水平意義和符號見表2。方差分析過程中發(fā)現，ρ和a兩個參數對模態(tài)頻率影響最小，故將兩者視為誤差項，用于對其余5個參數的顯著性水平進行計算。

表2 顯著性水平Table 2 Significance level

由表3的F值檢驗結果可知，req、Eeq和 ρeq對模態(tài)頻率影響程度最高，E次之，μeq影響力最弱，由此選定板等效厚度req、等效彈性模量Eeq、等效密度ρeq及梁彈性模量E為待修正量。

表3 各參數顯著性水平Table 3 Significance level of parameters

3.3 響應面構造及有效性檢驗

用F值檢驗確定的4個結構參數req、Eeq、ρeq和E構造響應面的完全二次多項式為

響應面模型有效性由相對RMSE和決定系數R2評價，如圖4所示。

圖4 響應面有效性評價Fig.4 Validity evaluation of the response surface

由圖4可知，衛(wèi)星前6階模態(tài)的響應面計算值與相同結構參數的有限元計算值之間的RMSE均趨于0，決定系數R2均趨于1，響應面模型有效，可用于進一步分析計算。各階模態(tài)擬合精度具體分析為：f1和f5擬合精度最高，原因在于對這兩階模態(tài)有影響的結構參數req、Eeq、ρeq和E都參與了響應面的構建；f3和f4擬合精度相對較差，原因在于對這兩階模態(tài)有影響的參數μeq未參與響應面構建。由此可見響應面構建中參數選擇的重要性，在工程實際中可以對某重要模態(tài)單獨分析其影響參數，構建相應響應面模型。

3.4 衛(wèi)星模型修正及有效性檢驗

以參數水平(1,1,1,1)即x1=[req1,Eeq1,ρeq1,E1]所建有限元模型為基準模型M1，以參數水平(2,2,2,2)即x2=[req2,Eeq2,ρeq2,E2]所建有限元模型為待修正模型M2，兩組結構參數的關系為x2=x1(1-20%)=x1(1+Δx)， 其 中 Δx=[Δreq,ΔEeq,Δρeq,ΔE]為待修正量的偏差量，亦即粒子群算法搜索的目標解。模型結構參數如表4所示。

表4 模型結構參數Table 4 Structure parameters of the FEM

響應面模型修正中的適應度函數由基準模型M1計算所得模態(tài)頻率f1和響應面模型S計算所得模態(tài)頻率f?之間的差值構成，其可表示為

式中：N=6為模態(tài)頻率階數，該適應度函數的意義為使得各階模態(tài)頻率相對誤差和最?。籪i1（i=1,2,...6）為定值（基準模型M1計算所得模態(tài)頻率）；（i=1,2,...6）為響應面模型計算值，其值隨待修正量x=[req,Eeq,ρeq,E2]的不斷變化（由粒子群算法迭代搜索獲得）而變化，至適應度函數值達到設定條件時，待修正量得到確定。

用帶混沌搜索機制的改進粒子群算法對攝動后結構參數進行修正，每個粒子維數為D=4（對應4個待修正量），粒子群種群規(guī)模為N=50，即用50個粒子進行隨機搜索，迭代次數設為Tmax=30。搜索空間范圍（粒子位置x和速度v）由以下分析確定：任取待修正參數xi，設定其變化范圍為原值（x1=[req1,Eeq1,ρeq1,E1]）上下20%浮動，響應面構造過程中進行參數歸一化得

式中，ximax=(1+20%)xi1；ximin=(1-20%)xi1；為各水平參數平均值，則xi=(xi1+xi2)/2=0.9xi1?？傻?0.25≤x?i≤0.75，由此，粒子位置x和速度v變動范圍-0.25≤x≤0.75，-0.25≤v≤0.75。

由于算法的隨機性，用改進粒子群算法重復進行20次搜索，取各次搜索的適應度函數平均值繪制平均適應度函數收斂曲線；同時選各次搜索所得最優(yōu)解的平均值作為修正量Δx?，將Δx?從歸一化后的參數空間映射回原設計空間即可得到待修正參數的真實攝動量Δx；將Δx代入響應面模型和待修正的有限元模型得到修正后的響應面模型S1和修正后有限元模型M3；計算修正后模型M3的前6階模態(tài)頻率與基準模型M1比較，評價修正后模型對測試頻段內模態(tài)頻率的復現能力，同時為了考察修正后模型的預測能力，計算M3的7～10階模態(tài)頻率并與基準模型M1對比，分析修正后模型對測試頻段外模態(tài)頻率的預測效果。

由表5可知，用修正后參數建立的響應面模型S1和相同參數建立的有限元模型M3計算所得模態(tài)頻率接近，證實了響應面模型代替有限元模型計算的可行性。由表6可知，修正后有限元模型M3的結構參數較待修正模型M2更接近基準模型M1，模型質量得到改善。同時由表5和圖5可知：用攝動后所得待修正模型M2計算獲得的模態(tài)頻率明顯偏離基準模型M1計算值，修正后有限元模型M3在測試頻段內計算所得的前6階模態(tài)頻率與基準模型誤差在1%以內（第3、4階頻率誤差較大為0.93%和1.02%，與第3、4階模態(tài)響應面構造中參數選擇相對誤差較大有關），測試頻段外誤差在3%～5%之間，大于測試頻段內誤差，但仍在可接受精度范圍內。因此證實了響應面方法修正后的有限元模型不但能以高精度復現測試頻段內模態(tài)頻率，而且能以可接受的精度對測試頻段外頻率進行預測。

圖6描述了重復執(zhí)行帶混沌搜索機制的改進粒子群算法20次迭代且每次迭代30個循環(huán)過程中的平均適應度函數收斂曲線。適應度值隨迭代進行不斷減小，意即修正后模型計算頻率與實測頻率之間差值隨著迭代搜索過程的進行不斷減小，粒子群在群體最優(yōu)粒子的引導下逐漸向全局最優(yōu)解靠近。將搜索到的最優(yōu)解代入原有限元模型，所得模型計算模態(tài)頻率能充分靠近實測值，證實了算法的有效性。

表5 響應面和有限元模型計算模態(tài)頻率Table 5 Modal frequencies of the RSM and FEM

表6 模型參數對比Table 6 Comparison of the model parameters

圖5 模態(tài)頻率對比Fig.5 Comparison of the modal frequencies

圖6 平均適應度函數收斂曲線Fig.6 Mean convergence curve of the fitness function

4 結論

1）響應面模型計算所得的模態(tài)頻率接近由相同結構參數的有限元模型計算所得頻率，證實了響應面方法的可行性和可靠性。響應面模型便于與粒子群算法等優(yōu)化算法相結合，用響應面模型代替有限元模型可以提高修正效率和模擬精度。

2）基于正交試驗設計方法的響應面模型構造時，用F值檢驗結果作為選擇待修正量的依據，待修正量的選取是決定響應面精度的重要因素，應通過改進試驗設計方法或對現有方法增加補充條件以提高待修正量選取的準確性。

3）響應面模型的局限性在于一組響應面只能模擬一組響應數據與結構參數間的關系，要獲得其他響應數據，必須再次通過多組試驗數據獲得另一組響應面，無法像有限元模型那樣，通過一組結構參數，可以進行多種分析運算。因而用響應面模型修正所得結構參數重新建立有限元模型進行分析，可提高修正效率和計算效率。

（References）

[1]朱宏平, 徐斌, 黃玉盈.結構動力模型修正方法的比較研究及評估[J].力學進展, 2002, 32(4): 513-524 Zhu Hongping, Xu Bin, Huang Yuying.Comparison and evaluation of analytical approaches to structural dynamic model correction[J].Advances in Mechanics, 2002, 32(4):513-524

[2]任偉新, 陳華斌.基于響應面的橋梁有限元模型修正[J].土木工程學報, 2008, 41(12): 73-78 Ren Weixin, Chen Huabin.Response-surface based on finite element model updating of bridge structure[J].China Civil Engineering Journal, 2008, 41(12): 73-78

[3]Sobieski I P, Kroo I M.Collaborative optimization using response surface estimation[J].AIAA Journal, 2000,38(10): 1931-1938

[4]郭勤濤, 張令彌, 費慶國.用于確定性仿真的響應面法及其試驗設計研究[J].航空學報, 2006, 27(1): 55-61 Guo Qintao, Zhang Lingmi, Fei Qingguo.Response surface method and its experimental design for deterministic computer simulation[J].Acta Aeronautica et Asrtonautica Sinica, 2006, 27(1): 55-61

[5]費慶國, 張令彌, 李愛群, 等.基于統(tǒng)計分析技術的有限元模型修正研究[J].振動與沖擊, 2005, 24(3): 23-26 Fei Qingguo, Zhang Lingmi, Li Aiqun, et al.Finite element model updating using statistics analysis[J].Journal of Vibration and Shock, 2005, 24(3): 23-26

[6]Allen M S, Cambeves J A.Comparison of uncertainty propagation/ response surface techniques for two aeroelastic systems, AIAA 2009-2269[R]

[7]秦玉靈, 孔憲仁, 羅文波.GA-PSO組合算法模型修正[J].航天器環(huán)境工程, 2009, 26(4): 383-385Qin Yuling, Kong Xianren, Luo Wenbo.The model updating by GA-PSO algorithm[J].Spacecraft Environment Engineering, 2009, 26(4): 383-385

[8]孫木楠, 史志俊.基于粒子群優(yōu)化算法的結構模型修改[J].振動工程學報, 2004, 17(3): 350-353 Sun Munan, Shi Zhijun.Structural model updating based on particle swarm optimization[J].Journal of Vibration Engineering, 2004, 17(3): 350-353

[9]Venter G, Sobieszczanski-Sobieski J.Particle swarm optimization[J].AIAA Journal, 2003, 41(8): 1583-1589

[10]林昆, 馮斌, 孫俊.混沌量子粒子群優(yōu)化算法[J].計算機工程與設計, 2008, 29(10), 2610-2612 Lin Kun, Feng Bin, Sun Jun.Chaos quantum-behaved particle swarm optimization algorithm[J].Computer Engineering and Design, 2008, 29(10): 2610-2612

[11]侯力, 王振雷, 錢鋒.基于混沌序列的自適應粒子群優(yōu)化算法[J].計算機工程, 2008, 34(18): 210-212 Hou Li, Wang Zhenlei, Qian Feng.Adapting particle swarm optimization algorithm based on chaotics series[J].Computer Engineering, 2008, 34(18): 210-212

[12]夏利娟, 金咸定, 汪庠寶.衛(wèi)星結構蜂窩夾層板的等效計算[J].上海交通大學學報, 2003, 37(7): 999-1001 Xia Lijuan, Jin Xianding, Wang Xiangbao.Equivalent analysis of honeycomb sandwich plates for satellite structure[J].Journal of Shanghai Jiaotong University,2003, 37(7): 999-1001