基于監(jiān)督降維和自適應(yīng)Kriging建模的高維

宋周洲 張涵寓 劉釗 朱平

摘要：

高維不確定性傳播目前面臨維度災(zāi)難和小樣本的問題，難以利用有限的樣本資源獲得高精度的分析結(jié)果，針對(duì)此問題，提出了一種基于監(jiān)督降維和自適應(yīng)Kriging模型的高維不確定性傳播方法。利用改進(jìn)充分降維方法將高維輸入投影到低維空間中，并利用Ladle估計(jì)器確定低維空間的維度。將降維投影矩陣嵌入Kriging核函數(shù)中以減少待估計(jì)超參數(shù)的數(shù)量，提高建模精度和效率。最后，創(chuàng)新性地定義了投影矩陣留一交叉驗(yàn)證誤差，并基于此提出了相應(yīng)的Kriging自適應(yīng)采樣策略，可以有效避免模型精度在自適應(yīng)采樣過程中發(fā)生較大波動(dòng)。數(shù)值算例與工程案例的結(jié)果表明，相比現(xiàn)有方法，所提方法能夠以較少的樣本點(diǎn)獲得高精度不確定性傳播結(jié)果，對(duì)復(fù)雜裝備結(jié)構(gòu)的不確定性分析和設(shè)計(jì)具有一定參考作用。

關(guān)鍵詞：高維不確定性傳播；監(jiān)督降維；Kriging模型；自適應(yīng)采樣

中圖分類號(hào)：O212;O302;TH122

DOI：10.3969/j.issn.1004132X.2024.05.001

開放科學(xué)（資源服務(wù)）標(biāo)識(shí)碼（OSID）：

A High-dimensional Uncertainty Propagation Method Based on Supervised

Dimension Reduction and Adaptive Kriging Modeling

SONG Zhouzhou1，2? ZHANG Hanyu1，2? LIU Zhao3? ZHU Ping1，2

1.State Key Laboratory of Mechanical System and Vibration，Shanghai Jiao Tong University，

Shanghai，200240

2.National Engineering Research Center of Automotive Power and Intelligent Control，

Shanghai Jiao Tong University，Shanghai，200240

3.School of Design，Shanghai Jiao Tong University，Shanghai，200240

Abstract： High-dimensional uncertainty propagation currently faced the curse of dimensionality， which made it difficult to utilize the limited sampling resources to obtain high-precision uncertainty analysis results. To address this problem， a high-dimensional uncertainty propagation method was proposed based on supervised dimension reduction and adaptive Kriging modeling. The high-dimensional inputs were projected into the low-dimensional space using the improved sufficient dimension reduction method， and the dimensionality of the low-dimensional space was determined by using the Ladle estimator. The projection matrix was embedded into the Kriging kernel function to reduce the number of hyperparameters to be estimated and improve the modeling accuracy and efficiency. Finally， the leave-one-out cross-validation error of the projection matrix was innovatively defined and the corresponding Kriging adaptive sampling strategy was proposed， which might effectively avoid large fluctuations of model accuracy in the adaptive sampling processes. The results of numerical and engineering examples show that， compared with the existing methods， the proposed method may obtain high-precision uncertainty propagation results with fewer sample points， which may provide references for the uncertainty analysis and design of complex structures.

Key words： high-dimensional uncertainty propagation; supervised dimension reduction; Kriging model; adaptive sampling

收稿日期：20231222

基金項(xiàng)目：國家自然科學(xué)基金（52375256，12302152）；上海市自然科學(xué)基金（23ZR1431600）

0? 引言

在實(shí)際中，結(jié)構(gòu)或系統(tǒng)的性能會(huì)受到諸如材料屬性波動(dòng)、加工誤差和隨機(jī)外載荷等眾多不確定性因素［1-3］的影響而表現(xiàn)出一定的隨機(jī)性。不確定性傳播旨在量化輸入不確定性如何通過系統(tǒng)傳遞到輸出［4］，是不確定分析與優(yōu)化的基礎(chǔ)。為了降低計(jì)算量，節(jié)省時(shí)間，通常采用代理模型來代替高保真度的仿真模型進(jìn)行不確定性傳播。常用的代理模型有徑向基函數(shù)［5］、支持向量機(jī)［6］、Kriging模型［7］、混沌多項(xiàng)式展開［8］、神經(jīng)網(wǎng)絡(luò)［9］等。然而，大多數(shù)代理模型都存在“維度災(zāi)難”的問題，也即構(gòu)建高精度代理模型所需的訓(xùn)練樣本數(shù)量隨著問題維度的增加而呈指數(shù)級(jí)增長［10］。對(duì)于需要昂貴的物理實(shí)驗(yàn)或高保真度仿真來獲取輸出響應(yīng)的復(fù)雜裝備結(jié)構(gòu)，通常難以獲得大量的訓(xùn)練樣本點(diǎn)，也就限制了代理模型乃至不確定性傳播的精度。

基于監(jiān)督降維和自適應(yīng)Kriging建模的高維不確定性傳播方法研究——宋周洲? 張涵寓? 劉? 釗等

中國機(jī)械工程 第35卷 第5期 2024年5月

近年來，基于代理模型的高維不確定性傳播方法受到了國內(nèi)外學(xué)者的廣泛關(guān)注。劉竟飛［11］提出一種基于貝葉斯深度學(xué)習(xí)的神經(jīng)網(wǎng)絡(luò)方法，能處理輸入變量服從任意概率分布的高維不確定性傳播問題。ZHANG等［12］提出了一種徑向基函數(shù)網(wǎng)絡(luò)模型，并將其用于機(jī)器人定位精度分析。赫萬鑫［13］提出了基于貝葉斯理論的自適應(yīng)稀疏多項(xiàng)式維度分解代理模型并用于水下航行體結(jié)構(gòu)動(dòng)力響應(yīng)不確定性分析。Kriging模型適用性強(qiáng)且能提供代理模型不確定性的度量，基于降維的高維Kriging建模方法得到了較多的研究。LIU等［14］首先利用Sammon映射進(jìn)行降維，然后構(gòu)建低維輸入到輸出的Kriging模型。Sammon映射屬于非監(jiān)督降維，由于未能考慮到輸出的信息，所構(gòu)建的代理模型精度可能不高［15］，因此，一些學(xué)者將監(jiān)督降維引入代理模型建模中。BOUHLEL等［16］利用偏最小二乘回歸進(jìn)行降維，并且將降維投影矩陣的元素嵌入到Kriging核函數(shù)中，相比直接在低維空間中構(gòu)建的Kriging模型具有更高的精度。CONSTANTINE等［17］提出了活躍子空間降維方法并將其用于構(gòu)建Kriging模型。ZHOU等［15］利用切片逆回歸方法進(jìn)行監(jiān)督降維，并在訓(xùn)練Kriging模型時(shí)利用降維矩陣構(gòu)建了新的核函數(shù)。SONG等［18］提出一種基于改進(jìn)充分降維的Kriging模型建模方法并用于高維不確定性傳播。

雖然高維Kriging模型建模方法得到了較為廣泛的研究，但現(xiàn)有方法通?；谝淮涡猿闃?，即先利用試驗(yàn)設(shè)計(jì)進(jìn)行抽樣，然后利用樣本點(diǎn)訓(xùn)練代理模型。一次性抽樣無法自適應(yīng)地調(diào)整樣本在空間中的分布，所需要的樣本數(shù)量通常偏大，這會(huì)浪費(fèi)一些昂貴的樣本資源。而自適應(yīng)抽樣方法能根據(jù)當(dāng)前代理模型提供的預(yù)測誤差信息自適應(yīng)地選擇最能提升模型精度的樣本，從而減少訓(xùn)練樣本數(shù)量，提高計(jì)算效率，因此，如何將自適應(yīng)抽樣方法與監(jiān)督降維方法有效結(jié)合，用盡可能少的樣本構(gòu)建出高精度的代理模型，對(duì)高維不確定性傳播具有重要意義。

本文針對(duì)復(fù)雜裝備結(jié)構(gòu)不確定性傳播面臨的維度災(zāi)難以及小樣本問題，研究了基于監(jiān)督降維的高維Kriging代理模型建模方法，提出了基于降維投影矩陣估計(jì)誤差的Kriging自適應(yīng)采樣策略，建立了高維不確定性傳播分析框架，可高效獲取高精度不確定性傳播結(jié)果。

1? 基于Kriging模型的不確定性傳播

1.1? 不確定性傳播

假設(shè)結(jié)構(gòu)或系統(tǒng)的輸出性能為Y，由函數(shù)Y=g（X）確定，其中X為p維輸入。不確定性傳播的目標(biāo)是在給定輸入?yún)?shù)概率密度函數(shù)X～fX（x）（x為隨機(jī)變量X的一個(gè)具體取值）的情況下求得輸出Y的統(tǒng)計(jì)信息。常用的統(tǒng)計(jì)信息包括Y的均值［19］

μY=∫g（x）fX（x）dx（1）

和方差［19］

σ2Y=∫（g（x）-μY）2fX（x）dx（2）

以及概率密度函數(shù)Y～fY（y）（y為隨機(jī)變量Y的一個(gè)具體取值）

解析求解式（1）和式（2）中的積分是非常困難的，通常利用蒙特卡羅模擬（Monte Carlo simulation， MCS）方法來近似估計(jì)式中的積分。MCS方法首先根據(jù)fX（x）生成大量獨(dú)立同分布的輸入樣本xi（x1，x2，…，xNMCS，其中NMCS為輸入樣本數(shù)量），隨后利用計(jì)算模型獲得輸入樣本對(duì)應(yīng)的響應(yīng)yi（y1，y2，…，yNMCS），則μY和σ2Y可分別通過下式進(jìn)行估計(jì)［4］：

μ^Y=1NMCS∑NMCSi=1yi（3）

σ^2Y=1NMCS∑NMCSi=1（yi-μ^Y）2（4）

然而，MCS方法需要大量樣本才能準(zhǔn)確估計(jì)出Y的不確定性信息，這對(duì)復(fù)雜裝備結(jié)構(gòu)而言是不可行的，因此，通常構(gòu)建系統(tǒng)輸出函數(shù)g（X）的代理模型g^（X），并用計(jì)算效率高的g^（X）代替g（X）進(jìn)行不確定性傳播。

1.2? Kriging模型

Kriging模型強(qiáng)大的表征能力以及可以給出代理模型不確定性的量化表達(dá)，使其成為了最常用的代理模型方法之一。普通Kriging模型［20］可以由下式表達(dá)：

y（x）=μ+Z（x）（5）

其中，μ為全局均值，Z（x）是一個(gè)表示局部變化的零均值高斯過程，其協(xié)方差函數(shù)可以寫為

Cov（Z（x），Z（x′））=σ2R（x，x′）（6）

式中，σ2為高斯過程的方差；R（·，·）表示相關(guān)性函數(shù)（也被稱作Kriging核函數(shù)）；x′為與x不同的另一個(gè)輸入矢量。

高斯相關(guān)性函數(shù)是使用最多的一類核函數(shù)，其形式為

R（x，x′）=exp（-∑pm=1θm（xm-x′m）2）（7）

式中，θm為Kriging模型待估計(jì)的第m個(gè)超參數(shù)；xm、x′m分別為x和x′的第m個(gè)分量。

給定訓(xùn)練集D={（xi，yi），i=1，2，…，N}，μ和σ2可通過極大似然估計(jì)方法估計(jì)如下：

μ^=1TR-1y1TR-11

σ^2=（y-μ^1）TR-1（y-μ^1）N（8）

其中，N為構(gòu)建Kriging模型的樣本點(diǎn)的數(shù)量；y為輸出向量，y=（y1，y2，…，yN）T；1為一個(gè)長度為N、元素均為1的向量；R為相關(guān)性矩陣，其元素為R（xi，xj）（i，j=1，2，…，N）。Kriging超參數(shù)θ={θ1，θ2，…，θp}可通過最大化中心對(duì)數(shù)似然函數(shù)獲得：

θ^=argmaxθ（-N2ln（σ^2）-12ln|R|）（9）

估計(jì)出Kriging超參數(shù)后，Kriging模型在未知點(diǎn)處的預(yù)測均值和預(yù)測方差可分別由下式獲得：

μ^y^（x）=μ^+rTR-1（y-μ^1）（10）

σ^2y^（x）=σ^2［1-rTR-1r+（1TR-1r-1）21TR-11］（11）

r=（R（x，x1），R（x，x2），…，R（x，xN））T

2? 基于監(jiān)督降維和自適應(yīng)采樣的高維不確定性傳播

2.1? 監(jiān)督降維

監(jiān)督降維的目標(biāo)是求得p×d維的降維投影矩陣W，將原p維的高維輸入X投影成d維的低維輸入WTX。 監(jiān)督降維主要分為兩步：①確定低維投影的方向，也就是如何求出W的列向量；②確定低維投影的數(shù)量，也就是如何確定低維空間的維度d。

確定低維投影方向的方法有偏最小二乘回歸［21］（partial least squares， PLS）、活躍子空間［17］（active subspace， AS）、切片逆回歸［22］（sliced inverse regression， SIR）、改進(jìn)充分降維［18］（improved sufficient dimension reduction， ISDR）等。PLS方法可能不適用于小樣本或輸入變量是非高斯的情形［23］。而AS方法則需要響應(yīng)的梯度信息，這在實(shí)際應(yīng)用中獲取難度較大。SIR方法需要對(duì)輸出數(shù)據(jù)進(jìn)行切片，通常難以確定最佳切片數(shù)量，并且由于自適應(yīng)采樣過程中樣本數(shù)量會(huì)動(dòng)態(tài)變化，導(dǎo)致SIR方法的切片數(shù)量也需要?jiǎng)討B(tài)調(diào)整，實(shí)施起來較為繁瑣。而ISDR方法不需要梯度信息，也不需要對(duì)數(shù)據(jù)進(jìn)行切片，因此本文選擇ISDR方法對(duì)高維數(shù)據(jù)進(jìn)行監(jiān)督降維。

ISDR屬于線性充分降維方法的一種，旨在尋找高維輸入的一個(gè)低維線性組合，該線性組合包含預(yù)測輸出所需要的全部信息。這在數(shù)學(xué)上可以描述為

YX|WTX（12）

其中，W為待估計(jì)的p×d維的投影矩陣，它將p維的高維輸入X投影到d維空間中；表示隨機(jī)變量的獨(dú)立性。式（12）的含義是，在給定低維輸入WTX的條件下，原高維輸入X與輸出Y相互獨(dú)立，也就是低維投影WTX包含了預(yù)測Y所需要的全部信息。在ISDR方法中，W可以按如下步驟計(jì)算得到［18］。

（1）對(duì)輸入數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化：

Z=Σ-1/2XX（X-E［X］）（13）

式中，ΣXX為輸入數(shù)據(jù)X的協(xié)方差矩陣；E［X］為X的數(shù)學(xué)期望；Z為標(biāo)準(zhǔn)化之后的輸入數(shù)據(jù)。

（2）對(duì)矩陣M=-E［Z′ZT|Y-Y′|］進(jìn)行估計(jì)，則有

M^=-1N（N-1）∑1≤i≠j≤NzizTj|yi-yj|（14）

（3）求出矩陣M^中前d大的特征值λ1≥λ2≥…≥λd所對(duì)應(yīng)的特征向量ξ1，ξ2，…，ξd。

（4）求出投影矩陣：

W=Σ-1/2XX（ξ1，ξ2，…，ξd）（15）

監(jiān)督降維的另一個(gè)關(guān)鍵問題是如何確定低維空間的維度d。常用的方法為：當(dāng)矩陣M中前d大的特征值之和占全部特征值之和的比例超過一閾值時(shí)（如99%），或當(dāng)矩陣M的特征值圖譜在第d個(gè)特征值處有顯著變化時(shí)，將維度確定為d。但是這類方法帶有一定的主觀性，缺乏數(shù)學(xué)上的嚴(yán)謹(jǐn)性。另一類方法則是從1開始逐漸增大低維空間的維度，對(duì)每一個(gè)維度構(gòu)建相應(yīng)的代理模型，并且用交叉驗(yàn)證求出代理模型的預(yù)測誤差，直至預(yù)測誤差收斂。但是這類方法需要反復(fù)構(gòu)建代理模型，計(jì)算效率低。注意到低維空間維度的估計(jì)問題等價(jià)于確定矩陣M的秩，而矩陣的秩可由其特征值以及特征向量的變化情況準(zhǔn)確估計(jì)得到，因此本文采用LUO等［24］提出的秩估計(jì)器——

Ladle估計(jì)器來估計(jì)矩陣M的秩。Ladle估計(jì)器計(jì)算方便，具有嚴(yán)格的數(shù)學(xué)定義，且LUO等［24］證明了Ladle估計(jì)器在大樣本下的收斂性質(zhì)，因此本研究用它來確定低維空間的維度d。有關(guān)Ladle估計(jì)器算法的實(shí)施細(xì)節(jié)可以參考文獻(xiàn)［24］。

2.2? 基于監(jiān)督降維的Kriging建模方法

估計(jì)出降維投影矩陣W后，即可用其構(gòu)建代理模型。由于直接將低維輸入WTX作為Kriging模型的輸入可能會(huì)造成信息損失從而導(dǎo)致模型預(yù)測精度偏低，因此將W的元素嵌入到Kriging相關(guān)性函數(shù)中構(gòu)建一個(gè)新的核函數(shù)［18］：

RW（x，x′）=exp（-∑dl=1∑pm=1θlw2ml（xm-x′m）2）（16）

其中，wml為位于矩陣W第m行、第l列的元素，θl為第l個(gè)新Kriging核函數(shù)的超參數(shù)?？芍虢稻S矩陣后，Kriging超參數(shù)的數(shù)量由式（7）中的p個(gè)減少到了式（16）中的d個(gè)，由于dp，因此利用改進(jìn)后的核函數(shù)訓(xùn)練Kriging模型會(huì)有更高的計(jì)算效率。

超參數(shù)θl同樣也可以通過最大化中心對(duì)數(shù)似然函數(shù)獲得。在進(jìn)行超參數(shù)優(yōu)化時(shí)，式（16）中θl的取值范圍是（0，+∞）。為了提高優(yōu)化性能，可以將θl進(jìn)行重參數(shù)化［25］：

θl=10φl? l=1，2，…，d（17）

重參數(shù)化后的超參數(shù)φl的取值范圍為全體實(shí)數(shù)。本文利用有限內(nèi)存Broyden-Fletcher-Goldfarb-

Shanno（BFGS）優(yōu)化算法［26］來對(duì)φl進(jìn)行優(yōu)化。

2.3? 自適應(yīng)采樣策略

有了高維Kriging模型建模方法后，就可將其與自適應(yīng)采樣相結(jié)合來進(jìn)一步減少訓(xùn)練模型所需的樣本量。目前已有許多面向全局Kriging建模的自適應(yīng)采樣策略，但是面向高維Kriging模型的自適應(yīng)采樣方法還不多見。由于基于監(jiān)督降維的Kriging建模方法引入了降維矩陣W，因此W的估計(jì)精度會(huì)直接影響到Kriging模型的精度。直觀上看，可以選取對(duì)W估計(jì)精度影響最大的樣本點(diǎn)作為新增樣本點(diǎn)。類似于代理模型的留一交叉驗(yàn)證誤差，可以定義投影矩陣在訓(xùn)練集中每一個(gè)樣本點(diǎn)處的留一估計(jì)誤差為

e（W）LOO（xi）=‖W-W-i‖F(xiàn)? i=1，2，…，N（18）

其中，W為利用訓(xùn)練集D={（xi，yi），i=1，2，…，N}中所有樣本點(diǎn)估計(jì)得到的投影矩陣，W-i為將樣本點(diǎn)（xi，yi）從D中剔除后估計(jì)得到的投影矩陣，‖·‖F(xiàn)表示矩陣的Frobenius范數(shù)。e（W）LOO（xi）越大，表明xi對(duì)投影矩陣估計(jì)精度的影響越大，則在xi附近選擇新增樣本點(diǎn)就有更大可能提高W的估計(jì)精度，因此可選擇在具有最大e（W）LOO的樣本點(diǎn)x*附近進(jìn)行采樣，x*由下式確定：

x*=argmaxx∈{x1，x2，…，xN}e（W）LOO（xi）（19）

這里參照LIU等［27］的工作，利用歐氏距離來定義“鄰近”關(guān)系，即如果未知點(diǎn)x到x*的歐氏距離小于x到訓(xùn)練集中其他任意點(diǎn)的歐氏距離，則認(rèn)為x是x*附近的點(diǎn)。

由于輸入空間中位于x*附近的點(diǎn)有無窮多個(gè)，故為了簡化采樣流程，通常根據(jù)輸入?yún)?shù)的不確定性分布，利用MCS方法獲得一個(gè)大樣本集S（樣本數(shù)量通常為105及以上），并用該樣本集中的點(diǎn)來代替整個(gè)輸入空間，因此，在確定具有最大e（W）LOO的x*后，直接從大樣本集S中找到x*附近的點(diǎn)，記為S*：

S*={x∈S|‖x-x*‖2≤‖x-xi‖2，xi≠x*}（20）

隨后從S*中選擇具有最大Kriging預(yù)測方差的點(diǎn)作為新增樣本點(diǎn)，其表達(dá)式如下：

xnew=argmaxx∈S*（σ^2y^（x））（21）

值得指出的是，之所以這里從S*而不是整個(gè)輸入空間S中選擇具有最大預(yù)測方差的點(diǎn)，是因?yàn)檫@可以保證在相鄰兩次自適應(yīng)采樣中，投影矩陣W不至于發(fā)生很大的波動(dòng)進(jìn)而導(dǎo)致Kriging模型精度產(chǎn)生大幅變動(dòng)，本文第3節(jié)將結(jié)合算例結(jié)果對(duì)此問題進(jìn)行更深入的討論。

2.4? 算法流程

本文所提出的基于監(jiān)督降維和自適應(yīng)采樣的高維不確定性傳播流程如圖1所示，其主要步驟如下：

（1）利用MCS方法獲得輸入空間的大樣本集S，利用最優(yōu)拉丁超立方抽樣獲得初始輸入樣本點(diǎn)并利用高保真度有限元仿真獲得相應(yīng)的輸出，形成訓(xùn)練集D。

（2）對(duì)訓(xùn)練集D中的數(shù)據(jù)進(jìn)行改進(jìn)充分降維，獲得降維投影矩陣W，將W嵌入Kriging核函數(shù)后訓(xùn)練得到高維Kriging代理模型。

（3）判斷新增樣本點(diǎn)是否達(dá)到規(guī)定值，若是，則跳轉(zhuǎn)至步驟（5），若不是，則跳轉(zhuǎn)至步驟（4）。

（4）進(jìn)行自適應(yīng)采樣。首先根據(jù)式（19）確定具有最大投影矩陣估計(jì)誤差的樣本點(diǎn)x*；再根據(jù)式（20）從S中選擇出x*附近的樣本點(diǎn)形成集合S*；隨后根據(jù)式（21）確定新增樣本點(diǎn)xnew；最后利用高保真度仿真模型獲取新樣本點(diǎn)對(duì)應(yīng)的響應(yīng)ynew，并更新訓(xùn)練集D=D∪（xnew，ynew），且跳轉(zhuǎn)至步驟（2）。

（5）結(jié)束自適應(yīng)采樣，并獲得高維Kriging模型。

（6）利用構(gòu)建好的高維Kriging模型進(jìn)行不確定性傳播。

3? 算例與應(yīng)用研究

3.1? 導(dǎo)熱問題不確定性傳播

本算例是一個(gè)二維矩形區(qū)域中的熱傳導(dǎo)問題［28］。如圖2所示，所考察的區(qū)域?yàn)镃=（-0.5，0.5）m×（-0.5，0.5）m。區(qū)域D的上邊界溫度恒定為0 ℃，并且C的初始溫度也為0 ℃。C的左、下、右邊界分別有恒定的熱流H1、H2、H3流入（熱流密度的符號(hào)以流入?yún)^(qū)域C為正）。C內(nèi)有一熱源A=（0.2，0.3）m×（0.2，0.3）m以恒定功率Q向外釋放熱量。受材料分布的場不確定性的影響，區(qū)域C的熱導(dǎo)率κ（z）用一個(gè)對(duì)數(shù)正態(tài)隨機(jī)場來描述，其表達(dá)式如下：

κ（z）=exp（aκ+bκg（z））（22）

其中，aκ、bκ為參數(shù)，z為區(qū)域中某點(diǎn)的坐標(biāo)，g（z）表示一個(gè)標(biāo)準(zhǔn)高斯隨機(jī)場，其自相關(guān)函數(shù)為

ρ（z1，z2）=exp（-‖z1-z2‖20.04）（23）

式（22）中aκ、bκ的取值使κ（z）的均值和標(biāo)準(zhǔn)差分別為μκ=1 W/（m·K）以及σκ=0.3 W/（m·K）。本算例所考察的輸出為區(qū)域B=（-0.2，-0.3）m×（-0.2，-0.3）m在熱流傳導(dǎo)1 s后的平均溫度，其表達(dá)式如下：

Y=1|B|∫z∈BT（z，t=1）dz（24）

其中，T（z，t）表示坐標(biāo)為z的點(diǎn)在t時(shí)刻的溫度，Y由有限元分析方法求出。表1給出了導(dǎo)熱問題隨機(jī)參數(shù)的信息。隨機(jī)場g（z）利用最優(yōu)線性估計(jì)展開［29］離散為53個(gè)隨機(jī)變量，即

g（z）=∑53i=1ξiλiφTicg（z）（25）

以確保截?cái)嗾`差小于1%。其中，ξi為相互獨(dú)立的標(biāo)準(zhǔn)正態(tài)變量，λi、φi分別為隨機(jī)場協(xié)方差矩陣的特征值和特征向量，cg（z）為協(xié)方差向量。在本算例中，由于隨機(jī)場κ（z）被離散成53個(gè)隨機(jī)變量，再加上表1中的4個(gè)隨機(jī)參數(shù)，本問題是一個(gè)57維的高維問題。

首先比較了利用不同方法構(gòu)建的代理模型建模誤差，這里利用歸一化的均方根誤差（normalized root mean square error， NRMSE）來量化模型誤差：

VNRMSE=1Ntest∑Ntesti=1（yi-y^iyi）2（26）

式中，Ntest為測試集樣本數(shù)量;yi為響應(yīng)的真實(shí)值;y^i為代理模型的預(yù)測值。

利用MCS方法額外生成1000個(gè)樣本作為測試集。所有方法運(yùn)行10次后的箱線圖見圖3，其中OK代表普通Kriging模型；A-OK代表自適應(yīng)普通Kriging建模方法，自適應(yīng)采樣策略為最大化期望預(yù)測誤差［27］；KISDR代表基于改進(jìn)充分降維的Kriging建模方法；A-KISDR代表基于最大化期望預(yù)測誤差［27］采樣準(zhǔn)則的改進(jìn)充分降維Kriging建模方法；A-KISDR-W則是本文所提出的方法。OK和KISDR方法的訓(xùn)練集的樣本數(shù)量為150，A-OK、A-KISDR和A-KISDR-W方法的初始訓(xùn)練集樣本數(shù)量為100，隨后利用自適應(yīng)采樣的新增樣本數(shù)量為50，以保證最后所有方法的訓(xùn)練樣本數(shù)量一致。從圖3中可以看出，OK方法和A-OK方法的建模精度相差不大，這說明對(duì)于高維問題，如果不引入降維技術(shù)，則自適應(yīng)采樣相比一次性采樣并沒有明顯優(yōu)勢。對(duì)于KISDR、A-KISDR和A-KISDR-W這三種基于降維的方法，自適應(yīng)采樣方法（A-KISDR-W和A-KISDR）的建模精度和穩(wěn)定性優(yōu)于一次性采樣方法（KISDR），并且基于投影矩陣估計(jì)誤差方法（A-KISDR-W）的建模精度比傳統(tǒng)自適應(yīng)采樣策略（A-KISDR）的建模精度更高，這說明了本文所提方法的優(yōu)勢。

圖4給出了A-OK、A-KISDR和A-KISDR-W這三種自適應(yīng)采樣方法預(yù)測誤差的迭代曲線，可以看出，A-KISDR方法相比于A-OK方法，誤差下降速度總體上更快。相比于A-KISDR方法，A-KISDR-W方法在迭代過程中的波動(dòng)更小。這是因?yàn)锳-KISDR-W方法每次選擇在投影矩陣估計(jì)誤差最大的區(qū)域內(nèi)進(jìn)行采樣，能夠有效防止投影矩陣出現(xiàn)較大波動(dòng)，從而避免代理模型精度出現(xiàn)大的波動(dòng)。

圖5給出了不同方法預(yù)測得到的區(qū)域B平均溫度的概率密度函數(shù)，概率密度函數(shù)是根據(jù)MCS方法生成的105個(gè)樣本點(diǎn)利用核密度估計(jì)得到的。由圖5可以看出，本文提出的A-KISDR-W方法估計(jì)得到的概率密度曲線與真實(shí)的概率密度曲線最為接近，這驗(yàn)證了所提方法在高維不確定性傳播上的優(yōu)勢。

3.2? 異種材料膠鉚混合連接結(jié)構(gòu)不確定性傳播

本算例考察的是碳纖維增強(qiáng)復(fù)合材料與鋁合金膠鉚混合連接結(jié)構(gòu)準(zhǔn)靜態(tài)拉伸強(qiáng)度的不確定性傳播。異種材料膠鉚混合連接結(jié)構(gòu)結(jié)合了兩種優(yōu)異輕質(zhì)高強(qiáng)材料，通過膠接、鉆孔、鉚接的工序制成，在新能源汽車輕量化中具有良好的應(yīng)用前景。由于汽車上大部分的結(jié)構(gòu)失效發(fā)生在連接位置，因此連接性能對(duì)車身強(qiáng)度、彎扭剛度、疲勞耐久等性能有重大影響。圖6給出了膠鉚混合連接結(jié)構(gòu)的示意圖。

由于膠鉚連接結(jié)構(gòu)屬于多材料多工藝連接結(jié)構(gòu)，影響其力學(xué)性能的參數(shù)也非常多，如圖7所示，既包括四種材料（碳纖維復(fù)合材料、膠、不銹鋼鉚釘、鋁合金）的力學(xué)性能參數(shù)，又包括膠接、鉆孔、鉚接等工藝參數(shù)，以及結(jié)構(gòu)的幾何參數(shù)等，因此，膠鉚連接結(jié)構(gòu)性能的不確定性傳播屬于典型的高維問題。本文主要考察連接結(jié)構(gòu)的30個(gè)不確定性參數(shù)對(duì)結(jié)構(gòu)的準(zhǔn)靜態(tài)拉伸強(qiáng)度不確定性的影響，30個(gè)參數(shù)的不確定性分布信息可參考ZHANG等［30］的工作。這里利用經(jīng)過試驗(yàn)驗(yàn)證的高保真度有限元仿真模型［31］來獲得結(jié)構(gòu)的強(qiáng)度，由于仿真模型非常耗時(shí)，因此需要利用本文所提方法構(gòu)建出結(jié)構(gòu)強(qiáng)度的代理模型。

表2列出了采用不同方法所構(gòu)建的代理模型的預(yù)測誤差，其中OK和KISDR方法的訓(xùn)練集樣本數(shù)量為150，A-OK、A-KISDR和A-KISDR-W方法的初始訓(xùn)練集樣本數(shù)量為100，隨后利用自適應(yīng)采樣新增50個(gè)樣本。利用拉丁超立方抽樣額外生成100個(gè)樣本作為測試集。從表2中可以看出，基于降維的方法的建模精度高于普通Kriging方法的建模精度，自適應(yīng)采樣方法的建模精度高于一次性采樣方法的建模精度。同時(shí)，本文所提出的采樣策略相比傳統(tǒng)的自適應(yīng)采樣策略具有更高的建模精度。

圖8給出了A-OK、A-KISDR和A-KISDR-W這三種自適應(yīng)采樣方法預(yù)測誤差的迭代曲線，可以看出，A-OK方法在此問題中精度波動(dòng)非常大，這可能是由于普通Kriging方法在樣本量較小時(shí)建模精度不穩(wěn)定導(dǎo)致的。而A-KISDR和A-KISDR-W方法相對(duì)來說波動(dòng)較小，也說明了基于降維的Kriging模型在小樣本時(shí)也可以獲得較為精確和穩(wěn)定的代理模型。相比于A-KISDR方法，A-KISDR-W方法在迭代過程中的波動(dòng)更小。這說明所提方法能夠有效防止投影矩陣出現(xiàn)較大波動(dòng)，從而避免代理模型精度出現(xiàn)大的波動(dòng)。圖9給出了不同方法預(yù)測得到的膠鉚連接結(jié)構(gòu)強(qiáng)度的概率密度函數(shù)，概率密度函數(shù)是根據(jù)MCS方法生成的105個(gè)樣本點(diǎn)利用核密度估計(jì)得到的。雖然不同方法預(yù)測得到的概率密度函數(shù)有細(xì)微區(qū)別，但都近似服從于正態(tài)分布。

4? 結(jié)論

本文針對(duì)基于代理模型的高維不確定性傳播中面臨的維度災(zāi)難問題，提出了一種基于改進(jìn)充分降維和自適應(yīng)Kriging建模的高維不確定性傳播方法，可以較少的樣本量獲得高精度不確定性傳播結(jié)果。所獲得的結(jié)論如下：

（1）利用改進(jìn)充分降維對(duì)高維輸入進(jìn)行降維，不需要響應(yīng)的梯度信息，也不需要對(duì)數(shù)據(jù)進(jìn)行切片，簡單易行，適用范圍廣。

（2）利用具有嚴(yán)格數(shù)學(xué)定義的Ladle估計(jì)器對(duì)低維輸入空間的維度進(jìn)行估計(jì)，既提高了維度估計(jì)的效率，又能確保結(jié)果穩(wěn)定可靠。

（3）本文創(chuàng)新性地提出了基于投影矩陣留一估計(jì)誤差的自適應(yīng)采樣策略，數(shù)值算例和工程應(yīng)用結(jié)果表明，所提出的策略相比傳統(tǒng)的Kriging自適應(yīng)采樣策略在迭代過程中精度波動(dòng)更小，并且在新增樣本數(shù)量相同時(shí)，能獲得更高精度的Kriging模型。

所提方法對(duì)復(fù)雜裝備結(jié)構(gòu)的不確定性分析和設(shè)計(jì)具有一定參考。

參考文獻(xiàn)：

［1］? 許燦. 面向復(fù)雜產(chǎn)品穩(wěn)健可靠性的多層次系統(tǒng)不確定性分析與優(yōu)化設(shè)計(jì)方法研究［D］. 上海：上海交通大學(xué)， 2021.

XU Can. Research on Uncertainty Analysis and Uncertainty-based Design Optimization of Multilevel Systems for Robustness and Reliability of Complex Products［D］. Shanghai：Shanghai Jiao Tong University， 2021.

［2］? 蔣琛， 邱浩波， 高亮. 隨機(jī)不確定性下的可靠性設(shè)計(jì)優(yōu)化研究進(jìn)展［J］. 中國機(jī)械工程， 2020， 31（2）： 190-205.

JIANG Chen， QIU Haobo， GAO Liang. Research Progresses in Reliability-based Design Optimization under Aleatory Uncertainties［J］. China Mechanical Engineering， 2020， 31（2）： 190-205.

［3］? 高進(jìn)， 崔海冰， 樊濤，等. 一種基于自適應(yīng)Kriging集成模型的結(jié)構(gòu)可靠性分析方法［J］. 中國機(jī)械工程， 2024， 35（1）： 83-92.

GAO Jin， CUI Haibing， FAN Tao， et al. A Structural Reliability Calculation Method Based on Adaptive Kriging Ensemble Model［J］. China Mechanical Engineering， 2024， 35（1）： 83-92.

［4］? LEE S H， CHEN W. A Comparative Study of Uncertainty Propagation Methods for Black-box-type Problems［J］. Structural and Multidisciplinary Optimization， 2009， 37： 239-253.

［5］? GUTMANN H M. A Radial Basis Function Method for Global Optimization［J］. Journal of Global Optimization， 2001， 19： 201-227.

［6］? SMOLA A J， SCHLKOPF B. A Tutorial on Support Vector Regression［J］. Statistics and Computing， 2004， 14： 199-222.

［7］? RASMUSSEN C E， WILLIAMS C K. Gaussian Processes for Machine Learning［M］. Cambridge： MIT Press， 2006.

［8］? XIU D， KARNIADAKIS G E. The Wiener-Askey Polynomial Chaos for Stochastic Differential Equations［J］. SIAM Journal on Scientific Computing， 2002， 24： 619-644.

［9］? BISHOP C M. Neural Networks for Pattern Recognition［M］. New York： Oxford University Press， 1995.

［10］? SHAN S， WANG G G. Metamodeling for High Dimensional Simulation-based Design Problems［J］. Journal of Mechanical Design， 2010， 132： 051009.

［11］? 劉竟飛. 基于貝葉斯深度學(xué)習(xí)的高維不確定性傳播方法研究［D］. 長沙：湖南大學(xué)， 2021.

LIU Jingfei. Study of Uncertainty Propagation Method via Bayesian Deep Learning for High Dimensional Problems［D］. Changsha：Hunan University， 2021.

［12］? ZHANG D， ZHANG N， YE N， et al. Hybrid Learning Algorithm of Radial Basis Function Networks for Reliability Analysis［J］. IEEE Transactions on Reliability， 2021， 70： 887-900.

［13］? 赫萬鑫. 基于高維模型表征的隨機(jī)不確定性分析高效算法研究［D］. 大連：大連理工大學(xué)， 2021.

HE Wanxin. Efficient Stochastic Uncertainty Analysis Algorithms Based on High-dimensional Model Representation［D］. Dalian： Dalian University of Technology， 2021.

［14］? LIU B， ZHANG Q， GIELEN G G E. A Gaussian Process Surrogate Model Assisted Evolutionary Algorithm for Medium Scale Expensive Optimization Problems［J］. IEEE Transactions on Evolutionary Computation， 2014， 18： 180-192.

［15］? ZHOU Y， LU Z. An Enhanced Kriging Surrogate Modeling Technique for High-dimensional Problems［J］. Mechanical Systems and Signal Processing， 2020， 140： 106687.

［16］? BOUHLEL M A， BARTOLI N， OTSMANE A， et al. Improving Kriging Surrogates of High-dimensional Design Models by Partial Least Squares Dimension Reduction［J］. Structural and Multidisciplinary Optimization， 2016， 53： 935-952.

［17］? CONSTANTINE P G， DOW E， WANG Q. Active Subspace Methods in Theory and Practice： Applications to Kriging Surfaces［J］. SIAM Journal on Scientific Computing， 2014， 36： A1500-A1524.

［18］? SONG Z， LIU Z， ZHANG H， et al. An Improved Sufficient Dimension Reduction-based Kriging Modeling Method for High-dimensional Evaluation-expensive Problems［J］. Computer Methods in Applied Mechanics and Engineering， 2024， 418： 116544.

［19］? TRIPATHY R， BILIONIS I， GONZALEZ M. Gaussian Processes with Built-in Dimensionality Reduction： Applications to High-dimensional Uncertainty Propagation［J］. Journal of Computational Physics， 2016， 321：191-223.

［20］? KLEIJNEN J P C. Kriging Metamodeling in Simulation：a Review［J］. European Journal of Operational Research， 2009， 192：707-716.

［21］? HELLAND I S. On the Structure of Partial Least Squares Regression［J］. Communications in Statistics-simulation and Computation， 1988， 17：581-607.

［22］? LI K C. Sliced Inverse Regression for Dimension Reduction［J］. Journal of the American Statistical Association， 1991， 86：316-327.

［23］? GOODHUE D L， LEWIS W， THOMPSON R. Does PLS Have Advantages for Small Sample Size or Non-normal Data？［J］. MIS Quarterly， 2012， 36：981-1001.

［24］? LUO W， LI B. Combining Eigenvalues and Variation of Eigenvectors for Order Determination［J］. Biometrika， 2016， 103：875-887.

［25］? MACDONALD B， RANJAN P， CHIPMAN H. GPfit：an R Package for Fitting a Gaussian Process Model to Deterministic Simulator Outputs［J］. Journal of Statistical Software， 2015， 64：1-23.

［26］? LIU D C， NOCEDAL J. On the Limited Memory BFGS Method for Large Scale Optimization［J］. Mathematical Programming， 1989， 45：503-528.

［27］? LIU H， CAI J， ONG Y S. An Adaptive Sampling Approach for Kriging Metamodeling by Maximizing Expected Prediction Error［J］. Computers & Chemical Engineering， 2017， 106：171-182.

［28］? KONAKLI K， SUDRET B. Polynomial Meta-models with Canonical Low-rank Approximations：Numerical Insights and Comparison to Sparse Polynomial Chaos Expansions［J］. Journal of Computational Physics， 2016， 321：1144-1169.

［29］? LI C C， DER KIUREGHIAN A. Optimal Discretization of Random Fields［J］. Journal of Engineering Mechanics， 1993， 119：1136-1154.

［30］? ZHANG H， ZHANG L， XU C， et al. Global Sensitivity Analysis of Mechanical Properties in Hybrid Single Lap Aluminum-CFRP （Plain Woven） Joints Based on Uncertainty Quantification［J］. Composite Structures， 2022， 280：114841.

［31］? ZHANG H， ZHANG L， LIU Z， et al. Research in Failure Behaviors of Hybrid Single Lap Aluminum-CFRP（Plain Woven） Joints［J］. Thin-Walled Structures， 2021， 161：107488.

（編輯? 胡佳慧）

作者簡介：

宋周洲，男，1998年生，博士研究生。研究方向主要為不確定性分析、可靠性設(shè)計(jì)等。E-mail：zhouzsong@sjtu.edu.cn。

朱? 平（通信作者），男，1966年生，教授、博士研究生導(dǎo)師。研究方向主要為輕量化設(shè)計(jì)與制造、汽車安全性與可靠性、大數(shù)據(jù)與數(shù)字孿生等。E-mail：pzhu@sjtu.edu.cn。