基于變量選擇的尖點(diǎn)突變模型的兩步構(gòu)建方法

張 明，付冬梅，程學(xué)群，楊丙坤，郝文魁，陳 云，邵立珍

1) 北京科技大學(xué)順德研究生院，佛山 528300 2) 北京科技大學(xué)自動(dòng)化學(xué)院，北京 100083 3) 北京科技大學(xué)北京市工業(yè)波譜成像工程技術(shù)研究中心，北京 100083 4) 北京科技大學(xué)新材料技術(shù)研究院，北京 100083 5) 北京科技大學(xué)國(guó)家材料腐蝕與防護(hù)科學(xué)數(shù)據(jù)中心，北京100083 6) 全球能源互聯(lián)網(wǎng)研究院有限公司先進(jìn)輸電技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 102209

在復(fù)雜的系統(tǒng)中，外界因素的變化可能會(huì)導(dǎo)致系統(tǒng)狀態(tài)的跳躍式變化，稱(chēng)為突變.Qiao 等[1]認(rèn)為深埋隧道圍巖的失穩(wěn)是一種突變現(xiàn)象，給礦井的安全生產(chǎn)帶來(lái)極大威脅，并通過(guò)分析得出圍巖失穩(wěn)發(fā)生在第3 至4 步開(kāi)挖過(guò)程中.Zhi 等[2]通過(guò)分析得到使腐蝕速率急劇變化的環(huán)境變量的閾值，當(dāng)環(huán)境變量超過(guò)該閾值時(shí)，腐蝕速率會(huì)發(fā)生突變.裴甲坤等[3]認(rèn)為化工事故是由危險(xiǎn)源和不安全因素引起的突發(fā)事件，二者的綜合影響導(dǎo)致系統(tǒng)的安全狀態(tài)發(fā)生突變.其他諸如股市崩盤(pán)[4-5]、人的心理狀態(tài)變化[6]、電力系統(tǒng)故障[7]等也屬于突變現(xiàn)象.此類(lèi)現(xiàn)象包含復(fù)雜的系統(tǒng)行為，既有連續(xù)性變化又有突發(fā)的不連續(xù)性變化，且影響因素往往眾多，給實(shí)際工程問(wèn)題的建模和解釋帶來(lái)困難.

對(duì)于小樣本工程數(shù)據(jù)，線(xiàn)性模型、灰色模型[8]是常用的方法；對(duì)于大樣本工程數(shù)據(jù)，人工神經(jīng)網(wǎng)絡(luò)[9]、隨機(jī)森林(Random forest,RF)[10]等機(jī)器學(xué)習(xí)模型通常可以獲得較好的建模效果.雖然機(jī)器學(xué)習(xí)模型具有強(qiáng)大的非線(xiàn)性映射能力，但無(wú)法解釋研究對(duì)象的突變現(xiàn)象.突變理論是用以解釋復(fù)雜系統(tǒng)中不連續(xù)性和質(zhì)變現(xiàn)象的數(shù)學(xué)理論，由法國(guó)數(shù)學(xué)家Thom[11]提出.假定一個(gè)系統(tǒng)的動(dòng)力學(xué)方程可以由一個(gè)光滑的勢(shì)函數(shù)導(dǎo)出，根據(jù)控制因子和狀態(tài)因子個(gè)數(shù)的不同，Thom 定義了7 種基本的突變模型，并推導(dǎo)出每一種模型勢(shì)函數(shù)的解析形式.由于形式簡(jiǎn)單、直觀(guān)，具有兩個(gè)控制因子和一個(gè)狀態(tài)因子的尖點(diǎn)突變模型應(yīng)用最為廣泛，模型參數(shù)的估計(jì)可由Cobb 提出的極大似然估計(jì)法(Maximum likelihood estimation,MLE)實(shí)現(xiàn)[12-13].作為解決工程領(lǐng)域中不連續(xù)性復(fù)雜問(wèn)題的一種數(shù)學(xué)工具，突變理論在經(jīng)濟(jì)學(xué)、生物學(xué)、物理學(xué)、心理學(xué)等領(lǐng)域應(yīng)用廣泛.

在以往的尖點(diǎn)突變模型中，組成控制因子的輸入變量往往依據(jù)經(jīng)驗(yàn)或已有的結(jié)論來(lái)確定.如，文獻(xiàn)[4-5]基于Zeeman[14]的理論基礎(chǔ)，將股票市場(chǎng)中基本面交易者和技術(shù)分析交易者的多維數(shù)據(jù)作為輸入變量構(gòu)建股票市場(chǎng)的尖點(diǎn)突變模型.這種建模方式受限于特定學(xué)科，不利于推廣，且輸入變量的實(shí)際價(jià)值難以判斷.在待選的輸入變量較多且突變機(jī)理不明確的情況下，如何利用少量的重要變量構(gòu)建尖點(diǎn)突變模型依然是一個(gè)難點(diǎn).常見(jiàn)的變量選擇方法分為過(guò)濾法、嵌入法、封裝法.其中，過(guò)濾法依據(jù)待選變量統(tǒng)計(jì)特性的各項(xiàng)指標(biāo)來(lái)選擇重要變量，如皮爾遜相關(guān)系數(shù)法、方差過(guò)濾法；嵌入法依據(jù)機(jī)器學(xué)習(xí)算法本身來(lái)分析待選變量的重要性，如RF 的排列變量重要性算法[10,15]；封裝法基于構(gòu)造的最終模型來(lái)選擇使模型性能達(dá)到最優(yōu)的變量子集，最終模型可以是支持向量機(jī)(Support vector regression,SVR)[16]、梯度提升回歸樹(shù)(Gradient boosted regression trees,GBRT)[17]等機(jī)器學(xué)習(xí)算法.過(guò)濾法的評(píng)價(jià)標(biāo)準(zhǔn)獨(dú)立于特定的學(xué)習(xí)算法，具有較好的通用性，但難以取得很好的建模效果；嵌入法、封裝法雖然可以取得較好的建模效果，但是這種基于單一模型的變量選擇算法存在特定的偏差，變量子集的選取依賴(lài)于特定模型，容易產(chǎn)生過(guò)擬合現(xiàn)象.采用集成方法，即通過(guò)組合不同方法的變量選擇結(jié)果來(lái)產(chǎn)生變量子集，既減輕了對(duì)特定模型的依賴(lài)性，又可以很好地提高結(jié)果的準(zhǔn)確性和穩(wěn)定性[18-20].

針對(duì)傳統(tǒng)尖點(diǎn)突變模型依據(jù)經(jīng)驗(yàn)建模的問(wèn)題，提出基于變量選擇的尖點(diǎn)突變模型的兩步構(gòu)建方法.該方法的通用性較強(qiáng)，可廣泛應(yīng)用于具有突變特征的系統(tǒng)的建模并能得到模型的數(shù)學(xué)解析式.建模過(guò)程分為兩步.第一步，以RF、GBRT、SVR 作為基學(xué)習(xí)器，利用多模型集成重要變量選擇算法(Multi-model ensemble important variable selection,MEIVS)來(lái)量化待選變量的重要性，提取得分之和超過(guò)總分90%的前n個(gè)待選變量作為后續(xù)建模的輸入變量；第二步，基于MLE 算法構(gòu)建尖點(diǎn)突變模型.本文首先介紹了尖點(diǎn)突變模型的原理、數(shù)據(jù)擬合方法以及突變特征，其次介紹了MEIVS 算法的實(shí)現(xiàn)流程，最后結(jié)合工程實(shí)例，驗(yàn)證了該方法的有效性.

1 基本原理

1.1 尖點(diǎn)突變模型與突變特征

1.1.1 尖點(diǎn)突變模型

突變理論描述了動(dòng)力學(xué)系統(tǒng)中控制因子和狀態(tài)因子之間的關(guān)系，在控制因子固定的情況下，系統(tǒng)始終尋求平衡狀態(tài)，直到達(dá)到勢(shì)函數(shù)的極小值或極大值為止.以動(dòng)力學(xué)系統(tǒng)表達(dá)式來(lái)描述系統(tǒng)的狀態(tài)因子z在控制因子a的影響下隨時(shí)間t的變化：

V(z;a)是系統(tǒng)的勢(shì)函數(shù).應(yīng)用最廣泛的尖點(diǎn)突變模型由2 個(gè)控制因子α、β和一個(gè)狀態(tài)因子z組成，其勢(shì)函數(shù)的規(guī)范形式是：

系統(tǒng)的平衡方程由(3)式確定，在無(wú)擾動(dòng)的情況下，系統(tǒng)的狀態(tài)不隨時(shí)間變化:

當(dāng)平衡點(diǎn)的勢(shì)函數(shù)V(z;α,β)是關(guān)于z的極小值時(shí)，平衡點(diǎn)是穩(wěn)定的，系統(tǒng)即使受到擾動(dòng)的影響，也會(huì)隨著時(shí)間t回到穩(wěn)定狀態(tài)；當(dāng)平衡點(diǎn)的勢(shì)函數(shù)V(z;α,β)是關(guān)于z的極大值時(shí)，平衡點(diǎn)是不穩(wěn)定的，系統(tǒng)在擾動(dòng)的影響下會(huì)偏離此平衡點(diǎn)，從而被穩(wěn)定的平衡點(diǎn)吸引.在不同的α和β值下系統(tǒng)平衡點(diǎn)的數(shù)目和性質(zhì)可以由Cardan 判別式δ判斷，表示為：

當(dāng)δ＞0 時(shí)，存在一個(gè)穩(wěn)定的平衡點(diǎn)；當(dāng)δ＜0 時(shí)，存在兩個(gè)穩(wěn)定的平衡點(diǎn)和一個(gè)不穩(wěn)定的平衡點(diǎn)；當(dāng)δ=0 時(shí)，存在一個(gè)穩(wěn)定的平衡點(diǎn)和一個(gè)不穩(wěn)定的平衡點(diǎn).圖1 給出了由平衡點(diǎn)的集合構(gòu)成的平衡曲面和由控制因子構(gòu)成的控制平面.平衡曲面的形狀像一個(gè)有“褶皺”的連續(xù)曲面，并且由上葉、中葉、下葉3 部分構(gòu)成，上葉和下葉部分對(duì)應(yīng)的平衡點(diǎn)是穩(wěn)定的，中葉部分對(duì)應(yīng)的平衡點(diǎn)是不穩(wěn)定的.控制平面是平衡曲面在z軸方向上的投影，中葉區(qū)域在控制平面上的投影稱(chēng)為尖點(diǎn)突變模型的分叉集.圖1 中，若控制因子α、β沿紅色軌跡A 變化，狀態(tài)因子z會(huì)在分叉集內(nèi)發(fā)生突變，從平衡曲面的下葉直接跳變到上葉而不經(jīng)過(guò)中葉；若控制因子α、β沿藍(lán)色軌跡B 變化，則狀態(tài)因子z不會(huì)發(fā)生突變.

圖1 尖點(diǎn)突變模型的平衡曲面和控制平面Fig.1 Equilibrium surface and control plane of the cusp catastrophe model

在實(shí)際應(yīng)用中，數(shù)據(jù)難免受到隨機(jī)噪聲影響，Cobb 和Zacks[12-13]通過(guò)引入隨機(jī)微分方程，以概率密度函數(shù)的形式描述了系統(tǒng)在α、β固定的條件下z的分布，表示為：

式中，ψ是歸一化常數(shù)，α、β分別為n維輸入變量{X1,···,Xn}的線(xiàn)性組合，z為輸出變量Y的線(xiàn)性變換，表示為：

參數(shù)θ={w0,w1,a0,···,an,b0,···,bn}由MLE 方法估計(jì).在給定N個(gè)觀(guān)測(cè)樣本的情況下，參數(shù)θ的對(duì)數(shù)似然函數(shù)如下：

最大化觀(guān)測(cè)樣本的對(duì)數(shù)似然函數(shù)得到參數(shù)θ的估計(jì)值，優(yōu)化搜索算法為帶邊界約束的Broyden-Fletcher-Goldfarb-Shanno 算法：

1.1.2 評(píng)價(jià)指標(biāo)

為了驗(yàn)證尖點(diǎn)突變模型的性能，將線(xiàn)性模型和非線(xiàn)性的Logistic 模型與尖點(diǎn)突變模型作對(duì)比[21].其中，Logistic 模型可以模擬研究對(duì)象的急劇變化，但沒(méi)有考慮不連續(xù)性變化.評(píng)價(jià)指標(biāo)采用可決系數(shù),R2、赤池信息準(zhǔn)則(Akaike information criterion,AIC)和貝葉斯信息準(zhǔn)則(Bayesian information criterion,BIC).當(dāng)R2越大時(shí)，模型的精度越高；AIC 和BIC 考慮了模型復(fù)雜度，當(dāng)AIC 和BIC 越小時(shí)，模型越好.

1.1.3 突變特征

在系統(tǒng)的勢(shì)函數(shù)未知的情況下，常常根據(jù)系統(tǒng)表現(xiàn)的外部性態(tài)來(lái)判斷系統(tǒng)是否存在突變，這些性態(tài)被稱(chēng)為突變特征[14,21].尖點(diǎn)突變有5 個(gè)特征：(1)多模態(tài)：系統(tǒng)中可能出現(xiàn)兩個(gè)不同的狀態(tài)；(2)不可達(dá)性：系統(tǒng)存在不穩(wěn)定的平衡態(tài)；(3)突跳：系統(tǒng)從一個(gè)勢(shì)函數(shù)極小值跳到另一個(gè)極小值；(4)發(fā)散：控制因子的微小變化可以導(dǎo)致?tīng)顟B(tài)因子的質(zhì)變；(5)滯后：當(dāng)物理過(guò)程可逆時(shí)，發(fā)生突變時(shí)對(duì)應(yīng)的控制參數(shù)位置可能不同.當(dāng)系統(tǒng)存在突變現(xiàn)象時(shí)，對(duì)外往往表現(xiàn)為其中的一個(gè)或幾個(gè)的組合.在實(shí)際應(yīng)用中，針對(duì)截面數(shù)據(jù)，應(yīng)首先檢查研究對(duì)象概率密度的雙峰性，雙峰性意味著系統(tǒng)可能存在多個(gè)狀態(tài)；針對(duì)時(shí)序數(shù)據(jù)，則應(yīng)首先檢查時(shí)間序列中的跳變現(xiàn)象[21].

1.2 多模型集成重要變量選擇算法

而在傳統(tǒng)的尖點(diǎn)突變模型的建模過(guò)程中，輸入變量的選取往往依賴(lài)于已有的實(shí)踐或經(jīng)驗(yàn)，這與目前數(shù)據(jù)規(guī)模的爆發(fā)式增長(zhǎng)相矛盾，不利于尖點(diǎn)突變模型的普及應(yīng)用.為了解決上述問(wèn)題，同時(shí)提高模型的精度、降低模型的復(fù)雜度，本文基于排列[23]的思想提出MEIVS 算法.

排列的思想借鑒于隨機(jī)森林的變量重要性度量方法，認(rèn)為模型會(huì)更依賴(lài)于重要的輸入變量做預(yù)測(cè).當(dāng)打亂某一變量在測(cè)試集上的觀(guān)測(cè)序列后，用新生成的數(shù)據(jù)做預(yù)測(cè)，更重要的輸入變量會(huì)使模型的精度損失更大.MEIVS 算法組合了RF、GBRT、SVR 3 種常用的機(jī)器學(xué)習(xí)算法，其中RF 和GBRT 都屬于決策樹(shù)的集成學(xué)習(xí)算法，但它們采用的計(jì)算策略不同；SVR 采用高斯核函數(shù).文獻(xiàn)[24-26]中對(duì)每種方法的機(jī)理都作了解釋.本文的損失函數(shù)采用的是均方根誤差(Root mean squared error,RMSE):

以樣本的80%作為訓(xùn)練集，20%作為測(cè)試集，使用Z-Score 標(biāo)準(zhǔn)化方法對(duì)輸入變量進(jìn)行處理，經(jīng)過(guò)處理的數(shù)據(jù)的均值為0，標(biāo)準(zhǔn)差為1.記m個(gè)待選變量的集合為{S1,···,Sm}，目標(biāo)是得到n個(gè)重要變量的集合{X1,···,Xn}作為尖點(diǎn)突變模型的輸入變量.算法步驟如下，流程圖如圖2 所示.

圖2 MEIVS 算法流程圖.(a) MEIVS 算法主流程;(b) 排列算法流程Fig.2 MEIVS algorithm flowchart: (a) main process steps of the MEIVS algorithm;(b) process steps of the permutation algorithm

步驟1 利用訓(xùn)練集訓(xùn)練RF、GBRT、SVR 模型，記為M1、M2、M3，對(duì)于所建立的每個(gè)模型Mi，分別基于置換算法計(jì)算變量重要性，即執(zhí)行步驟2、步驟3；

步驟2 計(jì)算模型Mi在測(cè)試集上的均方根誤差并記為，對(duì){S1,···,Sm}，依次執(zhí)行(1)～(3)：

(1) 打亂Sj在測(cè)試集上的觀(guān)測(cè)序列并重新計(jì)算模型的均方根誤差，由于涉及隨機(jī)性，此過(guò)程重復(fù)10 次，分別記為

(2) 計(jì)算Sj在測(cè)試集上的平均預(yù)測(cè)精度損失：

(3) 計(jì)算Sj在模型Mi上的排列重要性得分，當(dāng)時(shí)，將重要性得分記為0，該變量無(wú)用：

步驟3 計(jì)算Sj的標(biāo)準(zhǔn)化排列重要性得分:

步驟4 計(jì)算Sj在M1、M2、M3上的重要性總得分:

步驟5 按變量重要性得分{V1,···,Vn}降序排列待選變量{S1,···,Sm}，提取得分之和超過(guò)總分90%的前n個(gè)待選變量作為重要變量，記為{X1,···,Xn}.

1.3 基于變量選擇的尖點(diǎn)突變模型的兩步構(gòu)建方法

將MEIVS 方法與基于MLE 的尖點(diǎn)突變模型參數(shù)估計(jì)方法相結(jié)合，分兩步構(gòu)建尖點(diǎn)突變模型.第一步，利用MEIVS 來(lái)量化待選變量{S1,···,Sm}的重要性，提取重要變量{X1,···,Xn}；第二步，利用提取的n個(gè)重要變量，基于MLE 算法構(gòu)建尖點(diǎn)突變模型.

2 仿真結(jié)果和分析

以?xún)蓚€(gè)不同領(lǐng)域的、具有突變特征的數(shù)據(jù)集為例，驗(yàn)證了所提方法的有效性.其中，歐洲旅館住宿價(jià)格數(shù)據(jù)集[27]為截面數(shù)據(jù)集，來(lái)源于Kaggle平臺(tái)；北京大氣腐蝕數(shù)據(jù)集為時(shí)序數(shù)據(jù)集，來(lái)源于北京地區(qū)的大氣暴露實(shí)驗(yàn).

2.1 在截面數(shù)據(jù)集上的應(yīng)用—以歐洲旅館住宿價(jià)格數(shù)據(jù)集為例

Kaggle 平臺(tái)的歐洲旅館住宿價(jià)格數(shù)據(jù)集一共包含120 個(gè)樣本，每個(gè)樣本包括每日住宿價(jià)格(Price)、星級(jí)(Star)、離市中心距離(Distance)、評(píng)分(Rating)、房間數(shù)目(Room)、房間面積(Square)和所在城市(City).Price 為輸出變量，單位為每日花費(fèi)的可兌換馬克(KM·d-1)，其余為輸入變量，其中類(lèi)別變量City 以Price 的類(lèi)別均值來(lái)編碼.Price 概率密度的非參數(shù)估計(jì)如圖3，非參估計(jì)的核函數(shù)選用高斯核，帶寬設(shè)置為25，概率密度的雙峰性暗示了Price 可能會(huì)發(fā)生突變，因此適用于建立尖點(diǎn)突變模型.兩步構(gòu)建方法中第一步為提取重要變量.利用MEIVS 得到各個(gè)待選變量的重要性得分，如圖4.條形圖中橫軸表示影響每日住宿價(jià)格的待選變量，縱軸表示每個(gè)待選變量的重要性總得分，每個(gè)待選變量在各模型上的得分以不同的顏色區(qū)分，并且根據(jù)得分降序排列.依據(jù)MEIVS 算法中步驟(5)，Square、Rating、Star 和Room 為重要變量,設(shè)為X1、X2、X3、X4，每日住宿價(jià)格Price 設(shè)為Y.算法基于R 語(yǔ)言中的DALEX 程序包[28]實(shí)現(xiàn).

圖3 每日住宿價(jià)格的概率密度非參數(shù)估計(jì)Fig.3 Nonparametric estimation of the probability density of the daily accommodation price

圖4 歐洲旅館住宿價(jià)格數(shù)據(jù)集待選變量重要性得分Fig.4 Importance score of the variables to be selected in the European hotel accommodation price dataset

將MEIVS 提取的重要變量X1、X2、X3、X4作為輸入變量、每日住宿價(jià)格Y作為輸出變量建立尖點(diǎn)突變模型，為了消除變量間量綱的影響，用ZScore 標(biāo)準(zhǔn)化方法對(duì)原始輸入變量進(jìn)行處理.算法基于R 語(yǔ)言Cusp 程序包[22]實(shí)現(xiàn).利用MLE 算法和120 條樣本對(duì)參數(shù)θ={w0,w1,a0,a1,a2,a3,a4,b0,b1,b2,b3,b4}進(jìn)行估計(jì)，代入式(6)中，得到如下形式的尖點(diǎn)突變模型的平衡方程：

表1 展示了采用兩步構(gòu)建法建立的尖點(diǎn)突變模型與經(jīng)MEIVS 降維后構(gòu)建的線(xiàn)性模型、Logistic 模型的評(píng)價(jià)指標(biāo)，同時(shí)與傳統(tǒng)的直接建模方法、經(jīng)斯皮爾曼相關(guān)系數(shù)(Spearman's correlation coefficient,SCC)、最大互信息系數(shù)(Maximal information coefficient,MIC)、隨機(jī)森林變量重要性算法(Random forest variable importance measure,RFVIM)降維的建模方法作比較.其中，SCC 和MIC 剔除系數(shù)小于0.3 的弱相關(guān)變量，RFVIM 提取累計(jì)變量重要性達(dá)到90%的前n個(gè)變量.結(jié)果顯示，在考慮樣本量的情況下，更高的R2和更低的BIC 說(shuō)明基于兩步構(gòu)建法所構(gòu)建的尖點(diǎn)突變模型優(yōu)于未降維的傳統(tǒng)尖點(diǎn)突變模型以及經(jīng)SCC、MIC、RFVIM降維后所構(gòu)建的尖點(diǎn)突變模型.

表1歐洲旅館住宿價(jià)格數(shù)據(jù)集建模結(jié)果評(píng)價(jià)Table 1 Evaluation of the modeling results of the European hotel accommodation price dataset

圖5(a)展示了樣本在控制平面上的分布，其中散點(diǎn)的顏色代表經(jīng)過(guò)式(6)線(xiàn)性變換后旅館價(jià)格的數(shù)值大小.影響旅館價(jià)格的控制因子的變化軌跡從左到右穿過(guò)了分叉集，表明旅館的價(jià)格發(fā)生了突變.圖5(b)展示了樣本在平衡曲面上的分布，平衡曲面設(shè)置為半透明狀態(tài)，顏色較暗的散點(diǎn)位于平衡曲面下方.易觀(guān)察到在較低的價(jià)格范圍內(nèi)旅館價(jià)格的變化具有連續(xù)性，而從低價(jià)到高價(jià)的變化并不連續(xù).

圖5 歐洲旅館住宿價(jià)格數(shù)據(jù)在控制平面 (a) 和平衡曲面(b) 上的分布Fig.5 Distribution of the European hotel accommodation price dataset on the control plane (a) and equilibrium surface (b)

2.2 在時(shí)序數(shù)據(jù)集上的應(yīng)用—以北京大氣腐蝕數(shù)據(jù)為例

北京大氣腐蝕數(shù)據(jù)集一共包含719 個(gè)樣本，采集時(shí) 間為2018 年8 月5 日16 時(shí)至9 月6 日14 時(shí)，采集地點(diǎn)為北京，每個(gè)樣本包括大氣腐蝕監(jiān)測(cè)儀(Atmospheric corrosion monitor,ACM)采集得到的早期大氣腐蝕電偶電流(Galvanic current)、溫度(T)、相對(duì)濕度(RH)、降雨?duì)顟B(tài)(Rainfall)以及大氣環(huán)境中PM2.5、PM10、SO2、NO2、O3的濃度.電偶電流與腐蝕速率成正相關(guān)關(guān)系[29]，為了便于分析，取電偶電流的自然對(duì)數(shù)作為輸出變量.圖6 展示了電偶電流的時(shí)間序列表明腐蝕電偶電流波動(dòng)較大，表明時(shí)間序列中具有突變的特性，因此適用于建立尖點(diǎn)突變模型.

圖6 ACM 采集到的電偶電流時(shí)間序列Fig.6 Time series of the galvanic current collected by ACM

通過(guò)MEIVS 算法得到待選變量重要性得分如圖7,可見(jiàn)T、RH 和Rainfall 為影響早期大氣腐蝕的重要變量，設(shè)為X1、X2、X3，對(duì)數(shù)化腐蝕電偶電流設(shè)為Y.其他污染物濃度的影響是微弱的.

圖7 北京大氣腐蝕數(shù)據(jù)集待選變量重要性得分Fig.7 Importance score of the variables to be selected in the Beijing atmospheric corrosion dataset

以Z-Score 標(biāo)準(zhǔn)化后的X1、X2、X3作為輸入變量、Y作為輸出變量構(gòu)建尖點(diǎn)突變模型，利用MLE算法和719 條樣本對(duì)參數(shù)θ={w0,w1,a0,a1,a2,a3,b0,b1,b2,b3}進(jìn)行估計(jì)，得到的平衡方程如下：

此外，采用上文所述方法，表2 的模型評(píng)估結(jié)果顯示了兩步構(gòu)建法的優(yōu)越性.

表2 北京大氣腐蝕數(shù)據(jù)集建模結(jié)果評(píng)價(jià)Table 2 Evaluation of the modeling results of the Beijing atmospheric corrosion dataset

樣本在控制平面和平衡曲面的分布情況如圖8(a)、8(b)，圓點(diǎn)代表未降雨時(shí)的樣本，三角形代表降雨時(shí)的樣本.當(dāng)由溫度、相對(duì)濕度、降雨組成的控制因子進(jìn)入分叉集時(shí)，腐蝕電偶電流在平衡曲面的下葉和上葉之間跳躍.從圖8(a)觀(guān)測(cè)到，降雨會(huì)促使腐蝕系統(tǒng)中的電偶電流不能沿著原有的軌跡運(yùn)動(dòng)，而是突變到新的演變軌跡上.

圖8 北京大氣腐蝕數(shù)據(jù)在控制平面(a)和平衡曲面(b)上的分布Fig.8 Distribution of the Beijing atmospheric corrosion dataset on the control plane (a) and equilibrium surface (b)

3 結(jié)論

(1) 對(duì)于存在突變現(xiàn)象的系統(tǒng)，通過(guò)理論和數(shù)據(jù)相結(jié)合的方式建立尖點(diǎn)突變模型是一種有效的建模手段.

(2) 提出了基于變量選擇的尖點(diǎn)突變模型的兩步構(gòu)建方法.在具有突變特征的數(shù)據(jù)集上，相比于其他模型，利用本文所提方法構(gòu)建的尖點(diǎn)突變模型擬合效果更優(yōu).

(3) 結(jié)合樣本在控制平面和平衡曲面的分布圖，尖點(diǎn)突變模型可以解釋系統(tǒng)的突變行為.