余 波, 陶盈盈

(1. 合肥工業(yè)大學(xué) 土木與水利工程學(xué)院 工程力學(xué)系, 合肥 230009;2. 大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室, 遼寧 大連 116024)

0 引 言

天然氣和石油是當(dāng)今世界上重要的化石能源,管道運輸是其主要的運輸方式．由于長時間的運輸,管道內(nèi)部很容易發(fā)生腐蝕老化,進(jìn)而導(dǎo)致管壁破裂和油氣泄漏等災(zāi)難性事件的發(fā)生[1]．因此,管道缺陷檢測對保證其運輸安全尤為重要,識別管道內(nèi)壁被腐蝕的幾何形狀對精準(zhǔn)評估管道的腐蝕程度具有深刻的指導(dǎo)意義．

幾何形狀識別[2-4]主要通過系統(tǒng)外表面響應(yīng)數(shù)據(jù)進(jìn)行內(nèi)部幾何形狀估算,屬于反幾何問題．近年來,多種無損檢測技術(shù)[5-6]已被廣泛應(yīng)用于管道狀況的評估,其中漏磁檢測法[7]由于檢測范圍廣、靈敏度高、適應(yīng)性強以及零污染等優(yōu)勢備受關(guān)注．顯然,準(zhǔn)確快速地求解磁場正問題是確保幾何形狀高效識別的前提．在眾多正問題數(shù)值求解方法中,有限元法因其強大的適應(yīng)性和穩(wěn)定性,迄今為止已被廣泛地應(yīng)用于求解磁場問題[8-10],并獲得了準(zhǔn)確的結(jié)果．因此,本文采用有限元法求解磁場的相關(guān)響應(yīng)量,并以此為基礎(chǔ)進(jìn)行管道內(nèi)壁幾何形狀的識別．

選用合適的優(yōu)化方法是反幾何問題成功的另一關(guān)鍵點,優(yōu)化方法可分為局部搜索的梯度類優(yōu)化算法和全局搜索的演化類優(yōu)化算法．由于在幾何形狀識別過程中其形狀會隨迭代過程不斷更新,從而對目標(biāo)函數(shù)靈敏度的求解帶來一定的困難,且靈敏度的求解在一定程度上會增加反演計算的成本．為避免靈敏度的計算并實現(xiàn)全局搜索的目標(biāo),本文選用非梯度類算法進(jìn)行目標(biāo)函數(shù)的優(yōu)化．具有代表性的方法包括遺傳算法[11]、布谷鳥搜索算法[12]、蟻群算法[13-14]、粒子群算法[15-16]和灰狼優(yōu)化算法(grey wolf optimization, GWO)[17-18]等,其中GWO算法因其具有簡單高效的優(yōu)勢而廣受學(xué)者的關(guān)注．Mirjalili等[17]提出了GWO算法并將其應(yīng)用于解決工程設(shè)計問題．Kohli等[19]將混沌理論引入到GWO算法中,以加快其全局收斂速度．文獻(xiàn)[20]對目前GWO算法的相關(guān)改進(jìn)及應(yīng)用進(jìn)行了總結(jié),驗證了該方法的可行性和有效性．基于此,本文選取GWO算法來實現(xiàn)管道內(nèi)壁的幾何形狀識別．

在反演迭代過程中需多次實施正問題分析,尤其對于全局類演化搜索算法,調(diào)用正問題分析的次數(shù)或許可達(dá)成千上萬次．此外,對于大型復(fù)雜模型,因其高維數(shù)和復(fù)雜性,直接進(jìn)行分析相對困難且數(shù)值模擬耗時過長[21]．為此,對數(shù)值模型的規(guī)?；螂A數(shù)進(jìn)行有效的降階處理就顯得尤為重要．Hellstr?m等[22]利用本征正交分解法(proper orthogonal decomposition, POD)進(jìn)行了管道流動中結(jié)構(gòu)的識別．Eftekhar Azam等[23]結(jié)合該方法和人工神經(jīng)網(wǎng)絡(luò),提出了一種監(jiān)督學(xué)習(xí)方法,并用于檢測、定位和量化結(jié)構(gòu)損傷強度．POD是預(yù)測響應(yīng)最受歡迎的模型降階方法之一,該方法以有限的數(shù)據(jù)捕獲計算過程中的主要成分,可顯著降低自由度,提升計算效率．但該方法的精度受樣本數(shù)量與樣本間相關(guān)性等因素的影響,樣本的選擇至關(guān)重要．利用本文建立的變幾何樣本庫可有效處理這一難題．

常規(guī)POD的預(yù)測模型隨著幾何形狀的改變需要反復(fù)更新有限元剛度矩陣,這在一定程度上嚴(yán)重增加了POD模型的計算成本．因此,探索一種代理模型以避免在計算過程中因幾何改變反復(fù)求解有限元剛度矩陣是非常有必要的．常見的代理模型有多項式回歸(polynomial regression, PR)[24-25]、Kriging(KRG)[26-27]、徑向基函數(shù)(radial basis function, RBF)[28-30]和人工神經(jīng)網(wǎng)絡(luò)(artificial neural network, ANN)[31-32]等．其中,PR模型是通過樣本數(shù)據(jù)基于最小二乘原理來擬合近似多項式．由于PR的模型簡單、計算量小,且多項式的平滑能力能使帶有噪聲的函數(shù)快速收斂,因而成為最常用的代理模型之一．然而,PR在處理高度非線性問題時,由于高次多項式的使用會出現(xiàn)不穩(wěn)定現(xiàn)象．KRG是對區(qū)域變量求無偏內(nèi)插估計值的一種插值方法,預(yù)測精度主要依賴于初始采樣,可能會導(dǎo)致模型過早停止或過于局部收斂．ANN是一種模擬人腦分析和處理信息方式的人工智能算法．該方法從復(fù)雜的數(shù)據(jù)中學(xué)習(xí),確定輸入和輸出變量之間的關(guān)系,從而進(jìn)行預(yù)測分析．RBF模型以關(guān)于Euclidean距離或者其他類似度量的函數(shù)為基函數(shù)[33-34],通過線性加權(quán)來插值擬合數(shù)據(jù)．Jin等[35]利用多種性能評價標(biāo)準(zhǔn)和14個標(biāo)準(zhǔn)測試函數(shù)系統(tǒng)地比較研究了PR、RBF和KRG等代理模型的性能,得出在處理不同階數(shù)的非線性和問題規(guī)模時,RBF模型在準(zhǔn)確性和穩(wěn)定性等方面的表現(xiàn)最好．Jing等[36]提出了一種基于自適應(yīng)RBF與遺傳算法相結(jié)合的可靠性評估方法,以降低結(jié)構(gòu)可靠性分析的計算成本．Liu等[37]針對多目標(biāo)優(yōu)化計算量大的問題,利用RBF提出一種基于自適應(yīng)逼近模型的高效多目標(biāo)優(yōu)化方法,并通過算例驗證了該方法的有效性和實用性．綜上,若將POD和RBF的優(yōu)勢耦合,即可實現(xiàn)在降階分析正問題的同時避免剛度矩陣的重復(fù)計算．例如,Khatir等[38]基于擴(kuò)展等幾何分析,將斷裂力學(xué)試驗與數(shù)值模型相結(jié)合,利用POD-RBF識別了板結(jié)構(gòu)的單、多裂紋．Henneron等[39]應(yīng)用該方法能夠在合理的計算時間和良好的精度下得到有限元的近似解,從而模擬非線性靜磁器件．Wang等[40]基于仿真結(jié)果,采用POD-RBF的降階方法,預(yù)測了不同滑移幅值和不同微動磨損循環(huán)次數(shù)下試樣的磨損特性．研究表明,POD-RBF可高效準(zhǔn)確地預(yù)測相關(guān)響應(yīng)量,顯著地降低計算成本．截至目前,已出版的研究工作基本上是基于固定幾何樣本庫采用POD-RBF進(jìn)行響應(yīng)量的預(yù)測．為此,本文基于構(gòu)建的變幾何樣本庫建立了POD-RBF降階代理模型．該方法可期望在識別過程中避免因管道內(nèi)壁幾何形狀的改變而需反復(fù)求解剛度矩陣,在滿足擬合精度要求的前提下降低計算成本．

1 管道內(nèi)部幾何反演理論

1.1 靜磁場二維有限元基本理論

本文在靜磁場環(huán)境中進(jìn)行管道內(nèi)壁幾何形狀的識別,將其簡化為平面問題,對應(yīng)的平衡方程為[41]

(1)

其中x=(x,y),μr為相對磁導(dǎo)率,μ0為空氣磁導(dǎo)率,Az為磁勢,Jz為電流密度．

考慮兩類邊界條件:

(2)

如圖1所示,Γ=Γ1∪Γ2表示域Ω的邊界,m,n分別是Γ2和Γd的法向單位矢量,Bt為切向磁通密度．

圖1 求解域示意圖Fig. 1 The diagram of the solution domain

在磁導(dǎo)率不同的兩種媒質(zhì)界面Γd上,磁勢應(yīng)滿足連續(xù)性條件:

(3)

對式(1)應(yīng)用Galerkin有限元法可得

(4)

通過有限元離散可得

KAz=Jz,

(5)

其中K為總體剛度矩陣,Az為結(jié)點磁勢向量,Jz為電流載荷向量．

不同管道內(nèi)壁幾何對應(yīng)模型不同的表面磁勢信息,其中管道內(nèi)壁幾何由若干個幾何參數(shù)表示．在管道內(nèi)壁識別過程中,通過有限元計算獲得不同管道內(nèi)壁幾何相應(yīng)的磁勢信息,進(jìn)而生成樣本,并建立磁勢數(shù)據(jù)與幾何參數(shù)的關(guān)系．另外,利用管道內(nèi)壁真實幾何模型的磁勢信息,對管道內(nèi)壁幾何形狀進(jìn)行預(yù)測．

1.2 本征正交分解-徑向基函數(shù)

在采用POD-RBF求解之前,建立有效的樣本庫關(guān)系到問題的求解精度,因此選取合適的樣本點至關(guān)重要．拉丁超立方抽樣(Latin hypercube sampling, LHS)方法是目前使用較為廣泛的抽樣方法之一,該抽樣方法利用分層的思想,保證樣本點是從給定的設(shè)計空間內(nèi)均勻隨機抽取的．基于LHS生成POD-RBF降階代理模型所需樣本庫的參數(shù),有效避免了計算過程中因管道內(nèi)壁幾何形狀改變反復(fù)更新剛度矩陣,能在保證計算精度的同時提高計算效率．

絕大部分工作均選用某個方向上的磁感應(yīng)強度進(jìn)行研究,我們期望從數(shù)值上建立全面反應(yīng)有關(guān)磁感應(yīng)強度信息的數(shù)據(jù)與模型幾何的關(guān)系,故選用有限元分析中的磁勢作為響應(yīng)量．在相同邊界條件下,通過LHS生成對應(yīng)的不同管道內(nèi)壁幾何參數(shù),進(jìn)而利用有限元計算得到磁勢樣本矩陣:

ψ=[ζ(α1),ζ(α2),ζ(α3),…,ζ(αN)],

(6)

式中αi為控制管道內(nèi)壁幾何的參數(shù),ζ(αi)為n維磁勢列向量,N表示樣本個數(shù)．

對樣本矩陣ψ進(jìn)行奇異值分解:

ψ=CSVT,

(7)

式中C∈Rn×n為ψ的左奇異向量組成的矩陣;V∈RN×N為ψ的右奇異向量組成的矩陣;S∈Rn×N僅在主對角線上有值,稱為ψ的奇異值,可表示為si(i=1,2,…,min(n,N))．所有奇異值按降序排列,通過以下規(guī)則

(8)

實現(xiàn)k個奇異值和對應(yīng)奇異向量的自適應(yīng)截斷進(jìn)而可近似描述矩陣ψ,即

(9)

本文χ取99.999%．

(10)

(11)

其中i=1,2,…,N;q=1,2,…,k．

寫成矩陣形式為

σ=WΦ,

(12)

其中W為wqj組成的矩陣,Φ可表示為

(13)

為了提升計算效率這里選用inverse multiquadtric(InvM)核函數(shù)對Φ進(jìn)行近似展開,則Φ中任一元素φj(αi)可表示為

(14)

式中,c是取值大于零的平滑系數(shù),rji表達(dá)式如下:

rji=‖αi-αj‖2,

(15)

‖αi‖2為αi的2范數(shù)．由式(12)可得

W=σΦ-1．

(16)

通過式(10)、(12)、(16)可得任意參數(shù)αp對應(yīng)的磁勢響應(yīng)量:

(17)

由式(17)可以看出,通過變幾何樣本庫的構(gòu)建,在求解不同管道內(nèi)壁幾何對應(yīng)磁勢響應(yīng)量時,只需更新幾何對應(yīng)的基函數(shù)向量Φ(αp),這將在很大程度上節(jié)省計算成本．

1.3 目標(biāo)函數(shù)的優(yōu)化

通過確定管道內(nèi)壁幾何參數(shù)αp即可確定管壁的幾何形狀,建立與αp有關(guān)的目標(biāo)函數(shù)

(18)

受灰狼群體捕食行為的啟發(fā),GWO算法被提出．灰狼被認(rèn)為是食物鏈頂端的捕食者,大部分喜歡群居,一個群體平均有5～12只狼．在群體中,他們有著嚴(yán)格的等級制度,如圖2所示,從高到底依次是首領(lǐng)狼a,副首領(lǐng)狼b,普通狼c和底層狼ω．

圖2 灰狼的等級制度Fig. 2 The hierarchy of gray wolves

灰狼的狩獵是在首領(lǐng)狼a的帶領(lǐng)下,進(jìn)行圍攻捕食．捕食過程分為三個階段:

1) 通過氣味跟蹤、追逐、接近獵物;

2) 鎖定獵物位置后,進(jìn)行包圍;

3) 快速地攻擊獵物．

在狼群捕食過程中,圍剿獵物公式為

C=|CXp(t)-X(t)|,

(19)

C=2?1,

(20)

式中Xp(t)為獵物的位置,t是迭代步數(shù),X(t)為灰狼的位置,?1從[0,1]范圍內(nèi)隨機取值．灰狼的位置更新如下:

X(t+1)=Xp(t)-BD,

(21)

在更新過程中,保留當(dāng)前結(jié)果最好的前三組解,依次為a,b和c狼的位置．狼群中其他狼的位置根據(jù)a,b和c狼的位置進(jìn)行更新:

(22)

(23)

通過式(22)、(23)不斷更新前三個最優(yōu)灰狼的位置,直到滿足優(yōu)化標(biāo)準(zhǔn)條件．

2 數(shù) 值 算 例

2.1 基于POD-RBF對響應(yīng)量的預(yù)測

取管道任一橫截面的附近區(qū)域進(jìn)行分析,假定待識別的管道內(nèi)壁幾何分別為圓、橢圓和不規(guī)則形狀．圖3是該問題的簡化模型,幾何尺寸如圖所示,其中方塊被選為參考點(即后續(xù)算例中的測點),區(qū)域Ι、Ⅱ和Ⅲ分別為石油(或天然氣)、管道壁和空氣．在正方形上邊界(y=1)處施加磁勢Az=0.5 Wb/m,下邊界(y=-1)處施加Az=0 Wb/m．天然氣和石油由于抗磁性它們的相對磁導(dǎo)率和空氣相近,即相對磁導(dǎo)率μr均取為1．管道材料選用常用的X52鋼,對應(yīng)的B-H取值參考文獻(xiàn)[42]．

圖3 圓型管道內(nèi)壁Fig. 3 The inner wall of the circular pipeline

首先對POD-RBF模型的正確性進(jìn)行驗證,假定管道內(nèi)壁形狀為圓,該幾何對應(yīng)的參數(shù)αp里僅有一個元素,其值為0.533 m．樣本對應(yīng)的參數(shù)α通過LHS在0.5～0.6 m生成10組．用POD-RBF模型預(yù)測被對比點的磁勢．有限元計算采用9 856個四節(jié)點四邊形單元,10 037個節(jié)點,單元劃分如圖4所示,圖5展示了相應(yīng)磁場強度的大小．兩種方法計算結(jié)果誤差如圖6所示．

圖6(b)顯示POD-RBF計算預(yù)測的響應(yīng)量與有限元計算結(jié)果的最大相對誤差不超過0.06%,即該算法對求解管道內(nèi)部內(nèi)壁幾何正問題有效．

圖4 有限元網(wǎng)格 圖5 磁場強度大小Fig. 4 The FEM meshFig. 5 The magnetic field intensity

(a) 參考點磁勢的絕對誤差 (b) 參考點磁勢的相對誤差(a) Absolute errors of the magnetic potential at the reference points(b) Relative errors of the magnetic potential at the reference points圖6 參考點磁勢的誤差Fig. 6 The magnetic potential errors at the reference points

2.2 不同樣本數(shù)量對識別結(jié)果的影響

本小節(jié)對如圖3所示管道內(nèi)壁幾何為圓的問題進(jìn)行反演．利用LHS在0.5～0.6 m分別生成10組、20組和40組相應(yīng)的幾何參數(shù)α,再通過有限元分別建立N=10,20,40的樣本庫,進(jìn)一步探究不同樣本數(shù)量對反演結(jié)果的影響．

ε和tmax分別設(shè)定為10-10和100．從表 1對比分析的三組樣本方案結(jié)果顯示,不同數(shù)量的樣本均能達(dá)到較高的識別精度．如圖7和圖8所示,當(dāng)管道內(nèi)壁形狀為圓時,采用不同數(shù)量的樣本,識別結(jié)果均在16步內(nèi)收斂．使用40組樣本時,經(jīng)過6次迭代即可獲得更準(zhǔn)確的識別結(jié)果,但計算成本可能會隨著樣本數(shù)量的增加而增大．因此在復(fù)雜模型識別中,選擇合適數(shù)量的樣本至關(guān)重要．

圖7 不同樣本數(shù)量對應(yīng)的目標(biāo)函數(shù)值 圖8 不同樣本數(shù)量對應(yīng)的幾何參數(shù) Fig. 7 Objective function values with different sample sizes Fig. 8 Geometric parameters with different sample sizes

表1 識別結(jié)果

為了確定后續(xù)算例的樣本,采用繼承LHS[43]生成樣本,期望在保證精度的同時盡可能降低計算成本．經(jīng)測試,后續(xù)分別擬采用50組和40組樣本進(jìn)行橢圓形和不規(guī)則形內(nèi)壁幾何的識別．

同時樣本中對應(yīng)的幾何參數(shù)范圍也會影響識別結(jié)果,接下來我們討論在當(dāng)前樣本N=10下管壁幾何的識別情況．待識別的真實管壁幾何半徑在0.45～0.69 m內(nèi),識別結(jié)果如表 2所示．

表2 不同方案下的識別結(jié)果

從表中可知,真實內(nèi)壁幾何在0.50～0.61 m間識別結(jié)果的相對誤差均在1%內(nèi),在這范圍外的相對誤差隨著半徑的增大(或減小)而增加,當(dāng)半徑減小到0.45 m時,識別結(jié)果的相對誤差達(dá)到了10.892 1%．綜上,在有限的樣本下,待識別的幾何參數(shù)在所建立樣本庫相應(yīng)的幾何參數(shù)范圍內(nèi)識別結(jié)果更準(zhǔn)確．但當(dāng)樣本庫足夠大時,理論上可以識別出任意管壁的幾何．

2.3 基函數(shù)對反演結(jié)果的影響

如圖9所示,本小節(jié)假定管道內(nèi)壁形狀為橢圓,真實幾何對應(yīng)的參數(shù)αp=[0.625,0.42]Tm．通過InvM核函數(shù)和緊支撐的四階樣條函數(shù)[44]討論不同基函數(shù)對識別結(jié)果的影響．其中緊支撐的四階樣條函數(shù)可以表示為

圖9 橢圓型管道內(nèi)壁Fig. 9 The inner wall of the elliptical pipeline

(24)

樣本對應(yīng)的參數(shù)α由LHS在0.40～0.65 m生成50組．ε和tmax分別設(shè)定為10-8和300,兩種方法的迭代初始值均設(shè)置在0.9～1.0 m之間．

如圖10所示,使用緊支撐的四階樣條函數(shù)迭代總次數(shù)少于InvM核函數(shù),且兩種方法均能得到比較準(zhǔn)確的幾何．但在同樣的步數(shù)下,InvM核函數(shù)耗時更短．同時為兼顧反演成本,在后續(xù)算例中選用InvM核函數(shù)．

(a) Inverse multiquadtric核函數(shù)(a) Inverse multiquadtric kernel functions

2.4 不同誤差對反演結(jié)果的影響

為驗證本文方法的抗噪性,本小節(jié)在測點響應(yīng)上添加隨機誤差．算例采用如圖11所示的模型,管道內(nèi)壁的幾何形狀由12個參數(shù)控制,該組參數(shù)在平面內(nèi)均勻分布,從0°～360°每隔θ=30°選取,最后通過樣條曲線擬合成封閉圖形．樣本對應(yīng)的參數(shù)α由LHS在0.5～0.6 m內(nèi)生成40組, 真實幾何對應(yīng)的參數(shù)αp=[0.526,0.578,0.533,0.556,0.518,0.544,0.538,0.575,0.510,0.565,0.530,0.594]Tm．

圖11 不規(guī)則幾何形狀的管道內(nèi)壁Fig. 11 The inner wall of the pipeline with irregular geometry

考慮1%,2%和3%的隨機誤差,添加誤差后的測點磁勢可表示為

(25)

本算例假定ε和tmax分別為10-7和8 000．需要說明的是,有限元模型與POD-RBF降階代理模型的計算結(jié)果是存在一定偏差的．

表 3列出了在不同誤差水平下的參數(shù)識別結(jié)果．

圖12直接展現(xiàn)了本文方法識別帶有12個幾何參數(shù)問題的性能．值得注意的是,在δ=0時,識別出的幾何參數(shù)與真實參數(shù)存在微小差異,但擬合出的邊界與真實邊界幾乎完全重合．從圖13和圖14可以看到,在誤差水平分別為1%和2%時,識別結(jié)果與真實幾何仍吻合較好．甚至δ=3%時,圖15顯示本文方法仍能識別出幾何的基本輪廓．另外,從圖16更能直觀地觀察到,在不同的誤差水平下,本文方法可獲得較為理想的磁勢響應(yīng)場和識別結(jié)果．

圖12 δ=0時的識別結(jié)果 圖13 δ=1%時的識別結(jié)果 Fig. 12 The identified results with δ=0Fig. 13 The identified results with δ=1%

圖14 δ=2%時的識別結(jié)果 圖15 δ=3%時的識別結(jié)果 Fig. 14 The identified results with δ=2%Fig. 15 The identified results with δ=3%

(a) 真實幾何 (a) The true geometry

3 結(jié) 論

本文建立了變幾何樣本庫的POD-RBF降階代理模型,結(jié)合GWO算法構(gòu)建了一種新型的管道內(nèi)壁幾何識別框架．通過LHS生成樣本庫后,利用POD-RBF對響應(yīng)量進(jìn)行預(yù)測,可避免迭代過程中因幾何改變而反復(fù)更新剛度矩陣．GWO算法的引入實現(xiàn)了無需靈敏度計算的全局目標(biāo)函數(shù)優(yōu)化．?dāng)?shù)值算例表明,本文算法不僅能夠準(zhǔn)確預(yù)測正問題響應(yīng),同時可實現(xiàn)對管道內(nèi)壁幾何的準(zhǔn)確識別．即使對于帶有高維幾何參數(shù)且考慮噪聲問題,該算法仍具備強勁的識別性能,顯示出了良好的識別穩(wěn)定性．同時,本文的研究工作可為拓?fù)鋬?yōu)化、幾何識別和無損檢測等變幾何計算問題提供新思路．

致謝本文作者衷心感謝工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室開放基金(GZ21109)對本文的資助．