999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

一種基于證據理論的主動學習可靠性分析方法

2025-07-16 00:00:00張哲寶文禮姚中洋
湖南大學學報·自然科學版 2025年6期
關鍵詞:方法模型

中圖分類號:TB114.3 文獻標志碼:A

Abstract:For the reliability analysis problem characterized by a single failure mode,cognitive uncertainty, and“black-box”models,anactive learning reliability analysis method basedon evidence theory is proposed.This method eficientlyandaccurately determines thecredibilityandverisimilitudeof structures.It handles cognitive uncertain variables using evidence theory,initiates initial training sample construction for a Kriging model,and combines optimization methods with activelearning to search for optimal training samples across theentire input variable space.This approach refines the Kriging model chronicallywith optimal training samples,replacing the functionalfunction with the Kriging model to predict unknown points for credibilityand verisimilitude calculation of the structure.By integrating optimization methods with active learning,the method relaxes constraints on candidate samplelocationsduring traditional training samplesearch,thereby identifying training samples that beter enhance the Kriging model’'s correction effects and improve the effciencyand succssrate of Kriging model construction. Numerical examples demonstrate the method’s computational efectiveness and itsapplication to the reliability analysis of vehicle frontal collisions.

Key words: structural reliability;reliability analysis;evidence theory;black box problem;active learning kriging model

在實際工程問題中,由于測量條件受限、制造精度不足和服役環境多變等因素的影響,會產生各種不確定性.這些不確定性在多數情況下表現為均值附近的小范圍波動,但是,結構運行的安全性可能因為多源不確定性的影響而發生波動,甚至失效.因此,采用恰當的方法量化和控制不確定性對于確保結構安全性至關重要.基于對信息的掌握情況,不確定性可以分為兩類[1.一類是隨機不確定性,基于系統內部的隨機過程或外部隨機事件的影響而產生,通常通過概率密度函數表征.另一類是認知不確定性,主要受當前認知水平和測量工具等限制,沒有足夠的信息表征不確定變量的分布類型.目前用于處理認知不確定性的理論有:證據理論[3]、區間分析理論4模糊集理論5和可能性理論等.相比之下,證據理論使用靈活的識別框架和基本可信度分配函數表征不確定信息,是當前用于解決認知不確定性的主要方法之一.

證據理論于1967年由Dempster提出,經過其學生Shafer進一步研究與發展,也被稱為Dempster-Shafer(D-S)理論7].由于認知不確定性的存在,證據理論使用可信度函數(belieffunction,Bel)和似真度函數(plausibilityfunction,Pl)組成的信度區間共同量化不確定性,不同于概率理論使用單一的失效概率描述隨機不確定性.證據理論使用基本可信度分配(basicprobabilityassignment,BPA)函數量化認知不確定性,該函數是定義在集合上的可信度分配函數.當信息量足夠多以至于能夠獲得概率密度函數時,證據理論近似等效于概率理論;反之,當信息量減少到只有上下界時,證據理論近似等效于區間分析理論8.然而,基本可信度分配函數是定義于集合上的非連續函數,使得基于證據理論的可靠性分析需要對不確定域內每個集合上的功能函數進行極值分析,因此,近年來,如何提升基于證據理論的可靠性分析的計算效率成為該領域的研究重點之一[9-10].為此,同行學者提出了不同的求解策略,包括:頂點法、焦元削減法、概率等價法和代理模型法等.

頂點法:該方法最早由Dong等]提出.該方法考慮到在實際工程問題中,不確定性是以設計值為中心的小范圍波動,將其分配到焦元內部后波動范圍再次縮小.所以該方法以極限狀態函數在焦元的微小范圍內,非線性影響很小的前提,借助求解焦元頂點位置的響應值實現焦元類型的判斷[12].頂點法相比于傳統方法計算效率更高,結果魯棒性強13.但是頂點法需要計算所有焦元頂點的值,面對維度較高的證據可靠性分析問題,仍然存在沉重的計算負擔.

焦元削減法:該方法通過非概率可靠性指標構造了不確定域內的一個輔助區域,落入該區域的焦元不需要計算功能函數的極值,從而降低計算成本[14].在證據可靠性分析中,很多工程問題的可靠度通常很高 (Pr?99% ),可以斷定大部分的焦元位于安全域,安全域內的焦元不再需要極值分析,從而大幅減少計算成本[15].Mourelatos等[16]將順序劃分策略用于焦元和可靠域之間關系的識別,提升了計算效率.

概率等價法:類似于概率理論中的最可能失效點,姜潮等[7提出針對證據理論可靠性分析的最可能失效焦元,該焦元對證據理論可信度與似真度的計算具有重要貢獻度.Xiao等i8依據證據變量中BPA的面積與隨機變量的PDF相等,提出兩者之間的轉換方法,能夠實現在證據空間中穩定且準確地搜索到最可能失效焦元.Zhang等[以最大可能失效焦元為中心,使用一階泰勒級數展開和二階泰勒級數展開,提出證據理論可靠性分析的一次二階矩法和二次二階矩法,高效地實現了結構可靠度計算.

代理模型法:實際工程的“黑箱\"問題,結構功能函數調用通常只能借助有限元仿真或者計算流體力學等方法實現.但是,對于某些復雜模型,如高鐵或者飛機,單次仿真需要幾天甚至幾個月時間.使用傳統證據可靠性分析方法解決此類問題,所需的計算成本難以接受.面對耗時且昂貴的實驗或者仿真計算,很多學者通過構建代理模型替代原來的實驗或者有限元仿真,大幅減少計算成本同時也有效地實現了可靠性分析的目的.Bae等2用多點逼近法通過在設計空間中選擇一系列局部近似點構建極限狀態函數的代理模型,通過代理模型取代了昂貴的結構仿真.Zhang等2借助實驗設計技術,搜索重要控制點,并基于控制點構建高精度的徑向基函數響應面,高效地完成了可靠性計算.曹亮等22借助支持向量回歸構建隱示功能函數的近似模型. Yin 等[23]借助均勻性方法為每個證據變量創建概率密度函數,使用多項式混沌展開來近似模擬證據變量變化范圍內的聲學系統響應,然后基于多項式混沌展開,通過數值求解器高效地完成所有焦元響應的上下界計算.

Kriging模型是一種高精度的代理模型,通過顯式高斯隨機場逼近實際工程中的“黑箱”模型,能夠提供未知點的最佳線性估計和預測誤差[24.25].為了提高Kriging模型的精度,學者們引入了主動學習過程,使得具有主動學習過程的Kriging模型在可靠性分析領域被廣泛應用[26.27].Yang等[28]提出了證據理論框架下的主動學習Kriging的可靠性分析方法,該方法使用主動學習Kriging模型為功能函數提供正確符號預測焦元以判斷焦元類型,此外為了提高求解焦元極值的準確率和效率,基于卡羅需-庫恩-塔克條件提出了KKTO優化方法.Zhang等29以 U 學習函數為基礎,根據基本變量的BPA將原證據空間進行細分,然后確定功能函數與細分得到空間的交點,通過交點構建Kriging模型,再用拉丁超立方抽樣法(Latinhypercubesampling,LHS)生成額外的樣本點,以提高初始Kriging模型的準確性.Yang等[30]提出一種新的學習函數UET(uncertaintyestimation tech-nique),使用識別功能函數的下限或上限預測符號錯誤的概率最大的焦元,通過實時監控預測誤差,及時終止學習過程;陳澤權[31]從快速滿足全局收斂條件的角度出發,基于Krigingbeliever準則及重要性抽樣原理,提出了自適應結構可靠性分析的快速收斂策略.

傳統基于Kriging模型的證據可靠性分析方法在構建Kriging模型過程中,通常抽樣一定數量的樣本作為候選樣本,然后從確定的候選樣本中搜索訓練樣本[30.32],在文中稱其為樣本點法,該方法搜索到最佳訓練樣本的概率較小,導致模型構建效率下降,甚至失效.因此提出了一種基于證據理論的主動學習可靠性分析方法,將優化方法和主動學習過程相結合,從整個輸人變量空間中自適應地搜索訓練樣本,降低了傳統方法對訓練樣本的約束,大幅提高搜索到最佳訓練樣本的概率,進而達到提升Kriging模型構建效率和成功率的效果.本文的主要內容為:第1章闡述了證據可靠性分析的基本理論;第2章介紹了主動學習Kriging模型;第3章闡述所提方法的算法構造;第4章通過算例驗證了本文方法的可行性;第5章為文本的結論和未來期望.

1證據可靠性分析基本理論

在實際的結構可靠性分析中,通常需要使結構滿足一定的功能特性,比如:強度、撓度、硬度等.所需要滿足的功能特性即功能函數,用 g(X) 表示.向量 X={X1,X2,…,Xn} 表示與結構有關的多維證據變量的集合,在文中表示證據變量.

根據設計的結構是否滿足某個功能特性,可以將不確定性輸人變量組成的空間劃分為安全域 G= {g|g(X)gt;0} 和失效域 ,如圖1所示.

圖1安全域與失效域示意圖Fig.1 Safetydomainand failuredomaindiagram

在證據理論框架下進行可靠性分析,首先,定義每個變量 X. 定義識別框架 Θ={x1,x2,…,xi,…,xn}. xi 為識別框架 Θ 的一個事件或元素, n 為元素個數,i=1,2,…,n. 識別框架中識別的含義是:對于一個命題,從與之相關的所有可能結果中區分并唯一區分出正確的解答.識別框架中的所有基本元素構成一個冪集 基本可信度分配 m 用于描述命題的可信任程度,類似于概率理論中的概率密度函數.設Θ 為識別框架,則基本可信度分配 ?m 定義為從集合 到[0,1]的映射函數,即 然后,對證據變量的每一個焦元定義基本可信度分配.當信息來源于多個專家或者系統時,則需要使用信息融合的方式對多條信息進行融合,求解每個焦元的基本可信度分配[33].

多數情況下,結構的功能特性會受到多個輸入變量的影響,需要建立多維證據變量的聯合識別框架和聯合基本可信度分配[8.當認知不確定性變量之間相互獨立,對于 N 維輸入變量,可以使用笛卡兒乘積的方式求解,如式(2):

式中: Θ 為聯合識別框架 {AX1,AX2,…,AXN 分別表示第1個,第2個,…,第 N 個證據變量的焦元; 分別表示第1個,第2個,…,第 N 個輸入變量的識別框架.聯合焦元的求解:

AX=AX?X?1×AX?2×…×AX?N

因為 AX1,AX2,…,AXN 都是區間,所以 AX 是一個 N 維的空間.聯合基本可信度分配的求解:

獲得證據變量的聯合基本可信度分配以及安全域之后,便可以對結構進行可信度和似真度的求解[34].證據理論框架下,結構可靠度的求解公式如式(5):

根據式(5),在求解可信度和似真度之前要先確定焦元 AX 與安全域 G 的位置關系,表明焦元 AX 完全還是部分位于安全域內部,或者 AX 在安全域外.為了能夠準確地判斷焦元 Ax 與安全域的位置關系,需要先求解功能函數 g(X) 在每一個焦元 AX 上的極值.

目前求解極值的方法主要有頂點法和優化算法,鑒于大部分工程問題中,結構的極限狀態函數通常表現出較強的非線性,為了減小非線性對可靠度計算結果的影響,通常采用SQP(sequential quadraticprogramming)優化算法對每一個 ??AX 求解極值.

求解出每一個焦元的極大值和極小值后,根據Ax 的極值分析 Ax 與極限狀態面的位置關系.如圖2(a),以焦元 AX 為例, .Ax 的極小值大于0,表明 AX 全部位于安全域內,在計算安全性時, .Ax 需要同時計入 Bel(G) 和 Pl(G) 中;圖2(b)中, ?AX 的極大值小于0,表明 AX 完全位于失效域內,在計算安全性時, .Ax 不計入 Bel(G) 和 Pl(G) 中;圖2(c)中, .Ax 的極大值大于0,極小值小于0,由此可以推斷, .Ax 部分位于安全域內,部分位于失效域內,屬于邊界焦元, AX 只計入 Pl(G) 不計入

Bel(G)

通過以上方法,逐個完成對所有焦元的分類,然后代人式(5)計算,結合 Bel(G)?P(G)?Pl(G) ,得到結構可靠度 P(G) 的范圍[21].在極值的分析過程中往往會產生大量的計算成本,隨著維度的增加,焦元數量指數級增長,給可靠性分析帶來巨大挑戰.

g(X)lt;0g(X)gt;0 g(Xlt;0 Ax AX Ax \ g(Xgt;0 g(X)lt;0 g(X)gt;0 g(X)=0g(X)=0 g(X)=0gmin(X)lt;0gmin(X)gt;0 gmax(X)lt;0 gmax(X)gt;0(a)可靠焦元 (b)失效焦元 (c)邊界焦元圖2焦元與極限狀態函數的可能位置關系

2主動學習Kriging模型

Kriging模型是一種常用的空間插值方法,可以根據已知數據點的空間分布估計未知位置的數值[25].Kriging模型的預測功能是基于協方差函數來建立空間數據之間的相似性關系,借助數據之間的相關性與空間距離的關系,根據距離越近的數據之間的相關性越高原則,擬合已知數據點的空間半變異函數進行預測.

2.1 基本概念

給定一組有 n 個的設計點 其中 xi∈IRn ,對應 X 的響應值 Y=[y1,y2,…,yn]T ,其中 yi∈IRq, 設計點 X 以及響應值都滿足規范化條件[27].

G(X)=fT(x)β+Z(x)

式中 :f(x)=[f1(x),f2(x),…,fn(x)]T, ,表示均值函數,反映模型的回歸趨勢 I(x) 的值在本文中實際為單位向量 I ,即一個長度為 Ωn 的單位列向量; β= [β1,β2,…,βn]1 是系數向量; Z(x) 是具有零均值的平穩高斯隨機過程,可以表示為 N(0,σ2) ,描述了模型與其潛在趨勢的偏離程度.在實踐中,平均函數最常見的形式是常數或者是以 x 為輸人的線性函數,文中采用的是常數均值.

協方差函數的作用是捕捉 G(X) 與平均值的偏離程度.除此之外,協方差可以調節這些偏差的大小 、G(X) 的粗糙度和縱向偏離的尺度.關于各種協方差函數,統計學中經常使用馬氏協方差,在涉及確定性計算機模型近似時,高斯協方差函數的使用可以大大地降低計算的成本,故高斯協方差占主導地位.

在高斯隨機過程中,兩個樣本點之間的協方差函數表示如式(8):

式中: R 是描述 xi 與 xj 之間相關性的相關函數; σ2 表示隨機過程的方差.用高斯方法表示如式(9):

式中: n 表示變量的維數; xi,l 和 xj,l 分別是 xi 和 xj 的第l個變量; θl 表示相關參數,該參數控制著樣本點 xi 和xj 之間第 l 個變量的相關性.

基于現有的樣本點,可以計算得到回歸系數 β 以及隨機過程方差 σ2 的估計值:

式中:A是 k×1 的向量,里面所有的元素為1;Y為列向量; R 為相關矩陣,通過相關參數求得.相關參數的確定方式有交叉驗證法和最大似然估計法,在大部分文獻中借助的是最大似然估計法.通過最大似然估計計算如下:

式中: p 為訓練樣本的總數.通過足夠的訓練樣本得到滿足精度要求的Kriging模型之后,可以基于現有的模型對任意未知點進行預測.其中預測過程不僅可以得到未知點的值,而且能夠得到預測值的波動范圍.預測的最佳線性無偏估計和方差如下:

(14)式中: r0 為相關系數向量, r0= [R(x0,x1),R(x0,x2),…,R(x0,xk)]T ,預測的均值為μ=ATR-1r0-A. 預測點的方差可以用于評估預測值的波動范圍.方差越小代表波動越小,得到的預測誤差越小.

2.2學習準則

在大多數情況下,僅通過初始訓練樣本構建的初始Kriging模型,其精度通常較低.若直接使用初始Kriging模型進行可靠性分析,將不可避免地產生顯著誤差,進而引發工程問題的失敗.面對這類問題,研究人員通常會尋找一些適當的點來加人訓練集,以便更新和修正Kriging模型.在輸入樣本集中,尋找合適的樣本點加入訓練集中的準則,稱為加點準則.因為這是主動尋找樣本點的過程,所以稱之為主動學習.目前常見的加點準則有:高效的全局優化(efficientglobaloptimization,EGO)準則,使用期望改進作為搜索點的選擇標準,把期望改進最大的點作為新的訓練樣本;在使用Kriging模型計算函數的極值時,通過期望改進函數(expected improvement func-tion,EIF)準則構造學習函數,在能最大程度改進當前極值的地方增加訓練樣本;結合Kriging模型和蒙特卡羅的主動學習可靠性分析方法( U 準則),使用 U 學習函數在功能函數的符號誤判概率最大的地方增加訓練樣本[31].

3算法構造

本文提出的算法主要圍繞基于證據理論的主動學習Kriging模型展開.首先,利用拉丁超立方抽樣生成初始訓練樣本,確保樣本均勻分布在整個輸人變量空間中.接著,構建初始Kriging模型,并結合優化方法與主動學習過程,通過 U 學習函數識別需要修正的區域,利用內點法優化算法在整個輸入變量空間中搜索最佳訓練樣本,以修正和優化Kriging模型,然后,采用區間蒙特卡羅抽樣方法減輕高維問題的計算負擔,并對焦元進行分類,計算結構的可信度和似真度.整個算法流程包括證據理論處理認知不確定性、Kriging模型構建與更新以及可信度和似真度的計算,通過提高Kriging模型的構建效率和成功率,降低傳統方法中搜索訓練樣本的約束,從而有效解決可靠性分析問題.

3.1主動學習Kriging模型

主動學習Kriging模型的構建主要包括:確定初始訓練集、利用初始訓練集構建初始的Kriging模型、結合優化方法和主動學習過程搜索最佳訓練樣本、依據停止條件判斷Kriging模型的精度.構建主動學習Kriing模型的具體流程如下:

由于在“黑箱\"模型中,在進行采樣時無法判斷在Kriging模型中什么區域需要更多的樣本點,什么區域需要更少的樣本點.所以初始訓練樣本的采樣原則期望盡可能使其均勻地分布在整個輸入變量空間.故采用拉丁超立方抽樣,該抽樣方法常用于設計實驗和可靠性分析.優點是:能夠提供較好的抽樣均勻性和覆蓋性,有效捕捉多個變量之間的關系;在相同的樣本數量下,LHS可以提供更多的信息和數據密度.在利用LHS時,首先在 N 維的超立方空間[0,1]?N 中抽取 n 個樣本;然后將這 n 個樣本映射到輸入變量空間中抽取樣本,將其記為 Xin ,根據經驗,抽取的樣本數量 n=(N+1)(N+2)/2 ;最后利用 Xin 調用功能函數求解相應的響應值 ,由 Xin 組合得到訓練集

構建Kriging模型時采用高斯型相關函數(也稱 高斯核函數).

根據2.1節中Kriging模型的建立過程,通過訓練集 (Xin,gin) ,可以計算得到回歸系數 β 以及隨機過程方差 σ2 的估計值,從而得到Kriging模型 通過 可以完成對輸人變量空間內的任意點 x 的預測,獲得點 x 的均值 ,方差 ,其中

通過Kriging模型不僅可以獲得該點的預測值,還可以得到預測值的偏離程度等重要信息.可靠性分析是一個分類問題,可以通過未知焦元的符號判斷其類型,進而實現可靠度計算.Kriging模型的功能是為未知點處提供符號預測,并不需要輸出極限狀態函數的真實值,從而構建一個局部逼近極限狀態面的局部Kriging模型.

針對得到的初始Kriging模型,需要對其進行精度判斷,在Kriging模型精度不夠的區域加入訓練樣本,修正Kriging模型,使模型逼近真實的功能函數,達到要求的精度.難點在于如何搜索到最佳的訓練樣本.

根據訓練樣本的特點,可以分為3類:第1類為極限狀態面附近的點,這一類訓練樣本的特點是預測均值 在極限狀態面附近,在預測工程中,預測值有較大概率穿越極限狀態面,導致符號預測錯誤;第2類訓練樣本的不確定性程度高,即 過大,這類訓練樣本預測值 的分布分散,其預測值同樣有較大的可能性穿越極限狀態面;第3類訓練樣本同時具備前兩類訓練樣本的共同特征,僅僅考慮均值和方差無法將其準確地找出來.

搜索最佳訓練樣本策略如下:

通過拉丁超立方抽樣(LHS)生成初始訓練樣本,確保樣本均勻分布;利用Kriging模型中 U 學習函數判斷預測輸入變量的均值和方差的正確概率,識別需要修正的區域;然后通過停止條件控制Kriging模型循環更新的精度;最后結合優化算法內點法,從整個輸人變量空間中搜索使預測均值或方差最小的點,作為新的訓練樣本加入訓練集中,以修正和優化Kriging模型.

3.1.1 U 學習函數

U 學習函數最初是由Echard等2提出來的,可以用于間接度量Kriging模型預測值符號錯誤概率的大小,其公式如下:

式中: 代表點 x 的預測均值; 代表預測標準差.

可靠性分析是一個二分類的過程,可以通過焦元上極大值點和極小值點的正負判斷焦元位于安全域還是失效域內.所以Kriging模型符號預測的正確率對可靠性計算的準確性有很大影響.雖然借助 的大小,可以簡單地對符號預測錯誤的概率進行判斷,即 的值越大,符號預測錯誤的概率越大; 越小,符號預測錯誤的概率越小,但是這種方式無法完成第1類和第3類訓練樣本的判斷.為了保證 符號預測的正確率,依據 ,可以推導出 符號預測正確率概率的計算公式[31],如式(16):

式中: 為標準高斯變量的累積分布函數; sign(?) 是符號函數; 表示 g(x) 的預測值; 表示 g(x) 的符號預測; 表示Kriging模型預測的符號與真實功能函數符號相同的概率.式(16)成立的根本條件是任意點 x 處 ,且 不等于0.

U 學習函數用于判斷符號預測正確概率,如圖3所示,當預測均值 時,通過改變方差來獲得不同的概率密度函數. U(x) 反映在 x 點處預測值 波動范圍的大小, U(x) 越大, 取值的波動范圍就越小.當預測均值 為0.5、1和2時, U(x) 分別為4、2和1,對應的正態分布曲線分別為實線、虛線和點畫線.從分布曲線中可以明顯地看出實線的分布最集中,預測值 的波動范圍最小,對應的 U(x)=4 ,此時 ,當預測值 0時表示預測值與均值符號相同,根據正態分布的累積分布函數計算原理,圖3中實線與橫軸大于0的區域面積代表預測值 的概率,表明Kriging模型 在 x 點處符號預測正確的概率,此時預測正確的概率值為 p??(4) ,即 p≈1. 點畫線的分布最為分散,對應的預測值 的波動范圍最大,其U(x)=1 ,同樣的原理,點畫線對應的符號預測正確的概率為 p??(1) ,即 p?0.841B ,表明該條件下,符號預測錯誤的概率約為 16% 對于虛線, 的分布相對點畫線比較集中,它對應的 U(x)=2 ,其符號預測正確的概率為 p??(2) ,即 p?0.9772. 在相關文獻中,當 ?(U(x))=?(2) 時,可以認為有足夠高的概率Kriging模型能夠提供正確的符號預測.當符號預測為負號時,也可以做出以上同樣的推斷,

圖3點 x 處預測均值為2時,不同方差下正態分布圖Fig.3Normal distribution plots at point x with predictedmean2underdifferentvariances

3.1.2停止條件

在Kriging模型更新的循環過程中,停止條件用來控制模型的精度.借助 U 學習函數衡量模型的精度,當Kriging模型預測的 U(x) 的最小值滿足 U(x)?2 時,可以認為Kriging模型已經達到精度要求;當Kriging模型預測的 U(x) 的最小值滿足 U(x)lt;2 時,則繼續加入新的訓練樣本修正Kriging模型.

3.1.3內點法優化

內點法(interiorpointmethod)是一類求解凸優化問題的數值方法,通過在可行域內部進行迭代搜索進而逼近最優解[35-37].基本思想是將原始問題轉化為一個帶約束的優化問題,通過迭代的方式沿著可行域的內部路徑逐步逼近最優解.內點法的應用在于從整個輸入變量空間中搜索使 U(x) 最小的點.

min U(x)

式中: 表示Kriging模型; x 表示輸人變量空間 X 中的任意點;predictors是Kriging模型的預測函數; 表示在輸入變量為 x 時 的預測均值; 表示在輸入變量為 x 時 的預測誤差.

本文結合優化方法和主動學習過程,從整個輸入變量空間中搜索最佳訓練樣本,從而更好地提高了Kriging模型的構建效率以及成功率.

首先,根據不確定性確定證據輸人變量空間 X

式中: XiL 表示第 i 個變量的下限值; XiU 表示第 i 個變量的上限值.由求解空間 X 的中心點 Xc

然后將中心點 Xc 作為優化的初始點 X0 ,設置評估次數以及容差,調用MATLAB中的fmincon工具箱進行最優化求解,根據式(17)輸出 U(x) 的最小值.如果 U(x) 最小值 lt;2 ,該點即為下一個訓練樣本.

式中: xnext 為下一個訓練樣本,即最佳訓練樣本.根據以上主動學習過程,定位到最佳訓練樣本的位置.然后根據 xnext 的值進行有限元仿真計算求解得到 ,并將 (xnext,gnext(xnext)) 作為新的訓練樣本加入訓練集 (Xin,gin) 中,修正Kriging模型,直到滿足精度要求.

3.2區間蒙特卡羅抽樣

盡管使用Kriging模型能夠大幅地降低計算的成本和難度,但是在面對高維問題的時候,仍然會面臨沉重的計算負擔.對于復雜問題,構建的Kriging模型同樣會比較復雜,在未知點的預測解速度也會下降,要完成大規模焦元的預測,可能花費幾天的時間.實現大規模焦元極值的計算對計算機的性能要求也會大幅提升.因此,在計算可信度和似真度的過程中需要借助區間蒙特卡羅模擬(MonteCarlosimulation,

MCS)[30]來減輕計算的壓力.MCS的關鍵在于如何生成與基本可信度分配(BPA)相對應的區間樣本.對具有 N 個證據變量的輸入變量,用區間蒙特卡羅抽取 Nmc 個樣本,具體步驟如下:

步驟1:初始化 i=1,j=1

步驟2:在區間[0,1]上生成一個均勻分布的隨機數 u

步驟3:如圖4所示,對于第 j 個變量的第 i 次模擬,生成區間樣本 ,其中 如果滿足 其中 l=1,2,…,nj, 且 ,便可以將焦元 記為 Aji. 當抽取的樣本數量足夠大,在 Aji 內生成的樣本的概率等于m(Aji

步驟4:對步驟3進行循環計算,直到滿足 j= N,i=Nmc. (20號

圖4第j個變量的區間蒙特卡羅抽樣示意圖Fig.4TheintervalMonteCarlosamplingdiagramforvariable j

3.3可信度和似真度計算

在計算結構的可信度和似真度之前需要對焦元進行分類.其中Kriging模型為區間MCS抽取的焦元進行預測,考慮到功能函數的非線性,采用SQP優化算法對焦元的極值進行求解,由式(6)每一個焦元的極值可以表示為 (20其中 表示Kriging代理模型, Ai 表示焦元.

根據焦元極值進行分類,可以分為可靠焦元 G= ;邊界焦元 B=Ai{Ai [gminlt;0,gmaxgt;0]} 和失效焦元 0,gmaxlt;0] .完成焦元分類之后,根據式(5),可信度表示為 ;似真度表示為 Pl(G)= ,其中 Nmc 為使用區間蒙特卡洛從輸人變量空間中抽取的樣本個數.

3.4算法流程總結

算法的基本流程如圖5所示,主要分為3個模塊:第1個模塊在證據理論框架下對認知不確定性進行處理.第2個模塊Kriging模型的構建過程,其中Kriging模型的更新是難點,本文結合優化方法和主動學習過程,實現最佳訓練樣本的搜索;第3個模塊計算可信度和似真度,主要是進行IMCS抽樣和焦元類別的判斷,以及可信度和似真度的計算,

圖5基本算法流程Fig.5Thebasicalgorithmflowchart

4算例驗證

4.1 數值算例一

曲柄-滑塊機構[38通常用于工程機械中,其結構如圖6所示.在本次可靠性分析中,考慮曲柄滑塊機構的材料強度與最大應力之間的差值.結構的基本參數如下:桿 O1O2 的長度 Ψa ,桿 O2O3 的長度 b ,桿O2O3 為空心管,內徑為 d1 ,外徑為 d2,Oi 和 O3 之間垂直距離(偏心距)為 e,μ 為摩擦系數.在工作狀態時滑塊受到一個水平力,載荷為 P. 結構的許用應力 S= 2.1GPa ,極限狀態函數如下:

圖6曲柄-滑塊機構示意圖Fig.6Schematicdiagramofcrank-slidermechanism

曲柄-滑塊機構桿 O1O2 和桿 O2O3 長度 Φa 和 b 的BPA如表1所示,承受的載荷 P 及其偏心距 ρe 所對應的BPA如表2所示.

表1桿長a和b的BPA Tab.1 Basicprobabilityassignmentof rodsaandb
表2載荷 P 和偏心距 e 的BPA

基于變量信息,構建證據輸人變量空間 X. 每一變量的范圍為: X1=[94,106],X2=[290,310],X3= (204號 [220,280],X4=[100,150] ;通過LHS從 X 中抽取15個初始訓練樣本,代人功能函數 g(a,b,P,e) 中求解結構的響應值 ,得到初始訓練集 (Xin,gin) ;然后通過 (Xin,gin) 構建初始Kriging模型 ;接下來利用 提供預測,通過結合內點法和主動學習過程從X 中搜索最佳訓練樣本,判斷 U(x) 的最小值是否滿足停正條件 U(x)?2 ,若滿足停正條件,即可輸出Kriging模型,否則繼續更新Kriging模型直到所有U(x) 的最小值滿足 U(x)?2

為了驗證本文方法中新增訓練樣本的修正效果,將其與樣本點法進行比較.在Kriging模型更新過程中,樣本點法搜索到的第一個訓練樣本對應的U 函數值為0.001169,使用本文方法搜索到的第一個訓練樣本 U 函數的值為 4.94×10-5 ,在同樣的初始樣本下構建得到的Kriging模型,采用本文方法搜索到的第一個訓練樣本處,Kriging模型的預測錯誤率更大,在該位置處加入訓練樣本對Kriging模型的修正效果更好,表明本文方法搜索到的訓練樣本比樣本點法搜索到的訓練樣本修改效果更好.表3為樣本點法和本文方法構建Kriging模型使用的初始訓練樣本和新增訓練樣本數量.樣本點法和本文方法使用相同的初始訓練樣本條件下,樣本點法新增訓練樣本為14,一共用了29個訓練樣本,本文方法新增訓練樣本數為6,一共用了21個訓練樣本.相比之下,本文方法使用訓練樣本數比樣本點法有所減少,表明本文方法Kriging模型構建的效率得到了一定的提升.

表3構建Kriging模型的初始訓練樣本數和新增訓練樣本數Tab.3Initial trainingsample sizeand additional training sample size for constructing Kriging models個

為了驗證本文方法的計算效率,將其與優化法[9]和樣本點法進行對比.優化法[39]是一種傳統的證據可靠性分析方法,通常使用SQP優化算法對每一個焦元求解極值,進而對焦元進行分類,實現可信度和似真度的計算,該方法計算精度較高,通常作為其他方法的參考.使用SQP進行焦元極值計算時通常需要調用功能函數幾十次到幾百次,文中用 Φt 表示優化法調用功能函數的次數.在計算可信度和似真度時,采用了IMCS進行抽樣,樣本規模為 1×105 ,優化法調用功能函數次數為 Ωt×105 .樣本點法和本文方法使用Kriging模型代替功能函數,因此功能函數的調用次數為構建Kriging模型使用的樣本數.可以得到樣本點法和本文方法的功能函數調用次數分別為29和21.通過對比,本文方法對功能函數的調用次數最少,表明本文方法有較高的計算效率.

表4為優化法、樣本點法和本文方法的可靠性分析結果及其相對誤差.優化法計算的 Bel(G) 為0.8918,Pl(G) 為0.99783,以其為參考,樣本點法計算的 Bel(G) 為0.8923,相對誤差約為 0.056% Pl(G) 為0.9981,相對誤差約為 0.03% .本文方法計算的Bel(G) 為0.89185,相對誤差約為 0.006% Pl(G) 為0.99783,相對誤差為0.可見本文方法有較高的計算精度.通過這個算例可以看出本文方法可以用于解決可靠性分析問題.

表4可靠性分析結果及其相對誤差

4.2數值算例二

以13個獨立證據變量的組合梁結構的可靠性分析問題為例[40],驗證本文方法在高維高非線性問題中的適用性.組合梁的結構如圖7所示.

梁的楊氏模量為 Ew ,截面的寬度為A、高度為 B 、長度為 L. 梁下面安裝的鋁板,其楊氏模量為 Ea ,截面的寬度為 C 、高度為 D ,其長度和梁一樣,此外該鋁板固定在梁的底面.在梁的上方施加有外力 P1,P2 、P3,P4,P5 和 P6 ,分別對應梁的 L1,L2,L3,L4,L5,L6 位置,方向垂直于梁平面.許用的拉應力為S.梁承受的最大應力出現在中間的截面 M-M 中,為了保證組合梁的安全性,梁所承受的最大應力 σmax 必須小于梁的最大許用應力 s ,組合梁安全使用的功能函數如下:

G(X,Y)=S-σmax

其中組合梁承受的最大應力公式如下:

組合梁結構承受的載荷變量的BPA結構如表5和表6所示,梁的材料強度的BPA結構如表7所示,梁的截面尺寸變量的BPA結構如表8所示,以及鋁板截面尺寸的BPA結構如表9所示.

表5載荷變量 P1,P2,P3 的BPATab.5Basic probabilityassignment of P1,P2,P3

表10樣本點法和本文方法構建Kriging模型時的初始訓練樣本數和新增訓練樣本數.使用相同的初始訓練點得到初始Kriging模型,樣本點法沒有完成Kriging模型的更新,這是由于在更新過程中,樣本點法重復搜索到同一個訓練樣本,導致Kriging模型的更新過程無法繼續進行.本文方法的初始訓練樣本數為105,新增訓練樣本數為83,一共使用188個訓練樣本,對應的調用功能函數188次.表明在高維問題中,本文方法能夠成功構建Kriging模型并且有較高的計算效率.

表6載荷變量 P4,P5,P6 的BPATab.6Basic probability assignment of P4,P5,P6
表8梁截面尺寸的BPA
表9鋁板截面尺寸的BPA

表11為優化法、樣本點法和本文方法的計算結果.樣本點法因為沒有成功構建Kriging模型,所以沒有得到可靠性分析結果.優化法計算的 Bel(G) 為0.81283,Pl(G) 為0.99611.本文方法計算的 Bel(G) 為0.81720,相對誤差約為 0.538% ,Pl(G)為0.99611,相對誤差為0,可見本文方法有較高的計算精度.該算例表明本文方法可以用于高維可靠性分析問題的求解.

表10構建Kriging模型的初始訓練樣本和新增訓練樣本數 Tab.10 Initial training sample size and additional training sample size for constructingKriging models個
表11可靠性分析結果及其相對誤差Tab.11 Reliabilityanalysis resultsand their relative errors

4.3車輛正面碰撞可靠性分析

將本文方法應用到車輛正面碰撞的可靠性分析問題中,測試其在實際工程問題中的計算效果.在汽車事故中,很大的比例是由汽車正面碰撞[4造成的,所以對汽車正面碰撞的安全性分析對保護人的生命健康安全有很大的意義.據分析,車輛正面碰撞時危害主要來自車身的加速度.圖8為車輛正面碰撞示意圖[41],通過車輛的防撞梁、前軌、前軌罩和加強板結構實現吸收碰撞能量的功能,其中防撞梁的厚度為 t1 ,前軌的厚度為 t2 ,前軌罩的厚度為 t3 ,加強板的厚度為 t4 及材料強度為 E 和密度為 ρ ,它們均被看作是證據變量而且兩兩之間是相互獨立的,防撞梁的結構尺寸、材料強度和密度的焦元及BPA如表12和表13所示.考慮加速度的影響,建立了功能函數如下:

ga=amax-a(t1,t2,t3,t4,E,ρ)

為了提高計算的效率,使用二次響應面法和拉丁超立方設計構建了加速度的響應面a(t1,t2,t3,t4,E,ρ)

圖8車輛正面碰撞示意圖Fig.8 Diagram of frontal vehiclecollision
表13防撞梁尺寸及其材料強度和密度的BPA
表12防撞梁尺寸的BPA Tab.12 Basic probabilityassignmentof bumperbeamdimensions

表14為樣本點法和本文方法構建Kriging模型時的初始訓練樣本數和新增訓練樣本數.在相同的初始訓練樣本下,樣本點法的新增訓練樣本數為153,一共使用181個訓練樣本,即調用功能函數181次;本文方法新增訓練樣本數為30,一共使用58個訓練樣本,即調用功能函數58次,可見本文方法相比于樣本點法計算效率有一定的提升.

表15為使用優化法、樣本點法和本文方法的計算結果.優化法計算得到的 Bel(G) 為 0.92229,Pl(G) 為0.99678.樣本點法計算得到的 Bel(G) 為0.96111,相對誤差約為 4.209% Pl(G) 為0.96111,相對誤差約為 3.579% ,相對誤差較大,計算結果失效.本文方法計算得到的 Bel(G) 為0.96127,相對誤差約為0.067% Pl(G) 為0.99681,相對誤差約為 0.003% ,可見本文方法有較高的計算精度.本文方法對車輛正面碰撞的可靠性分析問題的成功求解,表明本文方法具有一定的工程價值.

表14構建Kriging模型的初始訓練樣本數和新增訓練樣本數
表15可靠性計算結果及相對誤差Tab.15Reliabilityanalysisresultsandtheirrelativeerrors

5結論

本文提出了一種基于證據理論的主動學習可靠性分析方法,該方法將優化方法和主動學習過程相結合,從整個輸入變量空間中搜索訓練樣本,確保每一次搜索得到的訓練樣本為最佳訓練樣本,提高了Kriging模型的構建效率和成功率.并通過數值算例驗證了該方法在解決可靠性分析問題的有效性,最后將該方法應用到車輛正面碰撞的可靠性分析問題中,成功完成了可靠性分析計算.可見,本文方法在解決可靠性分析問題時有較高計算效率和計算精度,在未來,希望本文方法能考慮變量之間的相關性,解決更多的實際工程問題,此外將其拓展到混合不確定性的可靠性分析問題中.

參考文獻

[1] FERSON S,GINZBURG L R.Different methods are needed to propagate ignorance and variability[J].Reliability Engineeringamp; System Safety,1996,54(2/3):133-144.

[2] HELTONJC.Uncertainty and sensitivity analysis in the presence ofstochastic and subjective uncertainty[J]. Journal of Statistical ComputationandSimulation,1997,57(1/2/3/4):3-76.

[3] DEMPSTERAP,LAIRDNM,RUBINDB. Maximumlikelihood fromincomplete data via the EM algorithm[J].Journal of the RoyalStatisticalSocietySeriesB:StatisticalMethodology,1977, 39(1):1-22.

[4] JIANGC,LIWX,HANX,etal.Structuralreliabilityanalysis based on random distributions with interval parameters[J]. Computersamp; Structures,2011,89(23/24):2292-2302.

[5] MOLLERB,GRAFW,BEER M. Fuzzy structural analysisusing α-level optimization[J].Computational Mechanics,2Ooo,26(6): 547-565.

[6] ZADEHL A.Fuzzy sets asa basis fora theory of possibility[J] FuzzySetsandSystems,1999,100:9-34.

[7] DEMPSTER AP.Upper and lower probabilities induced by a multivalued mapping[J].The Annals of Mathematical Statistics, 1967,38(2):325-339.

[8]張哲.基于證據理論的結構可靠性分析方法[D].長沙:湖南 大學,2016. ZHANG Z. Structural reliabilityanalysis method based on evidence theory[D].Changsha:Hunan University,2016.(in Chinese)

[9]YIN SW,YUDJ,YINH,et al.A new evidence-theory-based method for response analysis of acoustic system with epistemic uncertaintybyusing Jacobi expansion[J].ComputerMethods in Applied Mechanics and Engineering,2017,322: 419-440.

[10]劉鑫,龔敏,周振華,等.基于證據理論的機械結構高效可靠性 分析方法[J].中國機械工程,2020,31(17):2031-2037. LIU X,GONG M,ZHOU Z H,et al. An efficient mechanical structure reliability analysis method based on evidence theory[J]. China Mechanical Engineering,2020,31(17):2031-2037.(in Chinese)

[11] DONG W M,SHAH H C.Vertex method for computing functions of fuzzy variables[J].Fuzzy Setsand Systems,1987,24(1): 65-78.

[12]GAO RG,YIN SW,XIONG F. Response analysis and reliability-based design optimization of structural-acoustic system underevidence theory [J].Structural and Multidisciplinary Optimization,2019,59(3):959-975.

[13]BAEHR,GRANDHI RV,CANFIELD RA.An approximation approach foruncertainty quantification using evidence theory[J]. ReliabilityEngineeringamp;SystemSafety,2004,86(3):215-225.

[14]ZHANG Z,JIANGC,HAN X,et al. A response surface approach forstructural reliability analysis using evidence theory [J]. AdvancesinEngineeringSoftware,2014,69:37-45.

[15] ZHANG Z,RUAN X X,DUAN M F,et al. An efficient epistemic uncertainty analysis method using evidence theory[J].Computer Methods in Applied Mechanics and Engineering,2O18,339: 443-466.

[16] MOURELATOS Z P,ZHOU J. A design optimization method using evidence theory[J].Journal of Mechanical Design,2006, 128(4):901-908.

[17]姜潮,張哲,韓旭,等.一種基于證據理論的結構可靠性分析 方法[J].力學學報,2013,45(1):103-115. JIANG C, ZHANG Z,HAN X, et al. An evidence-theory-based reliability analysis method for uncertain structures [J]. Chinese Journal of Theoretical and Applied Mechanics,2013,45(1): 103-115.(in Chinese)

[18]XIAO M, XIONG H H,GAO L,et al.An efficient method for structural reliability analysisusingevidence theory[C]//2014 IEEE17th International Conference on Computational Science and Engineering.December19-21,2014,Chengdu,China. IEEE,2014:144-149.

[19] ZHANG Z, JIANG C,WANG G G,et al.First and second order approximate reliability analysis methods using evidence theory[J]. Reliability Engineeringamp; System Safety,2015,137:40-49.

[20]BAE H R,GRANDHI R V,CANFIELD R A.Epistemic uncertainty quantification techniques including evidence theory forlarge-scale structures[J].Computersamp; Structures,2004, 82(13/14): 1101-1112.

[21]ZHANG Z,JIANG C,WANG G G,et al. An efficient reliability analysis method for structures with epistemic uncertainty using evidence theory[C]//International Design Engineering Technical Conferences and Computers and Information in Engineering Conference.American Society of Mechanical Engineers,2014, 46322:V02BT03A055.

[22]曹亮,龔曙光,陳國強,等.基于支持向量機和證據理論的復雜 系統可靠性分析方法[J].機械設計,2024,41(5):131-137. CAOL,GONGSG,CHENGQ,etal.Method ofreliability analysis on complex systems based on support vector machine and evidence theory[J].Journal of Machine Design,2024,41(5): 131-137.(in Chinese)

[23]YIN S W,YU D J,LUO Z,et al. An arbitrary polynomial chaos expansion approach for response analysis of acoustic systems with epistemic uncertainty [J].ComputerMethods in Applied Mechanics and Engineering,2018,332:280-302.

[24]JONES D R,SCHONLAU M,WELCH W J.Efficient global optimization of expensive black-box functions[J]. Journal of Global Optimization,1998,13(4):455-492.

[25]KAYMAZ I. Application of Kriging method to structural reliability problems[J].Structural Safety,2005,27(2):133-151.

[26]ECHARD B,GAYTON N,LEMAIRE M.AK-MCS:an active learning reliability method combining Kriging and Monte Carlo simulation[J].Structural Safety,2011,33(2):145-154.

[27]LV Z Y,LU Z Z,WANG P.A new learning function for Kriging and its applications to solve reliability problems in engineering [J]. Computers amp; Mathematics with Applications,2015,70(5): 1182-1197.

[28]YANG X F,LIU Y S,GAO Y. Unified reliability analysis by active learning Kriging model combining with random-set based Monte Carlosimulationmethod[J].International Journal forNumerical Methods in Engineering,2016,108(11):1343-1361.

[29] ZHANG D Q,LIANG Y F,CAO L X,et al. Evidence-theorybasedreliability analysis through Kriging surrogate model[J]. Journal of Mechanical Design,2022,144(3):031701.

[30]YANG X F,LIU Z Q,CHENG X. An enhanced active learning Krigingmodel for evidence theory-based reliabilityanalysis[J]. Structural and Multidisciplinary Optimization,2O21,64(4): 2165-2181.

[31]陳澤權.基于Kriging模型的結構可靠性及全局靈敏度自適應 分析方法研空「n]長卷,吉林士學 2023 UHEN Z Q. Kesearcn on Adaptve analysis metnod Ior structural reliabilityand global sensitivity based on kriging model[D]. Changchun:Jilin University,2O23.(in Chinese)

[32]韋新鵬,姚中洋,寶文禮,等.一種基于主動學習克里金模型 的證據理論可靠性分析方法[J].機械工程學報,2024,60(2): 356-368. WEIXP,YAOZY,BAOWL,etal.Evidence-theory-based reliabilityanalysismethodusingactive-learningKriging model[J].Journal of Mechanical Engineering,2024,60(2): 356-368. (in Chinese)

[33]ZHANG Z,JIANG C.Evidence-theory-based structural reliabilityanalysiswith epistemic uncertainty:a review[J]. Structural and Multidisciplinary Optimization,2021,63(6):2935- 2953.

[34]于俊濤,鄧衛,王巨,等.基于近似移動矢量的證據理論可靠 性設計優化方法[J].湖南大學學報(自然科學版),48(8): 59-67 YUJT,DENGW,WANGJ,etal.Anevidence-theory-based design optimization method using approximate shifting vector.[J] Journal of Hunan University(Natural Sciences),48(8):59-67. (in Chinese)

[35]高峰,張連生.線性約束凸規劃內點法及其修正算法[J].運 籌學學報,1998,2(1):79-94. GAOF,ZHANG L S. Interior point methods for convex programming with linear constraint and their modified algorithms [J].Operation Research Transactions,1998,2(1):79-94.(in Chinese)

[36]EL-BAKRY A S,TAPIARA,TSUCHIYAT,etal.On the formulationand theory of the Newton interior-point method for nonlinear programming[J]. Journal of Optimization Theory and Applications,1996,89(3): 507-541.

[37]YAMASHITA H. A globally convergent primal-dual interior point method for constrained optimization[J].Optimization Methods and Software,1998,10(2):443-469.

[38]GUO J,DU X P.Sensitivity analysis with mixture of epistemic and aleatoryuncertainties[J].AIAAJournal,2007,45(9):2337- 2349.

[39]ROBINSONJA.Probability and statistics for engineers and scientists[J].Technometrics,1990,32(3):348-349.

[40]HUANG B Q,DU X P.Probabilistic uncertainty analysis by mean-value first order Saddlepoint Approximation[J].Reliability Engineeringamp;System Safety,2008,93(2):325-336.

[41]JIANG C,ZHANG W,WANG B,et al.Structural reliability analysis using a copula-function-based evidence theory model [J].Computers amp; Structures,2014.143: 19-31

猜你喜歡
方法模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
學習方法
3D打印中的模型分割與打包
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 亚洲av无码专区久久蜜芽| 中文字幕资源站| 91精品国产福利| 日韩AV手机在线观看蜜芽| 四虎精品黑人视频| 亚洲无码不卡网| 黄色网在线免费观看| 国产特级毛片| 亚洲精品成人片在线观看| 中文字幕有乳无码| 高h视频在线| 日本欧美成人免费| 亚洲一区第一页| 成人久久精品一区二区三区| 国产成人凹凸视频在线| 国产91麻豆免费观看| 欧美 亚洲 日韩 国产| 久久亚洲综合伊人| 亚洲v日韩v欧美在线观看| 91久久大香线蕉| 日韩a级片视频| 青青青草国产| 欧美一级片在线| 又大又硬又爽免费视频| 在线观看网站国产| 欧洲精品视频在线观看| 国产乱码精品一区二区三区中文 | 午夜在线不卡| 国产精品久久久久久久久kt| 六月婷婷精品视频在线观看 | 日韩在线第三页| 久久99热这里只有精品免费看| 久久青青草原亚洲av无码| 国产精品99r8在线观看| 伊人天堂网| 5388国产亚洲欧美在线观看| 久久人妻xunleige无码| 刘亦菲一区二区在线观看| 呦女亚洲一区精品| 日韩不卡高清视频| 亚洲中文字幕23页在线| 夜夜爽免费视频| 国产91透明丝袜美腿在线| 免费国产小视频在线观看| 国产美女在线免费观看| 中文字幕天无码久久精品视频免费| 成人在线观看一区| 国产18页| 伊人五月丁香综合AⅤ| 久久夜夜视频| 国产不卡网| 狠狠久久综合伊人不卡| 欧美日韩在线成人| 国产91熟女高潮一区二区| 青青青伊人色综合久久| 久久大香香蕉国产免费网站| 国产jizzjizz视频| 国产凹凸一区在线观看视频| 亚洲午夜福利在线| 91久久精品国产| 国产在线视频导航| 8090成人午夜精品| 欧美午夜理伦三级在线观看| 国产爽歪歪免费视频在线观看| 最新国产高清在线| 国产精品美乳| 欧美国产日韩在线| 日韩无码视频专区| 人妻中文久热无码丝袜| 国产日韩欧美成人| 91精品免费久久久| 国产精品午夜福利麻豆| 免费看一级毛片波多结衣| 在线不卡免费视频| 91精品综合| 欧美日本视频在线观看| 一级毛片免费不卡在线视频| 福利姬国产精品一区在线| 欧美日本视频在线观看| 香蕉视频在线观看www| 夜夜爽免费视频| 欧美日本在线一区二区三区|