李建宇,任 朝,張麗麗,魏凱杰
(1.天津科技大學(xué) 機械工程學(xué)院,天津市輕工與食品工程機械裝備集成設(shè)計與在線監(jiān)控重點實驗室,天津 300222;2.天津職業(yè)技術(shù)師范大學(xué) 理學(xué)院,天津 300222)
模型參量反求是一種利用試驗測量數(shù)據(jù)反求力學(xué)模型中的未知參量的方法。對于涉及多個物理參數(shù)的復(fù)雜力學(xué)模型,如復(fù)合材料和彈塑性等,模型參量反求技術(shù)是獲取模型參數(shù)的有效手段,近年來受到研究者的關(guān)注[1]。De wilde等[2]較早地提出結(jié)合數(shù)值計算與試驗測量來反求材料參量的方法。Ghouati等[3]采用集成有限元和優(yōu)化技術(shù)的方法來測量材料參數(shù)。國內(nèi)學(xué)者在模型反求算法方面也開展了系統(tǒng)而深入的工作[4-6]。
目前,數(shù)字圖像相關(guān)DIC(Digital Image Correlation)技術(shù)[7]的發(fā)展已經(jīng)使得對結(jié)構(gòu)變形信息的全場測量成為可能,不僅為直接觀測結(jié)構(gòu)的力學(xué)響應(yīng)提供了便利,也為利用變形測量信息進一步研究結(jié)構(gòu)和材料的物理性能提供了數(shù)據(jù)支撐。將DIC測量數(shù)據(jù)與模型參量反求技術(shù)結(jié)合并應(yīng)用于具體工程問題的研究受到關(guān)注,如馬源等[8]對風(fēng)電葉片的彈性參數(shù)進行了反求,楊思滿等[9]研究了車用復(fù)合材料的參數(shù)反演問題。
DIC測量技術(shù)的發(fā)展有效解決了結(jié)構(gòu)變形數(shù)據(jù)獲取的問題,但是如何更有效地使用DIC變形數(shù)據(jù)又成為新的問題。由于數(shù)字圖像包含的數(shù)據(jù)量極大,如果在模型參量反求中使用全部的測量數(shù)據(jù),不僅會導(dǎo)致計算耗時長及收斂速度慢等問題,而且也沒有必要。現(xiàn)有文獻還很少有關(guān)于如何有選擇性地使用部分DIC數(shù)據(jù)進行力學(xué)模型參量反求的研究。已有的工作多采用均勻抽樣的方法從測量數(shù)據(jù)中間隔抽取離散的數(shù)據(jù)點[10,11]進行計算,這種做法的問題是數(shù)據(jù)點之間的間距不好把握。倘若間隔過大,則抽取數(shù)據(jù)量少,無法確保較為全面地獲得測量場的主要信息,這在具有變形奇異的問題中較為突出,如非標準幾何形狀鋁板的變形[12]和剪切帶的滑移問題[13],數(shù)據(jù)量太少容易丟失局部信息;反之,若間隔過小,則抽取數(shù)據(jù)量增加,又為后續(xù)反求計算帶來困難。
在圖像識別和機器學(xué)習(xí)領(lǐng)域[14],常采用數(shù)據(jù)壓縮的手段來處理大數(shù)據(jù)問題。海量數(shù)據(jù)經(jīng)數(shù)據(jù)壓縮后,雖然存在一定的信息損失,但壓縮后的數(shù)據(jù)卻能夠保持原數(shù)據(jù)的主要特征。目前,常用的數(shù)據(jù)壓縮技術(shù)有主元分析法、局部保持投影法和拉普拉斯特征投影法等[14],其中主元分析法是最基本的線性數(shù)據(jù)壓縮方法。本文針對DIC測量數(shù)據(jù)量龐大和選點困難的問題,提出采用基于主元分析的數(shù)據(jù)降維方法對DIC數(shù)據(jù)進行壓縮處理,并在此基礎(chǔ)上給出一種模型參量反求方法。一方面,確保DIC變形場測量的重要特征不會因壓縮而丟失;另一方面,有助于顯著節(jié)省數(shù)據(jù)存儲和計算時間。
在各種非接觸全場光測力學(xué)方法中,DIC以其具有實驗設(shè)備及測量過程簡單、環(huán)境適應(yīng)性強、測量精度高和適用范圍廣等優(yōu)勢得到了廣泛應(yīng)用,已經(jīng)作為一種常用而有效的表面變形測量手段應(yīng)用于實驗固體力學(xué)領(lǐng)域[13,15]。
DIC測量的基本原理是匹配物體表面不同狀態(tài)下的數(shù)字化散斑幾何點,追蹤點的運動得到物體表面的變形信息。散斑場上的散斑是隨機分布的,將每個散斑點周圍的小區(qū)域稱為子區(qū)域。散斑場上以某一點為中心的子區(qū)域可作為該點位移的載體,通過搜索和分析該子區(qū)域的移動和變化,便可獲得該點的位移。如圖1所示,以特征點P為中心的參考區(qū)域S為子區(qū)域,其大小為m×n個像素,當(dāng)發(fā)生運動后,子區(qū)域S移動到子區(qū)域S1的位置,由此可以得到特征點P變形后的位置P1,從而得到P點的位移。對于高質(zhì)量的散斑圖,DIC的位移測量精度可達到0.01像素,圖像數(shù)據(jù)量在十萬級以上,極其龐大。
圖1 特征點P的位移
Fig.1 Displacement of characteristic pointP
主元分析PCA(Principal Component Analysis)的基本原理是從數(shù)據(jù)中識別主要特征,將高維數(shù)據(jù)通過線性變換投影到低維空間,再通過特征值分析,確定需要保留的主元個數(shù),實現(xiàn)數(shù)據(jù)的壓縮。主元分析經(jīng)常作為一種數(shù)據(jù)集預(yù)處理技術(shù),在數(shù)據(jù)應(yīng)用到其他算法之前使用,能夠去除冗余信息和噪聲,使數(shù)據(jù)變得更加簡單高效,提高計算效率。但是,將主元分析技術(shù)用于壓縮DIC測量數(shù)據(jù)并用于力學(xué)模型的參量反求的研究還鮮有報道。
給定一DIC測量位移場,包括m×n個測點,其中m為測點行數(shù),n為測點列數(shù)。所測量位移構(gòu)成二維矩陣Um n,如式(1)。
(1)
式中列向量
ui=[u1i…um i]T(i=1,…,n) (2)
利用主元分析方法對式(1)所示的位移場進行壓縮處理,具體步驟如下。
(1) 數(shù)據(jù)中心化處理
(3)
式中 mean為按列計算均值、repmat為將矩陣mean(Um n)復(fù)制m×1塊。
(2) 協(xié)方差矩陣特征值分解
(4)
式中Φ為特征向量,λ為特征值。
(3) 數(shù)據(jù)壓縮,截取前r (5) (6) 累計方差解釋率[16]: (7) r的取值通常要使得累計方差解釋率達80%以上,從而保證包含原始數(shù)據(jù)的主要信息特征。 利用上述PCA方法對DIC測量位移數(shù)據(jù)進行壓縮處理。試驗采用MTS E45萬能試驗機和美國CSI公司的DIC儀器。DIC儀器能最高處理55000個節(jié)點數(shù)據(jù),包括2個400萬像素的CCD相機、2個照明光源和Vic -3D測量系統(tǒng)。數(shù)據(jù)采集過程如圖2所示,設(shè)定測量區(qū)域,設(shè)置子區(qū)為 29 pixel×29 pixel,步長為7 pixel,在試件測量平面內(nèi)選擇208個數(shù)據(jù)點采集位移信息。 圖2 DIC數(shù)據(jù)采集過程 Fig.2 DIC data acquisition process 對于設(shè)定的測量區(qū)域,位移信息可以用一個二維矩陣表示,以鉛垂方向的位移分量為例,表面測點數(shù)據(jù)為26×8,利用主元分析對該數(shù)據(jù)矩陣進行壓縮處理,選擇第一主元(26×1)即能滿足精度要求。數(shù)據(jù)量由208個降低為26個。圖3為原始數(shù)據(jù)和第一主元表征數(shù)據(jù)的三維形貌對比。 綜上所述,位移數(shù)據(jù)通過PCA壓縮后具有以下特征。 (1) 存儲數(shù)據(jù)規(guī)模降低,原始數(shù)據(jù)維數(shù)由m×n壓縮為m×r。對于結(jié)構(gòu)變形問題,數(shù)據(jù)矩陣的主要信息集中在前幾個主元,r取值1~3即能滿足累計方差解釋率大于0.8的要求。 (2) 處理后的數(shù)據(jù)更趨于平滑,具有降噪的功能,可有效消除測量誤差的影響。 利用最小二乘法進行模型參量反求的基本數(shù)學(xué)模型為 (8) 式中A為待求的模型參量,uC(A)為模型計算位移值,uM為實際測量位移值。在上述模型中,測量數(shù)據(jù)的選擇對于模型參量的反求計算影響甚大。若選擇的數(shù)據(jù)點個數(shù)過多,容易導(dǎo)致計算耗時長和收斂速度慢的問題;若選取的數(shù)據(jù)點數(shù)量偏少,可能導(dǎo)致信息丟失和過擬合等問題。為此,本文提出一種基于數(shù)據(jù)PCA壓縮的模型參量反求方法。 圖3 原始數(shù)據(jù)和第一主元表征數(shù)據(jù)的三維形貌對比 Fig.3 Comparison of three -dimensional morphology between original data and first principal component characterization data 對數(shù)據(jù)進行PCA處理在實質(zhì)上是將原始數(shù)據(jù)在一組正交基向量{Φj}構(gòu)成的坐標系統(tǒng)中進行投影展開,即μi在基向量Φj上的投影坐標為zi j。而壓縮的含義為通過截取前幾階特征值較大的分量來近似逼近原始數(shù)據(jù)。因此,將直接采用測量數(shù)據(jù)轉(zhuǎn)為通過比較實測數(shù)據(jù)與模型計算數(shù)據(jù)在{Φj}上的投影坐標值{zi j}來實現(xiàn)模型參量的反求。進一步,經(jīng)PCA數(shù)據(jù)壓縮后,僅需比較前幾階坐標投影值,就可減少目標函數(shù)項數(shù)?;谠撍悸罚谧钚《朔ǖ目蚣芟?,構(gòu)建如下模型參量反求模型。 (9) 式中 zC r(A)=uC(A)Φ,uC(A)=zC r(A)ΦT 進入10月下旬,隨著南方最后一小部分地區(qū)秋季肥掃尾工作的順利結(jié)束,全國的秋季肥市場也徹底畫上了句號。從現(xiàn)階段開始,到11月冬儲政策發(fā)布之前,市場將會處在短暫的空檔期。目前復(fù)合肥主流出廠報價:45%氯基復(fù)合肥在2000-2200元/噸,45%硫基復(fù)合肥在2300-2500元/噸。 (10) zM r=uMΦ,uM=zM rΦT (11) Φ=[Φ1…Φr],Φj=[Φj1…Φj n]T (j=1,…,r) (12) (13) 由式(10,11),根據(jù)PCA截斷原理[16] (14) 式中 第一個等式為PCA全項正交展開,第二個約等號為PCA截取前r階特征。再利用{Φj}的正交性得 (15) 模型參量反求模型的求解,目前主要包括兩類算法,一類為基于梯度的局部優(yōu)化算法,如Gauss-Newton法[17]和L-M算法等;另一類為基于智能進化的全局優(yōu)化算法,如蟻群算法、遺傳算法和粒子群法等。本文采用Gauss-Newton法求解式(8,9),其中梯度值采用差分法計算。Gauss-Newton法對初值較為敏感,計算過程中所選的初值應(yīng)該在解的附近。考慮到本文研究是利用實驗數(shù)據(jù)反求模型參數(shù),在計算過程中選用所研究材料參數(shù)的標稱值為初值,這在一定程度上保證了初值與解相距不遠。圖4給出基于數(shù)據(jù)PCA壓縮的模型參量反求流程。 選取碳素結(jié)構(gòu)鋼Q235和纖維復(fù)合材料層合板作為試驗對象,應(yīng)用萬能試驗機進行單軸拉伸試驗和三點彎曲試驗。通過DIC測量試驗過程中的位移信息,反求材料的彈性模量和泊松比等參數(shù)[15,18]。 針對金屬材料拉伸試驗的模型參量反求,在有限元軟件中進行參量化建模,模型如圖5所示。設(shè)定初始材料參數(shù)彈性模量E和泊松比μ,邊界條件為一端固定,另一端施加位移載荷,單元類型為C3D8I(8節(jié)點三維實體,非協(xié)調(diào)單元),根據(jù)網(wǎng)格大小確定數(shù)據(jù)量,采集平面內(nèi)的全場數(shù)據(jù)點,獲取金屬材料在拉伸過程的位移信息。 在選擇原始數(shù)據(jù)量大且累計方差解釋率高的208個數(shù)據(jù)量時,圖6給出了直接反求法和數(shù)據(jù)壓縮反求法的結(jié)果對比,圖7展示了兩個材料參數(shù)的迭代過程??梢钥闯?,目標函數(shù)值下降較快,本文所提模型(9)的目標函數(shù)值和待求模型參量均逼近經(jīng)典模型(8)的目標函數(shù)值和待求參量值。 試驗表明,直接反求法和數(shù)據(jù)壓縮反求法的迭代過程呈現(xiàn)較大差異。從圖8可以看出,隨數(shù)據(jù)量的增大,本文方法的收斂迭代次數(shù)較直接反求法的迭代次數(shù)明顯減少。而且,本文方法的收斂迭代次數(shù)不會隨著數(shù)據(jù)量的增加而顯著增加,而是趨于穩(wěn)定。 圖4 模型參量反求流程 Fig.4 Flow chart of inverse of model parameters 圖5 Q235拉伸試件 Fig.5 Q235 tensile testspecimen 結(jié)果表明,經(jīng)多次迭代計算,數(shù)據(jù)PCA壓縮的模型參量計算反求方法實現(xiàn)了反求值快速收斂于參考值,相較直接反求法,兩者獲得待識別的材料參數(shù)都具有較高精度。 針對復(fù)合材料三點彎曲試驗的模型參量反求,加入需要反求的參數(shù)形式。在有限元軟件中的模型及參數(shù)如圖9所示,其中E11,E22和E33分別為11,22和33方向的彈性模量;G12,G13和G23分別為12,13和23方向的剪切模量;μ12,μ13和μ23為泊松比。設(shè)定初始材料參數(shù),邊界條件為下端兩側(cè)鉸支,上端施加位移載荷,單元類型為C3D8R(8節(jié)點三維實體,縮減積分單元),根據(jù)單元尺寸確定采集平面內(nèi)的全場數(shù)據(jù)量,獲取復(fù)合材料在三點彎曲過程的位移信息。 圖6 兩種方法的目標函數(shù)值對比(208個原始數(shù)據(jù)) Fig.6 Comparison of objective function by two methods(208 data) 圖7 2個材料參數(shù)迭代過程對比(208個原始數(shù)據(jù)) Fig.7 Comparison of iterative process of two material parameters (208 data) 圖8 數(shù)據(jù)量對迭代次數(shù)的影響 Fig.8 Effects of data volume on iteration times 在選擇原始數(shù)據(jù)量大且累計方差解釋率高的2667個數(shù)據(jù)下,圖10給出直接反求法和數(shù)據(jù)壓縮反求法的結(jié)果對比。圖11展示了其中四個材料參數(shù)的迭代過程。可以看出,目標函數(shù)值下降較快,本文模型(9)的目標函數(shù)值和待求模型參量均逼近經(jīng)典模型(8)的目標函數(shù)值和待求參量值。 試驗表明,直接反求法和數(shù)據(jù)壓縮反求法的迭代過程呈現(xiàn)較大差異。圖12給出了隨數(shù)據(jù)量增大,直接反求法與數(shù)據(jù)壓縮反求法的迭代次數(shù)對比??梢钥闯?,隨數(shù)據(jù)量的增大,本文方法的收斂迭代次數(shù)較直接反求法的迭代次數(shù)明顯減少。 圖9 復(fù)合材料三點彎曲試樣 Fig.9 Composite three -point bending testspecimen 結(jié)果表明,在識別九個復(fù)合材料參數(shù)時,基于PCA數(shù)據(jù)壓縮的計算反求方法對多個未知材料參數(shù)在較大區(qū)間下具有較強的識別能力,且反求的精度和速度均較為理想。 在DIC測量中,圖像噪聲引起的隨機誤差對系統(tǒng)的正常工作有著很大的影響[19],有必要對方法的抗噪聲性能進行分析。為此,本文直接在數(shù)值模型的基礎(chǔ)上考察噪聲對所提方法的影響。依據(jù)Vic-3D的測量誤差在100 μm范圍內(nèi),在理論計算位移數(shù)據(jù)中添加均值為0、方差σ=100 μm的高斯隨機噪聲,分析直接反求法和基于數(shù)據(jù)PCA壓縮的模型參量反求算法在引入噪聲下的反求效果。 圖10 兩種方法的目標函數(shù)值對比(2667個原始數(shù)據(jù)) Fig.10 Comparison of objective function by two methods(2667 data) 圖11 4個材料參數(shù)迭代過程對比(2667個原始數(shù)據(jù)) Fig.11 Comparison of iterative process of four material parameters (2667 data) 試驗表明,對比無噪聲結(jié)果,在引入噪聲的情況下,基于數(shù)據(jù)PCA壓縮的模型參量反求算法的結(jié)果精度變化更小。同時,對比直接反求,在有無噪聲情況下,本算法的迭代次數(shù)均明顯較小,呈現(xiàn)出更快的收斂速度、更強的穩(wěn)定性和抑噪能力。 圖12 數(shù)據(jù)量對迭代次數(shù)的影響 Fig.12 Effects of data volume on iteration times 本文提出一種基于變形場測量數(shù)據(jù)PCA壓縮的模型參量反求算法,以解決DIC測量數(shù)據(jù)量龐大導(dǎo)致的數(shù)據(jù)使用困難。算法采用測量平面內(nèi)全場性和無序性的數(shù)據(jù)點,通過PCA壓縮數(shù)據(jù),選擇滿足累計方差解釋率的主元信息進行反求,確保DIC位移測量的重要特征不會因壓縮而丟失。 針對具體的試驗?zāi)P?,研究表明,基于?shù)據(jù)PCA壓縮后的模型參量反求適用于多個參數(shù)的反求,并具有較高精度和較好穩(wěn)定性。特別是與直接反求法相比較,本文方法具有更快的收斂速度和更高的收斂精度,在處理噪聲方面,該方法具有降噪和去冗余的能力,能有效降低噪聲的影響。2.3 試驗分析
3 基于數(shù)據(jù)PCA壓縮的模型參量反求方法
4 算 例
4.1 金屬材料拉伸試驗
4.2 復(fù)合材料三點彎曲試驗
4.3 方法的抗噪性
5 結(jié) 論