潘海珠，葛海淼，蘇小紅，袁 琪

（1.齊齊哈爾大學(xué)計(jì)算機(jī)與控制工程學(xué)院，黑龍江齊齊哈爾 161006；2.哈爾濱工業(yè)大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，哈爾濱150001）

0 引言

由功能梯度材料（Functionally Graded Material，F(xiàn)GM）制備的受力部件在循環(huán)應(yīng)力或應(yīng)變作用下會(huì)產(chǎn)生局部結(jié)構(gòu)變化，形成疲勞裂紋并產(chǎn)生擴(kuò)展，對(duì)構(gòu)件及工程安全構(gòu)成嚴(yán)重威脅。裂紋尖端的應(yīng)力場(chǎng)強(qiáng)度，即應(yīng)力強(qiáng)度因子（Stress Intensity Factor，SIF），是表征材料斷裂的重要參量，它和材料屬性、應(yīng)力大小和裂紋幾何尺寸等因素密切相關(guān)。由于FGM 屬性呈現(xiàn)梯度變化，使得傳統(tǒng)的理論解析求解方法較難求解裂紋問題。隨著信息技術(shù)不斷發(fā)展，許多學(xué)者嘗試運(yùn)用人工智能方法與傳統(tǒng)理論求解模型相結(jié)合，來解決復(fù)雜的斷裂問題［1-5］。其中，支持向量回歸（Support Vector Regression，SVR）是支持向量機(jī)（Support Vector Machine，SVM）在解決回歸和預(yù)測(cè)問題中的擴(kuò)展［6］，作為一種高效的機(jī)器學(xué)習(xí)方法，采用最小化經(jīng)驗(yàn)誤差和最大化幾何邊緣策略，常用來解決小樣本、非線性問題。但是，SVR 方法在內(nèi)核參數(shù)選擇上存在一定困難，易使模型陷入局部最優(yōu)，文獻(xiàn)［7-12］中提出了基于SVR的改進(jìn)算法進(jìn)行參數(shù)優(yōu)化，可提高預(yù)測(cè)模型的準(zhǔn)確性和魯棒性。改進(jìn)的SVR 預(yù)測(cè)模型研究可以作為傳統(tǒng)斷裂力學(xué)研究的新手段，用于快速得到力學(xué)本構(gòu)模型和控制方程的函數(shù)表達(dá)形式?；诟怕首儺惢依莾?yōu)化算法的SVR是支持向量回歸的一種改進(jìn)類型，通過變異的灰狼優(yōu)化算法對(duì)SVR 的幾個(gè)主要參數(shù)進(jìn)行優(yōu)化，實(shí)現(xiàn)小樣本情況下的SIFs預(yù)測(cè)。

本文以含有任意方向裂紋的一般屬性FGM 為研究對(duì)象，分析應(yīng)力強(qiáng)度因子與FGM屬性、裂紋長度、裂紋角度等的關(guān)系，提出了SVR結(jié)合概率變異的灰狼優(yōu)化算法（Probability Mutation Grey Wolf Optimization Algorithm，MP-GWO），通過優(yōu)化SVR 參數(shù)，運(yùn)用分層指數(shù)模型方法獲取數(shù)據(jù)，訓(xùn)練一個(gè)以裂紋幾何參數(shù)及FGM屬性等特征值為輸入，應(yīng)力強(qiáng)度因子為輸出的預(yù)測(cè)模型。在傳統(tǒng)力學(xué)建模方法獲取數(shù)據(jù)基礎(chǔ)上，提出一種新的基于機(jī)器學(xué)習(xí)方法的斷裂力學(xué)預(yù)測(cè)模型，并通過有效實(shí)驗(yàn)驗(yàn)證該模型的適用性和精確性。本文的目的在于運(yùn)用分層指數(shù)模型獲取離散應(yīng)力強(qiáng)度因子，訓(xùn)練一個(gè)支持向量回歸模型，提出了概率變異灰狼優(yōu)化（MP-GWO）算法對(duì)模型參數(shù)進(jìn)行優(yōu)化，將這些離散應(yīng)力強(qiáng)度因子在高維空間擬合成一個(gè)超平面，實(shí)現(xiàn)了運(yùn)用機(jī)器學(xué)習(xí)方法擬合出隱式的本構(gòu)關(guān)系，對(duì)不同材料參數(shù)及裂紋幾何參數(shù)條件下的SIFs預(yù)測(cè)；對(duì)比傳統(tǒng)的灰狼優(yōu)化算法，提出的改進(jìn)算法，即基于MP-GWO算法的支持向量回歸模型能更好抑制算法陷入局部最優(yōu)，得到的應(yīng)力強(qiáng)度因子預(yù)測(cè)模型泛化能力更優(yōu)。

1 MP-GWO-SVR原理

1.1 SVR基本原理

設(shè)樣本集｛（xi，yi），i=1，2，…，n｝，其中，xi∈RD為數(shù)據(jù)樣本的輸入值；yi∈R為樣本的輸出值；D為樣本維度。對(duì)于非線性的SVR問題，可以引入映射函數(shù)φ（·），將輸入數(shù)據(jù)映射到高維特征空間F中，定義回歸預(yù)測(cè)模型

式中：w 為權(quán)重向量；b為常數(shù)；f（x）為回歸預(yù)測(cè)模型輸出。引入松弛變量（，ξi），定義f（xi）與yi的損失函數(shù)

式中：C是懲罰因子；ε 是最大允許誤差。目標(biāo)函數(shù)為最小化損失函數(shù)J，引入拉格朗日乘數(shù)及核函數(shù)可以將目標(biāo)函數(shù)轉(zhuǎn)化為如下對(duì)偶形式：

式中：αi和為拉格朗日乘數(shù)；K（xi，xj）=φ（xi）φ（xj）為核函數(shù)，定義徑向基核函數(shù)

式中，σ為內(nèi)核參數(shù)［7］。通過求得的αi和，可以得到SVR表達(dá)式

式中：xr和xs是任意的支持向量；SVs 為支持向量的數(shù)量。

在SVR 應(yīng)用中，參數(shù)懲罰因子C、RBF 核函數(shù)參數(shù)σ 和最大允許誤差ε 對(duì)模型的學(xué)習(xí)能力和泛化能力有重要影響，因此對(duì)這些參數(shù)進(jìn)行優(yōu)化尤為重要，本文提出使用MP-GWO算法優(yōu)化SVR的3 個(gè)參數(shù)。

1.2 概率變異灰狼優(yōu)化（MP-GWO）算法

（1）灰狼優(yōu)化（Grey Wolf Optimization，GWO）算法。GWO算法［13］是一種模擬狼群狩獵行為的群體智能優(yōu)化算法。GWO 算法在搜索空間隨機(jī)產(chǎn)生灰狼種群，將狼群劃分為不同等級(jí)的α 狼、β 狼和γ 狼，剩余個(gè)體為ω狼；通過包圍獵物、攻擊獵物行為實(shí)現(xiàn)對(duì)目標(biāo)獵物的獵殺，獲得待優(yōu)化問題的最優(yōu)解。

GWO的包圍行為首先確定灰狼與獵物之間的距離，更新灰狼的位置，其數(shù)學(xué)表達(dá)式如下：

式中：t表示當(dāng)前的迭代次數(shù)；D表示灰狼與獵物的距離；XP（t）為第t代獵物的位置；X（t）為第t代單只灰狼的位置向量；A和C為隨機(jī)常系數(shù)，

r1和r2為［0，1］間的隨機(jī)數(shù)；數(shù)值ɑ隨迭代次數(shù)增加從2 遞減至0，表達(dá)式為

l為當(dāng)前迭代次數(shù)；nmax為最大迭代次數(shù)。

GWO的狩獵行為在α狼、β狼和γ 狼的領(lǐng)導(dǎo)下進(jìn)行，其數(shù)學(xué)表達(dá)式為：

式中：Xα、Xβ、Xγ分別為頭狼α 狼、β 狼和γ 狼位置；Dα、Dβ、Dγ分別為狼群剩余個(gè)體ω狼與α狼、β 狼和γ狼的距離；X1、X2、X3分別用于計(jì)算3 個(gè)頭狼的運(yùn)動(dòng)軌跡；X（t+1）為灰狼更新后的位置。GWO 算法通過迭代過程更新狼群位置，逐漸逼近獵物，完成對(duì)獵物的捕食。

GWO算法的全局搜索能力體現(xiàn)在隨機(jī)數(shù)A和C，但該算法在初始化群體時(shí)就傾向于搜索可能存在的最優(yōu)解，這種搜索策略易導(dǎo)致結(jié)果局部最優(yōu)［14-17］。為改進(jìn)GWO算法的全局搜索能力，并進(jìn)一步增強(qiáng)算法的魯棒性，本文提出一種變異策略的灰狼尋優(yōu)算法。

（2）MP-GWO 算法。MP-GWO 算法在生成狼群位置X（t+1）后，以遞減概率對(duì)狼群位置進(jìn)行變異計(jì)算。變異計(jì)算分為兩種，分別為交叉變異和普通變異。狼群中的個(gè)體根據(jù)變異概率產(chǎn)生某一種變異。

交叉變異根據(jù)當(dāng)前最優(yōu)解Xα和當(dāng)前狼的位置X（t+1）決定ω狼的位置，交叉變異具體表達(dá)式如下：

式中：r3，r5為［0，1］間的隨機(jī)數(shù)；Mp為變異概率閾值。

普通變異根據(jù)當(dāng)前狼位置X（t+1）和搜索空間范圍決定ω狼的位置，普通變異的表達(dá)式如下：

式中：r4為［0，1］間的隨機(jī)數(shù)；ub和lb分別為搜索空間上下界。式（19）、（20）限制條件中的隨機(jī)數(shù)r5用于控制交叉變異或普通變異是否發(fā)生。

為保證狼群個(gè)體只能產(chǎn)生一種變異，變異概率閾值為［0，0.5］之間的值。為保證算法收斂，變異概率閾值應(yīng)隨迭代次數(shù)增加而減小，變異概率閾值Mp表達(dá)式如下

式中：p為變異概率常數(shù)，為［0，0.5］間的值。

MP-GWO算法的具體步驟如下：

（1）初始化GWO算法的參數(shù)。狼群中狼的數(shù)量Nw，變異概率常數(shù)p，搜索空間上下限ub和lb，搜索空間維度D，最大迭代次數(shù)nmax。

（2）計(jì)算灰狼個(gè)體的適應(yīng)度，根據(jù)適應(yīng)度值的大小進(jìn)行排序，前3 名灰狼位置設(shè)置為Xα，Xβ，Xγ。

（3）根據(jù)式（9）～（11）、（21）更新A，C，ɑ，Mp的值。

（4）遍歷每個(gè)狼個(gè)體，根據(jù)式（12）～（18）更新狼的當(dāng)前位置；當(dāng)r5∈［0，Mp］，其中r5為［0，1］間隨機(jī)數(shù)，根據(jù)式（19）計(jì)算當(dāng)前狼位置的交叉變異值，更新位置；當(dāng)r5∈［MP，2MP］，根據(jù)式（20）計(jì)算當(dāng)前狼位置的普通變異值，更新位置。

（5）判斷每一個(gè)狼個(gè)體的位置是否超出搜索空間上下限ub和lb，如位置值大于ub，則將當(dāng)前狼位置值設(shè)為ub；如位置值小于lb，則將當(dāng)前狼位置值設(shè)為lb。

（6）判斷是否達(dá)到最大迭代次數(shù)，若達(dá)到則跳出算法，輸出最優(yōu)值Xα；否則，設(shè)置l=l+1，跳轉(zhuǎn)到步驟（2）。

2 應(yīng)力強(qiáng)度因子計(jì)算

圖1 為含有任意方向裂紋具有一般材料屬性功能梯度材料示意圖。材料無限長且厚度為l，所含裂紋與材料梯度變化方向的角度為α。如圖1 所示，在x′Oy′坐標(biāo)系下，裂紋兩個(gè)端點(diǎn)坐標(biāo)分別為（ɑ，0），（b，0）。由于任意屬性FGM的裂紋問題難于解析求解，采用分層指數(shù)模型來求解一般材料屬性FGM斷裂問題，假設(shè)FGM實(shí)際屬性可以用指數(shù)函數(shù)、線性函數(shù)、冪函數(shù)等一般連續(xù)性函數(shù)來描述［18］。充分利用指數(shù)型材料的控制偏微分方程易求解這一特點(diǎn)，通過傅里葉變換、微分算子等方法求解應(yīng)力場(chǎng)和位移場(chǎng)；結(jié)合裂紋問題的邊界條件和位移條件，求解奇異積分方程組，求得混合型應(yīng)力強(qiáng)度因子k值如下：

圖1 含任意方向裂紋的功能梯度材料示意圖

式中：κ、δ 為材料參數(shù)；n為分層指數(shù)模型劃分的層數(shù)；f1、f2為求解奇異方程組引入的可求解函數(shù)。運(yùn)用該方法求得的應(yīng)力強(qiáng)度因子用于訓(xùn)練SVR預(yù)測(cè)模型。

3 MP-GWO-SVR應(yīng)力強(qiáng)度因子預(yù)測(cè)模型

由于SVR的性能主要受懲罰因子C、RBF核函數(shù)參數(shù)σ和最大允許誤差ε影響，使用MP-GWO算法實(shí)現(xiàn)對(duì)SVR 3 個(gè)參數(shù)進(jìn)行優(yōu)化。生成的SIFs 預(yù)測(cè)模型計(jì)算步驟如下：

（1）樣本數(shù)據(jù)初始化。材料屬性變化形式、模量比、裂紋相對(duì)長度和裂紋角度為輸入值，裂紋尖端應(yīng)力強(qiáng)度因子為輸出值，組成樣本集和標(biāo)簽集（Sample_all，label_all）?？紤]3 種材料屬性變化形式，即材料屬性變化分別用指數(shù)函數(shù)、線性函數(shù)和冪函數(shù)來描述；材料上下表面彈性模量比取值為（2，2.5，3，4，5）；裂紋相對(duì)長度的取值為（0.1，0.15，0.2，0.25，0.3，0.35，0.4，0.45，0.5）；裂紋角度的取值為（0π，0.05π，0.1π，0.15π，0.2π，0.25π，0.3π，0.35π，0.4π，0.45π，0.5π）；樣本總數(shù)為202 個(gè)；隨機(jī)選擇1/2 樣本作為訓(xùn)練樣本（Sample_train），共101 個(gè)樣本，其標(biāo)簽為label_train，其他樣本作為測(cè)試樣本（Sample_test），共101 個(gè)樣本。

（2）參數(shù)初始化。初始化狼群數(shù)量Nw=10，其中頭狼α、β和γ占狼群數(shù)量為3，剩余個(gè)體ω 狼占狼群數(shù)量為7；變異概率常數(shù)p=0.5；搜索空間上下限ub=［100，100，0.2］和lb=［0.01，0.01，0.001］；搜索空間維度D=3；nmax=10。

（3）使用Sample_train 訓(xùn)練SVR，使用MP-GWO算法尋找SVR的最佳參數(shù)C、σ和ε。

（4）使用C、σ 和ε 實(shí)現(xiàn)SVR 的建模，得到SVR模型model_SVR。

（5）使用model_SVR預(yù)測(cè)Sample_test的值，得到測(cè)試樣本的預(yù)測(cè)值label_test_SVR。

4 實(shí)驗(yàn)結(jié)果及分析

以冪函數(shù)表征材料屬性變化形式，材料上下表面彈性模量比值為2.5 的特征數(shù)據(jù)為例，將SVR模型預(yù)測(cè)結(jié)果，即無量綱化應(yīng)力強(qiáng)度因子kI（ɑ）、kI（b）、kII（ɑ）和kII（b）分別標(biāo)注為Label1、Label2、Label3和Label4，圖2～5 分別給出不同裂紋參數(shù)下MP-GWOSVR模型與GWO-SVR模型的SIFs預(yù)測(cè)值，以及運(yùn)用分層法求解SIFs（真實(shí)值）的變化曲線。

圖2 基于GWO-SVR的SIFs預(yù)測(cè)模型預(yù)測(cè)結(jié)果

對(duì)比圖2（a）和圖4（a）可以看到，對(duì)于不同長度的裂紋，當(dāng)裂紋角度較大時(shí)（如0.4π，0.45π，0.5π），應(yīng)力強(qiáng)度因子kI（ɑ）較真實(shí)值有較大的差別。比較圖3（a）和圖4（a），可以看到MP-GWO-SVR 模型預(yù)測(cè)SIFs曲線較分層模型求解SIFs 曲線變化趨勢(shì)相似度都很高；尤其注意到當(dāng)裂紋角度較大時(shí)（如0.4π，0.45π，0.5π），兩種模型給出曲線的變化趨勢(shì)幾乎一致。這種現(xiàn)象說明MP-GWO-SVR 模型較GWO-SVR模型在Label1上有更好的擬合度。對(duì)比圖2（b）、圖3（b）和圖4（b），可以看到MP-GWO-SVR 模型得到了與分層模型求解SIFs曲線相似度很高的擬合曲線，而GWO-SVR模型對(duì)于裂紋角度較大時(shí)預(yù)測(cè)結(jié)果與真實(shí)值偏離較大，這種現(xiàn)象說明MP-GWO-SVR 模型在Label2上可得到較為理想的擬合結(jié)果。對(duì)比圖2（c）、圖3（c）和圖4（c），可以看到情況有所不同。在裂紋相對(duì)長度為0.1 和0.2 時(shí)，GWO-SVR 和MP-GWOSVR兩種模型預(yù)測(cè)SIFs曲線都有變形，但是從圖中可以看到，MP-GWO-SVR模型的預(yù)測(cè)曲線變形程度明顯小于GWO-SVR 模型；當(dāng)裂紋相對(duì)長度為0.3 時(shí)，GWO-SVR模型和MP-GWO-SVR 模型給出的SIFs 預(yù)測(cè)曲線擬合度較高；當(dāng)裂紋相對(duì)長度為0.4 和0.5 時(shí)，GWO-SVR模型的SIFs預(yù)測(cè)值隨裂紋角度變化出現(xiàn)明顯波動(dòng)，相比真實(shí)值曲線有明顯的變形，而MP-GWOSVR模型的SIFs預(yù)測(cè)曲線與分層模型求解SIFs 曲線形狀基本一致。這種情況說明對(duì)Label3，GWO-SVR和MP-GWO-SVR模型對(duì)SIFs 擬合結(jié)果都有的誤差，但MP-GWO-SVR模型的預(yù)測(cè)誤差相對(duì)較小。對(duì)比圖2（d）、圖3（d）和圖4（d），可以看到GWO-SVR 模型給出的SIFs擬合曲線較分層模型求解SIFs 曲線有較明顯的變形，而MP-GWO-SVR模型給出的擬合曲線較分層模型求解SIFs曲線擬合度很高，說明MP-GWO-SVR模型在Label4上擬合能力優(yōu)于GWO-SVR模型。以上分析結(jié)果與圖5 中的SIFs樣本點(diǎn)對(duì)比結(jié)果一致。

圖3 基于MP-GWO-SVR的SIFs預(yù)測(cè)模型預(yù)測(cè)結(jié)果

圖4 分層模型求解SIFs

圖5 GWO-SVR、MP-GWO-SVR預(yù)測(cè)值與真實(shí)值作對(duì)比

為了進(jìn)一步比較MP-GWO算法和GWO算法的性能，本文采用均方誤差（Mean Square Error，MSE）和誤差范圍百分比（Percentage of Error Range，PER）兩種度量標(biāo)準(zhǔn)分別對(duì)基于兩種算法的SVR 生成SIFs 預(yù)測(cè)模型進(jìn)行預(yù)測(cè)精度評(píng)價(jià)，其定義如下：

式中：n為樣本數(shù)；f（xi）為樣本的預(yù)測(cè)值；yi為樣本的真實(shí)值（即運(yùn)用分層模型得到的解）；f（X）為預(yù)測(cè)值向量；Y為真實(shí)值向量；max｛·｝為求目標(biāo)集合的最大值；min｛·｝為求目標(biāo)集合的最小值。

選取其中兩組數(shù)據(jù)為例，分別計(jì)算其均方誤差及誤差范圍百分比，結(jié)果如表1、2 所示。由表1 可見，基于算法GWO 的SVR 模型在Label1、Label2、Label3和Label4上的均方誤差在4.583 5 ×10-4～3.939 1 ×10-5之間，其中Label1和Label2的均方誤差較小，Label3和Label4的均方誤差較大；基于算法MP-GWO的SVR模型在各標(biāo)簽上的均方誤差范圍為7.556 2 ×10-5～5.364 4 ×10-6，其中Label1的均方誤差較小，Label2、Label3和Label4的均方誤差較大。對(duì)比算法GWO和MP-GWO，可以看到MP-GWO-SVR 模型的均方誤差在各個(gè)Label 上都小于GWO-SVR 模型。這種現(xiàn)象說明MP-GWO 算法構(gòu)建的預(yù)測(cè)模型有較高的泛化能力。從表2 可以看到，基于算法MP-GWO的SVR模型在Label1、Label2和Label3上的誤差范圍都小于GWO-SVR模型，兩者在Label4上有相似的誤差范圍。這種現(xiàn)象表明，MP-GWO 算法較GWO 算法相比，其SVR模型不易出現(xiàn)偏離實(shí)際值過大的預(yù)測(cè)值，有更好的穩(wěn)定性。

表1 兩種算法構(gòu)建預(yù)測(cè)模型的均方誤差對(duì)比

表2 兩種算法構(gòu)建預(yù)測(cè)模型的誤差范圍對(duì)比 %

5 結(jié)語

本文以含任意方向裂紋的功能梯度材料斷裂問題為研究對(duì)象，運(yùn)用機(jī)器學(xué)習(xí)方法，用小樣本數(shù)據(jù)訓(xùn)練一個(gè)支持向量回歸模型，實(shí)現(xiàn)對(duì)含有任意方向裂紋的一般屬性功能梯度材料應(yīng)力強(qiáng)度因子的預(yù)測(cè)分析。本文在傳統(tǒng)的基于灰狼優(yōu)化算法的SVR 模型基礎(chǔ)上提出了一種改進(jìn)算法，即概率變異灰狼優(yōu)化算法，構(gòu)建SVR模型預(yù)測(cè)SIFs。該算法優(yōu)化SVR參數(shù)，獲取模型的預(yù)測(cè)精度比基于傳統(tǒng)灰狼優(yōu)化算法支持的預(yù)測(cè)模型有顯著提高，從得到的誤差范圍來看，概率變異灰狼優(yōu)化算法SVR預(yù)測(cè)模型的預(yù)測(cè)結(jié)果與樣本數(shù)據(jù)的擬合度更好。本文在傳統(tǒng)斷裂力學(xué)研究方法基礎(chǔ)上運(yùn)用智能優(yōu)化算法，為斷裂力學(xué)問題研究提供了一種可擴(kuò)展的研究模式，該方法將在今后的斷裂力學(xué)理論與實(shí)驗(yàn)研究中將發(fā)揮重大作用。

·名人名言·

對(duì)真理的追求要比對(duì)真理的占有更可貴。

——愛因斯坦