基于AR-GP模型的結(jié)構(gòu)損傷識(shí)別方法

唐啟智，辛景舟，周建庭，付 雷，周濱楓

(1.重慶交通大學(xué) 省部共建山區(qū)橋梁及隧道工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，重慶 400074；2.貴州橋梁建設(shè)集團(tuán)有限責(zé)任公司，貴陽(yáng) 550000；3.重慶大學(xué) 土木工程學(xué)院，重慶 400074)

隨著服役年限的增加，工程結(jié)構(gòu)物在運(yùn)營(yíng)過(guò)程中由于荷載效應(yīng)、材料劣化、環(huán)境侵蝕等因素的作用，其使用性能不斷惡化。識(shí)別結(jié)構(gòu)的早期損傷，掌握其健康狀況，有著重要的科學(xué)研究意義與工程實(shí)用價(jià)值[1]。

結(jié)構(gòu)損傷將導(dǎo)致結(jié)構(gòu)固有特性的改變，從而引起結(jié)構(gòu)響應(yīng)的變化。結(jié)構(gòu)損傷識(shí)別從結(jié)構(gòu)響應(yīng)出發(fā)，以期獲得結(jié)構(gòu)的損傷狀態(tài)，其核心環(huán)節(jié)即是損傷敏感性特征(damage sensitive feature,DSF)的提取[2]，常用的提取方法[3]包括頻域方法和時(shí)域方法。頻域方法通過(guò)實(shí)驗(yàn)?zāi)B(tài)分析，以頻率變化率[4]、模態(tài)曲率[5]、模態(tài)應(yīng)變能[6]等作為DSF進(jìn)行損傷識(shí)別。時(shí)域方法在時(shí)間域上對(duì)測(cè)試數(shù)據(jù)進(jìn)行分析，結(jié)合時(shí)間序列分析[7]、經(jīng)驗(yàn)?zāi)B(tài)分解[8]等方法來(lái)提取DSF。與頻域方法相比，時(shí)域方法由于不需要進(jìn)行繁瑣的頻譜分析，越來(lái)越受到廣泛關(guān)注[9～10]。

自回歸(autoregressive,AR)模型作為時(shí)域方法中的重要組成部分，常被用于回歸擬合結(jié)構(gòu)加速度響應(yīng)數(shù)據(jù)，并從中提取DSF。Lynch[11]通過(guò)對(duì)AR模型方程進(jìn)行Z變換，證明了AR模型系數(shù)與結(jié)構(gòu)固有特性相關(guān)。Nair等[12]以AR模型前三階系數(shù)構(gòu)造DSF，使用t檢驗(yàn)進(jìn)行了損傷識(shí)別及定位。Datteo等[13]通過(guò)全局敏感性分析和不確定分析，得出了AR模型系數(shù)對(duì)結(jié)構(gòu)剛度與質(zhì)量變化較為敏感，對(duì)阻尼與荷載的變化不敏感。杜永峰等以結(jié)構(gòu)損傷前后AR模型殘差方差之比作為損傷指標(biāo)，通過(guò)數(shù)值算例定性的實(shí)現(xiàn)了結(jié)構(gòu)損傷識(shí)別及定位。朱軍華等[14]以AR模型殘差的偏度和峰度作為DSF，取得了一定的識(shí)別效果?？梢姡瑢R模型用于DSF的提取展現(xiàn)出了極大的應(yīng)用潛力與科研價(jià)值。

獲取DSF后，結(jié)構(gòu)損傷識(shí)別就歸結(jié)為模式識(shí)別的問(wèn)題。機(jī)器學(xué)習(xí)作為模式識(shí)別的主要方法，包括神經(jīng)網(wǎng)絡(luò)[15]、支持向量機(jī)[16]、高斯過(guò)程等。高斯過(guò)程(Gaussian process,GP)是一種常用的監(jiān)督學(xué)習(xí)方法，旨在解決回歸問(wèn)題和概率分類問(wèn)題，包括高斯過(guò)程分類(Gaussian process classification,GPC)和高斯過(guò)程回歸(Gaussian process regression,GPR)。Sawant等[17]基于GP提出了一種分層算法，實(shí)現(xiàn)了人臉年齡估計(jì)。舒堅(jiān)等[18]利用GPR進(jìn)行了鏈路質(zhì)量預(yù)測(cè)，取得了比動(dòng)態(tài)貝葉斯網(wǎng)絡(luò)更精確的結(jié)果。朱彬等[19]提出了一種基于GPR的響應(yīng)面擬合方法，并實(shí)現(xiàn)了邊坡可靠度分析。程虎等[20]利用GPR實(shí)現(xiàn)了成像系統(tǒng)參數(shù)的在線優(yōu)化。由此可以看出，GP在生物識(shí)別、狀態(tài)預(yù)測(cè)、參數(shù)優(yōu)化等方面有較多應(yīng)用，但在損傷識(shí)別領(lǐng)域卻很少見。損傷識(shí)別分三個(gè)階段：一是識(shí)別損傷；二是定位損傷；三是識(shí)別損傷程度。實(shí)際過(guò)程中由于建模誤差、測(cè)量噪聲的影響，觀測(cè)值近似隨機(jī)過(guò)程，結(jié)合GP具有概率意義的輸出特性，將GPC用于損傷定位、GPR用于識(shí)別損傷程度，能使輸出結(jié)果更具信服力，因此將GP運(yùn)用于損傷識(shí)別凸顯了巨大的工程實(shí)用價(jià)值。

考慮到損傷識(shí)別的復(fù)雜性與既有方法的不足，本文提出了一種基于AR-GP模型的損傷識(shí)別方法。首先，利用AR模型回歸擬合結(jié)構(gòu)加速度響應(yīng)數(shù)據(jù)；然后，根據(jù)AR模型的殘差、系數(shù)分別構(gòu)造DSF；最后，借助GPC與GPR分別識(shí)別結(jié)構(gòu)損傷位置與損傷程度。通過(guò)某鋼筋混凝土模型拱的數(shù)值算例，驗(yàn)證了所提出的損傷識(shí)別方法。

1 自回歸模型

1.1 建模方法

如果一個(gè)隨機(jī)過(guò)程中的觀測(cè)量Yt可以表示為

Yt=c+φ1Yt-1+φ2Yt-2+…+φpYt-p+εt

(1)

AR模型建模步驟如下：(1)數(shù)據(jù)平穩(wěn)性檢驗(yàn)；(2)數(shù)據(jù)標(biāo)準(zhǔn)化處理；(3)模型識(shí)別；(4)參數(shù)估計(jì)。

數(shù)據(jù)平穩(wěn)性檢驗(yàn)即為單位根檢驗(yàn)，本文采用ADF檢驗(yàn)法[21]。ADF檢驗(yàn)包括三個(gè)檢驗(yàn)?zāi)Ｐ?/p>

(2)

三個(gè)檢驗(yàn)?zāi)Ｐ椭笑?、β、δ為最小二乘回歸系數(shù)。實(shí)際檢驗(yàn)過(guò)程中從模型3～模型1依次檢驗(yàn)，若有任何一個(gè)模型拒絕零假設(shè)，則認(rèn)為序列不存在單位根，即序列平穩(wěn)。各個(gè)模型原假設(shè)如下

(3)

數(shù)據(jù)標(biāo)準(zhǔn)化處理的目的是為了消除噪聲及荷載對(duì)測(cè)試數(shù)據(jù)的影響，計(jì)算公式如下

(4)

式中，u、σ分別為響應(yīng)數(shù)據(jù)的均值和標(biāo)準(zhǔn)差。

模型識(shí)別即確定AR模型滯后階數(shù)p，根據(jù)AR(p)模型偏自相關(guān)函數(shù)p階截尾的特點(diǎn)，可通過(guò)計(jì)算偏自相關(guān)函數(shù)來(lái)確定p。但實(shí)際過(guò)程中理論偏自相關(guān)函數(shù)難以獲取，故本文采用樣本偏自相關(guān)函數(shù)[22]與AIC準(zhǔn)則[23]聯(lián)合來(lái)確定p。

(5)

(6)

(7)

確定好AR模型的滯后階數(shù)p之后，選取最小二乘法進(jìn)行模型參數(shù)估計(jì)，基本思想為尋找使得模型殘差平方和最小的一組參數(shù)值，即：

(8)

1.2 基于殘差的損傷敏感性特征

建立好AR模型之后，即可提取DSF進(jìn)行損傷識(shí)別工作。綜合國(guó)內(nèi)外研究成果，常用的DSF大致分兩類：基于殘差的DSF和基于自回歸系數(shù)的DSF。本文分別建立了基于殘差和自回歸系數(shù)的DSF，并對(duì)其預(yù)測(cè)結(jié)果進(jìn)行了對(duì)比分析。

針對(duì)損傷識(shí)別過(guò)程中，定位損傷位置困難、多損傷狀態(tài)不易識(shí)別的問(wèn)題，引入了表征結(jié)構(gòu)損傷位置信息與損傷狀態(tài)信息的參數(shù)L1、L2

(9)

(10)

式中：i代表傳感器編號(hào)，i=1,2,…,d；DSFi表示第i個(gè)傳感器的DSF。

結(jié)構(gòu)發(fā)生損傷后，對(duì)應(yīng)AR模型殘差標(biāo)準(zhǔn)差(residual standard deviation,RSD)會(huì)發(fā)生變化，于是構(gòu)造基于殘差標(biāo)準(zhǔn)差的DSF

DSFi=σεi

式中,σεi表示第i個(gè)傳感器處殘差標(biāo)準(zhǔn)差，i=1,2,…,d。

在損傷定位時(shí)，引入?yún)?shù)L1、L2，定義基于殘差標(biāo)準(zhǔn)差的GPC輸入向量RSD1

RSD1=(σε1,σε2,…,σεd,L1,L2)

在識(shí)別損傷程度時(shí)，定義基于殘差標(biāo)準(zhǔn)差的GPR輸入向量RSD2

RSD2=(σε1,σε2,…,σεd)

1.3 基于自回歸系數(shù)的損傷敏感性特征

根據(jù)文獻(xiàn)[11]，AR模型系數(shù)與結(jié)構(gòu)固有特性存在一定的聯(lián)系，因此基于自回歸系數(shù)(autoregressive coefficient,ARC)的DSF表述如下

其中，φj,j=1,2,…,p為AR(p)模型系數(shù)；φi為第i個(gè)傳感器處的DSF，i=1,2,…,d。

同樣在損傷定位時(shí)，引入?yún)?shù)L1、L2，定義基于自回歸系數(shù)的GPC輸入向量ARC1

ARC1=(φ1,φ2,…,φd,L1,L2)

在識(shí)別損傷程度時(shí)，定義基于自回歸系數(shù)的GPR輸入向量ARC2

ARC2=(φ1,φ2,…,φd)

2 高斯過(guò)程

2.1 高斯過(guò)程回歸

GP是符合多維聯(lián)合正態(tài)分布隨機(jī)過(guò)程的總稱。對(duì)于給定的訓(xùn)練數(shù)據(jù)集

D={(X,y)|X={x1,x2,…,xN}T∈RN×d,y∈RN}

其中X、y分別為輸入和輸出變量，且N、d分別為輸入變量個(gè)數(shù)及輸入變量維數(shù)，y用于描述輸入變量的隨機(jī)過(guò)程，且服從GP，其性質(zhì)由均值函數(shù)E(y)、協(xié)方差函數(shù)K決定：

y～N(E(y),K)

其中K也叫做格拉姆矩陣，其元素為：

Kij=K(xi,xj),(i,j=1,2,…N)

高斯過(guò)程回歸(GPR)是基于貝葉斯理論和統(tǒng)計(jì)學(xué)習(xí)理論發(fā)展起來(lái)的一種新的機(jī)器學(xué)習(xí)方法，適用于處理高維數(shù)、小樣本、非線性等復(fù)雜回歸問(wèn)題，近幾年來(lái)在國(guó)內(nèi)外發(fā)展迅速[25]，其本質(zhì)是基于貝葉斯推斷，在給定輸入的情況下，給出輸出的條件分布，其基本步驟包括訓(xùn)練與預(yù)測(cè)。

(11)

GPR的關(guān)鍵是協(xié)方差函數(shù)K，即核函數(shù)的選取，包括線性核、多項(xiàng)式核等，使用頻率最高的是平方指數(shù)核，即：

(12)

(13)

其中k是K(xi,x*)(i=1,2,…,N)的N×1維向量，r=K(x*,x*)。

根據(jù)貝葉斯推斷，可得預(yù)測(cè)值y*的后驗(yàn)分布為

2.2 高斯過(guò)程分類

當(dāng)GPR中的連續(xù)型輸出變量變?yōu)殡x散的類別向量時(shí)，GPR就演化為高斯過(guò)程分類(Gaussian process classification,GPC)。

GPC的基本思想是首先賦予隱函數(shù)f(x)的高斯過(guò)程先驗(yàn)分布；然后引入反應(yīng)函數(shù)獲取先驗(yàn)分布π(x)=p(y=1|x)=σ(f(x))；最后輸出類別概率最高的分類結(jié)果。

對(duì)于給定的輸入x*，f*=f(x*)的分布以及類別預(yù)測(cè)概率分別按式(14)、(15)計(jì)算：

(14)

(15)

由于輸入變量非連續(xù)，似然函數(shù)p(y|f)服從非高斯分布，式(14)、(15)無(wú)法獲取解析解，常采用拉普拉斯近似方法求解[26]。

2.3 損傷識(shí)別流程

基于AR-GP模型的損傷識(shí)別方法的大致流程為：

(1) 獲取結(jié)構(gòu)在各個(gè)測(cè)點(diǎn)的加速度響應(yīng)數(shù)據(jù)；

(2) 根據(jù)1.1節(jié)所提建模方法，建立AR模型；

(3) 構(gòu)造DSF：計(jì)算RSD1、RSD2，ARC1、ARC2；

(4) 損傷識(shí)別：首先，以RSD1、ARC1為輸入向量，損傷位置為分類標(biāo)簽，利用GPC實(shí)現(xiàn)損傷定位。然后，以RSD2、ARC2為輸入向量，損傷程度為輸出，利用GPR識(shí)別損傷程度。算法流程如圖1所示。

圖1 AR-GP算法流程圖Fig.1 Flow chart of AR-GP algorithm

3 實(shí)例驗(yàn)證

3.1 鋼筋混凝土拱模型

本文所采用的數(shù)值模擬對(duì)象為某鋼筋混凝土模型試驗(yàn)拱，拱肋沿跨徑方向劃分成9個(gè)節(jié)段，共布置了7個(gè)加速度測(cè)點(diǎn)，其具體尺寸參數(shù)如圖2所示。

圖2 模型拱基本構(gòu)造Fig.2 Basic structure of model arch

利用ANSYS通用有限元軟件建立結(jié)構(gòu)有限元模型，采用分離式建模方法?；炷敛捎肧OLID65單元，加載塊采用SOLID45單元，鋼筋采用LINK8單元。有限元模型共劃分為2 888個(gè)單元，如圖3所示。

圖3 ANSYS有限元模型Fig.3 Finite element model of ANSYS

邊界條件以面約束的形式施加，約束兩拱腳所有平動(dòng)及轉(zhuǎn)動(dòng)位移；配重通過(guò)加載塊的自重來(lái)施加，并根據(jù)其實(shí)際重量及加載塊體積換算等效密度。隨后，進(jìn)行了瞬態(tài)分析。

瞬態(tài)分析采用完全法，結(jié)構(gòu)阻尼采用瑞利阻尼，質(zhì)量阻尼系數(shù)α與剛度阻尼系數(shù)β按下式計(jì)算

(16)

式中:ωi和ωj分別為結(jié)構(gòu)第i和第j階固有頻率，可通過(guò)模態(tài)分析獲??；ξi和ξj為相對(duì)應(yīng)于第i和第j階振型的阻尼比；一般取i=1，j=2，ξi=ξj=0.05。

在拱肋7L/16處施加白噪聲激勵(lì)，采樣頻率為1 000 Hz，采樣時(shí)間為0.5 s。

3.2 損傷工況

結(jié)構(gòu)損傷顯著影響結(jié)構(gòu)的剛度，剛度退化通過(guò)彈性模量的折減來(lái)實(shí)現(xiàn)。本文考慮了36種損傷工況，包括單損傷、多損傷，以及不同損傷程度的情況，其中有24種(A1～A24)訓(xùn)練工況以及12種(B1～B12)測(cè)試工況，損傷工況見表1。

表1 損傷工況Tab.1 Damage scenarios

為模擬實(shí)際的測(cè)試環(huán)境，對(duì)有限元提取的加速度響應(yīng)數(shù)據(jù)施加了一定程度的白噪聲干擾，按下式計(jì)算[27]

anoise=a+RMS(a)×Nlevel×Nunit

(17)

式中：anoise、a分別為含有噪聲與未含噪聲的加速度響應(yīng)數(shù)據(jù)；RMS(a)為a的均方根；Nunit為高斯白噪聲；Nlevel為所施加的噪聲水平，本文選為5%和10%。

3.3 結(jié)果分析

根據(jù)加速度響應(yīng)建立AR模型，然后計(jì)算各個(gè)工況、各個(gè)測(cè)點(diǎn)所對(duì)應(yīng)的RSD1、RSD2、ARC1、ARC2，核函數(shù)采用平方指數(shù)函數(shù)與白噪聲函數(shù)相加的形式，基于訓(xùn)練樣本，利用sklearn機(jī)器學(xué)習(xí)庫(kù)對(duì)GPC、GPR進(jìn)行訓(xùn)練，然后將測(cè)試樣本送入已經(jīng)訓(xùn)練好的GPC、GPR進(jìn)行損傷定位及識(shí)別損傷程度。

3.3.1 損傷定位結(jié)果

基于殘差標(biāo)準(zhǔn)差和自回歸系數(shù)的損傷定位結(jié)果見表2。

從表2的損傷定位結(jié)果可知，隨著噪聲水平的增加，識(shí)別精度逐漸降低，但總體而言，基于殘差標(biāo)準(zhǔn)差的識(shí)別方法具有較好的識(shí)別精度以及抗噪性，且對(duì)于單損傷和多損傷狀態(tài)有較好的區(qū)分能力；而基于自回歸系數(shù)的識(shí)別方法則表現(xiàn)出較差的識(shí)別精度與抗噪性。

表2 損傷定位結(jié)果Tab.2 Results of damage localization

3.3.2 損傷程度識(shí)別結(jié)果

基于殘差標(biāo)準(zhǔn)差和自回歸系數(shù)的損傷程度識(shí)別結(jié)果分別見表3、表4。表中:u表示識(shí)別結(jié)果的均值；σ表示識(shí)別結(jié)果的標(biāo)準(zhǔn)差；RE表示識(shí)別結(jié)果的相對(duì)誤差；DC表示識(shí)別結(jié)果的離散系數(shù)。

表4 基于自回歸系數(shù)的損傷程度識(shí)別結(jié)果Tab.4 Results of damage severity identification based on autoregressive coefficient

從表3的損傷程度識(shí)別結(jié)果可以看出，隨著噪聲水平的增加，RE與DC的均值變化不大，且處于較低的水準(zhǔn)，在10%噪聲污染的情況下，識(shí)別結(jié)果的RE與DC值也僅僅為6.52%和0.19。DC越接近于0，表明識(shí)別結(jié)果波動(dòng)性越小、越可靠，故可以認(rèn)為基于殘差標(biāo)準(zhǔn)差的方法有很高的識(shí)別精度及可靠度，且抗噪性能好。

表3 基于殘差標(biāo)準(zhǔn)差的損傷程度識(shí)別結(jié)果Tab.3 Results of damage severity identification based on residual standard deviation

從表4可知，與基于殘差標(biāo)準(zhǔn)差的識(shí)別結(jié)果相比，基于自回歸系數(shù)的識(shí)別精度及可靠度則相對(duì)較低，抗噪性能差。

圖4是所有識(shí)別結(jié)果的RE與DC散點(diǎn)圖，從圖中可以看出RE與DC呈現(xiàn)出較強(qiáng)的正相關(guān)性。在實(shí)際的監(jiān)測(cè)過(guò)程中，由于無(wú)法知曉預(yù)測(cè)結(jié)果的真值，也就不能獲取預(yù)測(cè)結(jié)果的RE。此時(shí)，可以根據(jù)DC判斷預(yù)測(cè)結(jié)果的精確度以及波動(dòng)程度，由此做出預(yù)測(cè)結(jié)果是否可接受的判斷，從而提高了識(shí)別可信度。

圖4 RE與DC關(guān)系Fig.4 The relationship between RE and DC

另外，利用GPR輸出結(jié)果概率分布的特點(diǎn)，在實(shí)際結(jié)構(gòu)的監(jiān)測(cè)過(guò)程中可以實(shí)現(xiàn)損傷預(yù)警，如圖5所示，其基本思想為：如果預(yù)測(cè)結(jié)果有1-α的概率超過(guò)給定的損傷程度，則發(fā)出預(yù)警。圖中:η為損傷程度;ηt為損傷閾值;ηα為根據(jù)特定問(wèn)題所規(guī)定的下α分位點(diǎn)。

圖5 損傷預(yù)警Fig.5 Damage alarming

4 結(jié) 論

本文提出了一種基于AR-GP模型的全新的結(jié)構(gòu)損傷識(shí)別方法。該方法利用AR模型回歸擬合結(jié)構(gòu)加速度響應(yīng)數(shù)據(jù)，并從中提取DSF，首次引入損傷位置信息與損傷狀態(tài)信息參數(shù)L1、L2，基于GPC實(shí)現(xiàn)了結(jié)構(gòu)多損傷定位，借助于GPR實(shí)現(xiàn)了損傷程度的概率結(jié)果輸出。通過(guò)36種不同損傷工況的數(shù)值算例，驗(yàn)證了所提方法，并對(duì)基于殘差標(biāo)準(zhǔn)差和自回歸系數(shù)的識(shí)別結(jié)果進(jìn)行了對(duì)比分析，得出了以下結(jié)論：

(1) 在損傷定位方面，所提識(shí)別方法能夠較好的區(qū)分單損傷和多損傷狀態(tài)，同時(shí)基于殘差標(biāo)準(zhǔn)差的定位結(jié)果較基于自回歸系數(shù)的結(jié)果有更好的識(shí)別精度與抗噪性能。在識(shí)別損傷程度方面，基于殘差標(biāo)準(zhǔn)差的識(shí)別結(jié)果表現(xiàn)出了良好的識(shí)別精度、可靠度與抗噪能力，在10%噪聲污染的情況下，識(shí)別結(jié)果的RE與DC值僅為6.52%和0.19，而基于自回歸系數(shù)的識(shí)別結(jié)果的精度與抗噪能力則較差。

(2) 損傷程度識(shí)別結(jié)果的RE與DC展示出了較強(qiáng)的正相關(guān)性，在實(shí)際的監(jiān)測(cè)過(guò)程中，可以由DC判斷預(yù)測(cè)結(jié)果的精確度以及波動(dòng)程度，能夠提高識(shí)別可靠性。

(3) 利用GPR輸出預(yù)測(cè)結(jié)果均值與方差的特性，建立相應(yīng)的損傷閾值，從而可以實(shí)現(xiàn)損傷預(yù)警。值得一提的是，可通過(guò)可靠度理論，按一定的失效概率構(gòu)造損傷閾值，這也是未來(lái)的研究方向。