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

基于數據驅動的結構鋼表面應力磁巴克豪森噪聲表征方法

2023-06-27 04:58:16崔西明邱志鵬魏嘉張弛宋凱李喆王樹鵬
航空學報 2023年8期
關鍵詞:信號模型

崔西明,邱志鵬,魏嘉,張弛,宋凱,*,李喆,王樹鵬

1.南昌航空大學 無損檢測技術教育部重點實驗室,南昌 330063

2.中國航發(fā)沈陽黎明航空發(fā)動機有限責任公司,沈陽 110043

飛機起落架、機翼主梁、平尾大軸、直升機旋翼軸、接頭和對接螺栓等鐵磁結構鋼航空構件在服役過程中經常承受局部應力集中,容易產生因超過疲勞強度而導致的裂紋,威脅航空構件的安全運行[1]。通過定量評估300M 超高強度鋼等鐵磁構件表面各部位殘余應力狀態(tài)對早期疲勞損傷進行預測,對預防結構鋼失效引起的重大飛行安全事故具有重要意義[2-3]。

磁巴克豪森噪聲(Magnetic Barkhausen Noise, MBN)檢測技術廣泛用于鐵磁材料微觀結構和殘余應力的定量評估,與現有的常見應力測量方法(如X 射線法、金屬磁記憶法、電阻應變片法)相比,具有便攜、可靠性高、響應速度快的優(yōu)點[4-5]。然而,磁巴克豪森噪聲除了受材料內部殘余應力影響[6-8],對材料表面硬度、內部晶粒尺寸、微觀組織變化也很敏感。因此,很難用確定的解析物理模型對磁巴克豪森噪聲隨應力的非線性變化規(guī)律進行準確描述[9-10]。

基于數據驅動的非線性映射算法,可在不確定測量系統(tǒng)精確解析模型的情況下完成輸入輸出數據的擬合。根據通用近似定理,反向傳播(Back Propagation,BP)神經網絡可實現以任意精度無限逼近任意非線性函數,具有根據磁巴克豪森噪聲輸入實現應力高精度預測的可行性[11-13]。BP 神經網絡是一種多層前饋網絡或誤差反向傳播網絡[14],由輸入層、隱藏層和輸出層3層結構組成,以網絡誤差平方為目標函數,采用梯度下降法計算目標函數的最小值,解決了多層網絡中隱含層節(jié)點連接權值調整的問題,并具有自學習能力、非線性映射能力、對任意函數的逼近能力、并行計算能力和容錯能力等[15]。但BP神經網絡用于鐵磁材料表面應力的定量預測時[16],存在易受噪聲影響、特征選取困難、模型復雜度受輸入節(jié)點個數影響等問題。

國內外學者針對相關問題在MBN 特征量選取、信號降噪、數據降維以及預測模型構建等方面開展了深入研究。張翔等[17]利用連續(xù)小波變換提取出了磁巴克豪森噪聲信號各小波分解量的能量,并結合全局閾值部分重構降噪法獲得了更好的降噪效果。王平等[18]提出采用小波分解獲得更多BP 神經網絡輸入量,加入小波分析后,神經網絡相對誤差相比于之前減小了約90%,符合應力檢測要求。但小波變換在分析信號時只能對低頻段范圍內的信號進一步分解,Coifman團隊[19]為獲得更加全面的信號分析結果,細化分析信號高頻段,克服了小波變換只能進一步分解低頻信號的固有缺陷,利用小波包變換對風力發(fā)電機組齒輪箱振動信號進行分解。吳杰等[20]通過改進磁化器的設計,結合BP 神經網絡修正了提離狀態(tài)下的MBN 信號均方根值,大幅減小了提離效應的影響。張傳棟等[21]采用BP 神經網絡對鋼軸表面硬度進行預測,并優(yōu)化了鋼軸表面MBN 檢測傳感器,實現了鋼軸表面硬度的無損檢測。雷瑛[22]和Pandit[23]等采用數理統(tǒng)計對鋼表面殘余應力的極差和方差進行分析,利用線性回歸法建立了磨削加工表面的殘余應力模型,分析了磨削階段的機械應力變化特征。董雷[24]搭建了一套電磁融合檢測系統(tǒng),并通過電磁融合檢測系統(tǒng)建立了BP 神經網絡模型,對比了單一特征值和多特征值的擬合效果,表明了多特征值在應力表征和缺陷深度檢測中都取得了較好的效果。國內外學者對應力表征框架的設計關注相對較少,因此構建一種可供參考的表征算法框架對提升應力預測精度具有重要意義。

提出了一種基于數據驅動框架建立MBN 噪聲和應力映射關系的預測模型,采用小波包變換提取磁巴克豪森噪聲的小波包變換系數作為表征鐵磁材料表面應力的特征向量;采用奇異值分解進行特征向量的數據降維,減小特征向量的冗余性;利用BP 神經網絡建立奇異值向量與應力的映射模型,以實現鐵磁材料表面應力的定量評估。

1 磁巴克豪森噪聲檢測系統(tǒng)

磁巴克豪森噪聲是鐵磁材料動態(tài)磁化時磁疇壁跳躍過程中產生的一種脈沖信號,如圖1 所示,圖中M為磁化強度,H為外加磁場強度,Mr為剩磁,Hc為矯頑力,Hs為磁化強度達到飽和時對應的外加磁場強度。鐵磁材料在磁化過程中,應力會影響微觀組織對磁疇壁的釘扎作用,進而影響磁疇壁的跳躍,對磁巴克豪森噪聲信號產生影響。因此磁巴克豪森噪聲檢測法成為評估鐵磁材料應力分布情況的一種重要方法,根據應力對鐵磁材料內部磁疇的影響,可實現應力的定量評估。

圖1 MBN 響應示意圖Fig.1 MBN response diagram

基于MBN 檢測原理,搭建了MBN 應力檢測系統(tǒng),如圖2 所示。系統(tǒng)由硬件和軟件兩部分組成:硬件部分包括激勵模塊、MBN 傳感器、信號調理模塊、AD 采集模塊等;軟件部分由基于Lab-VIEW 設計的具有數據采集和數據分析處理功能的上位機組成。激勵模塊由AFG3102C 任意波形函數信號發(fā)生器以及HSA4101 型高速雙極功率放大器組成,激勵信號采用頻率為45 Hz 的正弦波信號,信號發(fā)生器輸出的原始信號峰峰值為300 mV;AD620 電壓放大器和NF3628 可編程濾波器串聯(lián)組成信號調理電路,接收增益設置為125 倍;AD 采集模塊采用NI 公司的USB-6366 型高速同步采集卡,采樣頻率為2 MS/s。

圖2 MBN 檢測系統(tǒng)示意圖Fig.2 Schematic diagram of MBN detection system

由圖2 可知,MBN 傳感器由U 型磁軛、激勵線圈、檢測線圈3 個部分組成。其中,U 型磁軛材料采用的是高磁導率、低電導率的錳鋅鐵氧體,磁極端面為正方形平面,兩磁極中心點間距為19 mm。激勵線圈、檢測線圈均選用漆包線,分別纏繞在U 型磁軛、磁棒上,檢測線圈擺放于U 型磁軛下方正中心位置。當正弦信號通過激勵線圈時,U 型磁軛兩磁極之間會產生近似均勻的交變磁場,從而在U 型磁軛與被測試塊中形成磁回路,以產生毫伏級的微弱MBN 信號。檢測線圈將拾取到的MBN 信號通過電纜導線傳輸到信號調理模塊進行濾波、放大等處理,使其便于由AD采集模塊送入上位機實現信號顯示、數據存儲等功能。圖3 為濾波前后的MBN 信號波形對比。

圖3 MBN 信號濾波前后對比示意圖Fig.3 Comparison diagram of MBN signal before and after filtering

2 BP 神經網絡應力預測模型

2.1 磁巴克豪森噪聲信號特征向量提取

統(tǒng)計特征量主要反映了MBN 信號的時域特征,為融合MBN 信號的頻域信息,對MBN 信號進行時頻分析。小波包變換是一種具有較高時-頻分辨率的時頻分析方法,既可對信號低頻部分進行分解,還可對高頻部分進行分解?;谛〔ò儞Q的信號分解既無冗余,也無疏漏,對包含大量中、高頻信息的MBN 信號能進行時頻局部化分析,可用于MBN 信號時頻特征提取。

對于任意MBN 信號f(t)的小波包變換為

利用小波包變換得到MBN 信號各頻帶階次。小波包系數,圖4 為小波包變換中樹分解的示意圖,圖中LPF(Low Pass Filter)為低通濾波器,HPF(High Pass Filter)為高通濾波器,L1 為原始信號中的低頻信息,H1 為原始信號中的高頻信息。小波包樹的每個節(jié)點都有對應的小波包系數,決定了頻率的大小,攜帶了MBN 信號的局部時頻信息。根據計算得到的小波包分解系數,可構建小波包系數矩陣A:

圖4 小波包分解樹Fig.4 Wavelet packet decomposition tree

式中:n=2j為j層小波包分解得到的頻帶總數;m為每個頻帶的小波包系數長度。

利用3 層小波包變換對MBN 信號進行分解,導出不同分解尺度上的MBN 信號小波包變換系數,并按尺度順序排列成特征向量供識別使用。

系數矩陣A反映了MBN 信號在時頻域上的特征,但由于矩陣維度較大,需進一步采用奇異值分解,對系數矩陣A進行降維。設A∈Rm×n,Rm×n為實矩陣,則存在正交矩陣U∈Rm×n和正交矩陣V∈Rm×n,使得式(3)成立:

式中:S=diag[σ1,σ2,…,σp],p=min(m,n)。矩陣S對角線上非零的元素即為矩陣A的奇異值,且σ1≥σ2≥…≥σp≥0。矩陣A由包含應力狀態(tài)信息的MBN 信號在時-頻域的可逆分解系數組成,因此由矩陣S對角線上非零元素組成的MBN 奇異值向量σ=[σ1,σ2,…,σp]與結構鋼表面應力之間存在聯(lián)系。

2.2 鐵磁材料的應力預測

圖5 為BP 神經網絡應力預測模型的建立過程。采用小波包變換對MBN 原始信號進行分解,以小波包變換系數作為時頻特征向量;利用不同尺度的小波包變換系數向量構造系數矩陣,并通過系數矩陣的奇異值分解實現小波包變換系數特征向量的降維;以奇異值向量作為BP 神經網絡的輸入樣本,實際應力值作為輸出樣本,設定初始參數進行訓練優(yōu)化,建立精確度較高的應力預測模型。

圖5 基于數據驅動的應力分布預測流程圖Fig.5 Flow chart of stress distribution prediction based on data-driven method

網絡層數和各層的神經元數、學習速率、學習算法、期望誤差、激活函數等是BP 神經網絡的關鍵參數。層數過多會使網絡過于復雜繁瑣,導致網絡訓練時長增加,出現過擬合趨勢,因此選擇由輸入層、隱含層和輸出層組成的3 層神經網絡。根據輸入數據的維度與輸出期望,設置輸入層神經元數為6,輸出層神經元數為1。學習速率決定權值的變化量,學習速率過大會導致系統(tǒng)不穩(wěn)定,反之會使訓練時間過長,收斂速度也隨之下降,但可使誤差維持在較低水平。因此通常采用較小的學習速率維持網絡的穩(wěn)定性,選取學習速率值為0.05。在網絡收斂速度快且無振蕩產生的前提下,訓練目標誤差應取較小值。改變目標誤差進行訓練,得到最佳目標誤差為0.000 65,最大訓練次數為5 000 次。

典型的神經元模型的輸出可表示為

式中:xt為第t個神經元的輸入;wt為第t個神經元的連接權重;θ為閾值;φ為激活函數或作用函數。BP 神經網絡通過反向傳播不斷調整每個MBN 特征向量的權重,使最佳誤差值控制在預期范圍內。

BP 神經網絡的輸入層、隱藏層、輸出層連接方式如圖6 所示,圖中Z[l]為第l層神經元輸入矩陣,為第l層中第u個神經元的激活函數輸出,可表示為

圖6 3 層BP 神經網絡結構圖Fig.6 Three-layer BP neural network structure diagram

式中:為網絡第(l-1)層中第v個神經元指向第u個神經元的連接權重;為第l層中第u個神經元的偏差;為第(l-1)層中第v個神經元的激活函數輸出。

對于多個MBN 奇異值向量同時輸入的情況,將式(5)換成對應的奇異值向量,輸入向量列數為m,每一列表示一個輸入樣本,多樣本輸入可表示為

式中:為第l層神經元連接權重矩陣;為l層神經元激活函數輸出矩陣;為第l層神經元偏差為第(l-1)層第m列神經元激活函數輸出。式(8)中每一列表示一個降維后的MBN 特征向量,對應一個MBN 奇異值向量。

3 試驗研究與結果分析

3.1 懸臂梁梯度應力標定試驗

對圖7(a)所示的結構鋼試塊進行懸臂梁梯度應力加載試驗,如圖7(b)所示。采用夾持工裝固定試塊左端,在試塊右端施加載荷,拾取懸臂梁上表面梯度應力狀態(tài)下的MBN 信號,每隔2 mm 測量并存儲一次MBN 信號。

圖7 懸臂梁梯度應力標定試驗Fig.7 Cantilever beam gradient stress calibration test

采用懸臂梁應力計算公式獲取每個數據采集點對應的應力值:

式中:s為試塊橫截面積;a為試塊厚度;MF為LF力矩,L、F分別為圖7(b)中數據采集點與載荷的距離、所施加的載荷。

由表4可知,本實驗所選用的模型顯著(P<0.0001),失擬項P=0.505>0.05,不顯著;模型的校正系數RAdj2=0.9911,相關系數R2=0.9678,說明該模型擬合程度良好,試驗誤差小,可以用此來分析和預測殼聚糖-明膠復合澄清劑的澄清作用。由表4回歸方程系數顯著性檢驗可知,殼聚糖對20%vol紅棗白蘭地透過率的影響極顯著(P<0.01),明膠對20%vol紅棗白蘭地透過率的影響極顯著(P<0.01);二次項AB、A2和B2也對透過率影響極顯著(P<0.01)。這表明各試驗因素對透過率的影響呈二次關系,且因素之間存在交互作用。

選取不同應力水平下的磁巴克豪森噪聲信號時域波形信號進行分析,如圖8 所示,σ為梯度應力最大值。在時域內,MBN 信號強度隨拉應力減小而減??;在頻域內,MBN 信號能量主要集中在8~75 kHz,該頻帶能量隨結構鋼表面拉應力減小而減小。綜上所述,當結構鋼表面拉應力發(fā)生變化時,MBN 信號在時、頻域內均發(fā)生明顯變化,可實現梯度應力變化規(guī)律的定性描述,但時頻分析無法定量描述這一規(guī)律。

圖8 不同應力水平下的MBN 信號時頻特征分析Fig.8 Time-frequency characteristics analysis of MBN signals at different stress levels

3.2 數據降維

圖8 的時域波形特征和頻譜分布特征隨應力水平的變化規(guī)律表明,MBN 信號的時頻特征可作為定量表征應力的特征指標。在時頻分析中,小波包變換系數反映了MBN 信號在時頻域的局部能量特征,可作為定量表征應力的特征指標。

通過式(1)的小波包分解獲取梯度應力的小波包變換系數特征向量,如圖9 所示,進一步按照式(2)構造小波包系數矩陣A。小波包系數矩陣的行向量反映MBN 信號在某一頻帶上的能量波動,激勵信號的信號能量主要集中在頻帶1~頻帶4,這表明頻帶1~頻帶4 的能量衰減較弱,能有效保留MBN 信號中包含的結構鋼表面應力特征信息。

圖9 小波包系數矩陣Fig.9 Wavelet packet coefficient matrix

小波包系數矩陣反映了MBN 信號的時頻冗余特征,直接作為特征向量將使得模型復雜度較高。采用奇異值分解對系數矩陣A進行數據降維,降低模型復雜度,可得到不同應力水平對應的低維度奇異值向量,不同應力水平奇異值向量變化規(guī)律如圖10 所示。

圖10 不同應力水平MBN 奇異值向量Fig.10 MBN singular value vector with different stress levels

圖10 表明MBN 奇異值隨著結構鋼表面拉應力的增大而增大,說明MBN 奇異值與梯度拉應力具有一定的相關性,且特征向量的維數減小為6。因此通過觀察不同應力水平下MBN 信號對應的奇異值向量曲線,可定性評估結構鋼表面應力水平,但實現定量預測還需根據奇異值向量建立定量預測模型。

3.3 非線性映射定量預測模型

數據降維后的MBN 特征向量與對應的實際應力值組成BP 神經網絡訓練集。根據預測效果對BP 神經網絡的各項參數,尤其是隱藏層節(jié)點數進行優(yōu)化。優(yōu)化前的隱藏層節(jié)點數太少,網絡無法建立復雜邊界,不能識別未學習的樣本,容錯性差,導致網絡輸出應力與實際應力的誤差較大,因此需對神經網絡的參數進行優(yōu)化。

在BP 神經網絡中,隱藏層節(jié)點數可根據經驗公式來確定:

圖11 為BP 神經網絡誤差收斂過程。迭代次數為43 次時,網絡最佳訓練誤差為0.000 51,小于既定訓練目標誤差,表明該網絡學習效率較高,達到了既定訓練標準。利用懸臂梁數據集對BP神經網絡訓練完成后,對BP 神經網絡優(yōu)化前后進行對比,如圖12 所示。圖12 表明優(yōu)化前的網絡輸出應力與對應的實際應力偏差較大,離散性較大,最大誤差達到了15 MPa,預測準確性較差;優(yōu)化后的最大誤差為1.5 MPa,預測效果較好。

圖11 迭代43 次時的最佳訓練誤差Fig.11 Best training error after iterating 43 times

圖12 BP 神經網絡優(yōu)化前后對比示意圖Fig.12 Comparison diagram of BP neural network before and after optimization

進一步采用MATLAB 內置BP 神經網絡工具箱繪制了BP 神經網絡應力預測模型優(yōu)化前后的回歸線,可直觀表示網絡訓練效果,如圖13 所示。與優(yōu)化前相比,優(yōu)化后的絕大多數訓練樣本的預測誤差在2%以內,回歸線的斜率更趨近于1,與圖14 所示的COMSOL 有限元仿真云圖進行對比,可看出神經網絡輸出應力變化趨勢與仿真結果一致,表明該模型具有較高的預測精度。

圖13 懸臂梁數據集的訓練效果Fig.13 Training effect of cantilever data set

圖14 懸臂梁梯度應力COMSOL 仿真云圖Fig.14 COMSOL simulation nephogram of gradient stress of cantilever beam

3.4 三點彎曲模型預測精度測試試驗

試驗采用點彎曲應力加載平臺施加垂直向下的載荷于結構鋼試塊中心區(qū)域,圖15 為三點彎曲應力加載平臺。探頭沿試塊下表面掃查,從數據采集起點開始每隔2 mm 測量并存儲1 次MBN信號,測量時保持傳感器與試件表面穩(wěn)定良好接觸。利用三點彎曲應力計算公式可得出試塊表面實際應力值:

圖15 三點彎曲應力加載平臺Fig.15 Three-point bending stress loading platform

式中:Mr為試塊上某點的彎矩;Iz為試塊上某點對軸的慣性矩(z軸垂直于紙面方向);d為載荷施加點到數據采集點的距離。

提取經數據降維處理后的MBN 奇異值向量作為輸入量,導入訓練完成的BP 神經網絡,獲取三點彎曲狀態(tài)下不同測量點對應的神經網絡輸出應力,并采用熱力圖進行三點彎曲下結構鋼試塊下表面應力分布成像,如圖16 所示。熱力圖表明三點彎曲狀態(tài)下結構鋼試塊下表面應力呈壓-拉-壓分布,試塊中心區(qū)域拉應力最大,從中心區(qū)域往試塊兩端呈梯度下降趨勢,這與式(11)的規(guī)律是一致的。

圖16 神經網絡輸出的三點彎曲試塊下表面應力分布成像結果Fig.16 Imaging results of stress distribution on lower surface of three-point bending specimen output by neural network

優(yōu)化后的BP 神經網絡模型對測試集樣本的預測結果如表1 和圖17 所示。優(yōu)化后的應力預測值與實際應力值的相對誤差為0.58% ~16.30%,大多數測試樣本的預測誤差均在10%誤差區(qū)間內,預測精度較高。

表1 三點彎曲試塊下表面應力BP 神經網絡預測模型測試效果Table 1 Test effect of BP neural network prediction model for lower surface stress of three-point bending test block

圖17 BP 神經網絡模型三點彎曲試塊下表面應力預測效果Fig.17 Prediction effect of BP neural network model on lower surface stress of threepoint bending specimen

提取MBN 信號均方根、均值、標準差、偏度、峰值平均值、峰值最大值等常規(guī)統(tǒng)計特征量作為BP神經網絡測試集輸入樣本,導出預測值并與MBN奇異值向量作為測試集輸入樣本的模型預測值以及實際值進行對比分析,結果如圖18所示。

圖18 常規(guī)統(tǒng)計特征量與MBN 奇異值預測效果對比Fig.18 Comparison of conventional eigenvalues and MBN singular value prediction effects

圖18 表明以常規(guī)統(tǒng)計特征量建立的模型預測值偏差較大,離散型較大。與常規(guī)統(tǒng)計特征量相比,以MBN 奇異值向量建立的模型預測值與實際值吻合度更高,精度更高。綜上所述,小波包變換系數構造的MBN 奇異值向量能有效反映結構鋼應力值變化的特征信息。以降維后的MBN 奇異值向量為BP 神經網絡的輸入樣本,通過優(yōu)化神經網絡結構參數實現鐵磁材料表面應力的高精度預測。同時,該方法在進一步豐富試驗數據、增加訓練樣本集數量時,能有效減小誤差,提高應力預測的可靠性。

3.5 模型對比分析

采用三點彎曲試驗的MBN 奇異值向量作為自變量,應力值作為因變量,建立多元線性回歸模型,與BP 神經網絡模型的三點彎曲試塊下表面應力識別效果進行對比。多元線性回歸模型自變量個數為6,因變量個數為1。設置好參數后,分別用多元線性回歸和BP 神經網絡對三點彎曲試塊下表面應力分布情況進行預測,應力預測值和實際值對比結果如圖19 所示。

圖19 模型預測效果對比Fig.19 Comparison of prediction effects of models

由圖19 的預測誤差對比可知,多元線性回歸模型的預測結果中,絕大多數樣本的預測值偏離實際值較大,離散性較大;BP 神經網絡模型的預測值與實際值吻合度較高,離散性較小。因此與多元線性回歸模型相比,BP 神經網絡模型預測結果更為精確,更適用于鐵磁材料表面應力的定量評估。其主要原因為多元線性回歸對應力與奇異值向量的復雜映射關系擬合能力較差,而3 層BP 神經網絡可實現以任意精度無限逼近任意函數。

4 結 論

1)通過小波包變換將MBN 信號時域波形轉換到時、頻域,提取MBN 信號時頻冗余特征,可定性反映不同應力水平。

2)利用不同頻帶的小波包變換系數特征向量構造了小波包變換系數矩陣。通過奇異值分解,減小了MBN 特征向量的維數,降低了BP 神經網絡模型復雜度。

3)開展了懸臂梁梯度應力標定試驗,采集了懸臂梁上表面梯度應力的MBN 信號,根據不同測量點對應的應力大小,建立了BP 神經網絡的訓練樣本集。

4)通過三點彎曲狀態(tài)下結構鋼試塊下表面梯度應力的MBN 信號采集試驗,建立了BP 神經網絡測試集,測試了數據驅動模型的預測精度。

5)采用小波包系數作為時頻冗余特征,經奇異值分解降維后,建立BP 神經網絡應力預測模型,使得該模型復雜度低、響應速度快,大部分樣本預測誤差低于10%,實現了基于數據驅動的鐵磁構件表面應力的定量評估。在實際工程應用中,結合機械臂自動化掃查、傳感器自適應仿形優(yōu)化、信號處理等技術,可實現鐵磁構件表面應力分布高精度成像。

猜你喜歡
信號模型
一半模型
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
重要模型『一線三等角』
完形填空二則
重尾非線性自回歸模型自加權M-估計的漸近分布
孩子停止長個的信號
3D打印中的模型分割與打包
基于LabVIEW的力加載信號采集與PID控制
一種基于極大似然估計的信號盲抽取算法
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 国产无套粉嫩白浆| 亚瑟天堂久久一区二区影院| 超清人妻系列无码专区| 香蕉eeww99国产在线观看| 91网址在线播放| 欧美成人aⅴ| 97视频在线精品国自产拍| 亚洲av成人无码网站在线观看| 国产一级妓女av网站| 精品无码日韩国产不卡av | 成人在线观看不卡| 久久精品国产一区二区小说| 免费观看男人免费桶女人视频| 成人午夜网址| 91国内视频在线观看| 人禽伦免费交视频网页播放| 日本高清视频在线www色| 亚洲AV无码久久精品色欲| 久久久四虎成人永久免费网站| 色网站在线免费观看| 中文字幕亚洲精品2页| 在线毛片免费| 91亚洲免费| 亚洲AⅤ综合在线欧美一区| 99草精品视频| 亚洲日本一本dvd高清| 免费看一级毛片波多结衣| 日韩精品少妇无码受不了| 免费看美女自慰的网站| 99久久精品国产自免费| 日本三区视频| 国产精品福利在线观看无码卡| 亚洲性网站| 日韩欧美高清视频| 亚洲天堂视频在线观看| 久久精品无码专区免费| 乱码国产乱码精品精在线播放| 国产成人亚洲精品无码电影| 国产成年女人特黄特色毛片免| 欧洲极品无码一区二区三区| 日韩国产精品无码一区二区三区 | 最新国产成人剧情在线播放| 99视频在线观看免费| 91在线免费公开视频| 久久综合亚洲色一区二区三区| 日韩福利在线观看| 亚洲午夜国产片在线观看| 波多野结衣一区二区三区88| 日日摸夜夜爽无码| 91视频国产高清| 四虎成人在线视频| 亚洲日韩久久综合中文字幕| 亚洲日韩国产精品无码专区| 免费国产高清视频| 欧美日韩午夜| 国产丝袜一区二区三区视频免下载| 视频二区亚洲精品| 狠狠色狠狠色综合久久第一次| 久久久国产精品免费视频| 麻豆AV网站免费进入| 亚洲精品第1页| 国产在线一区视频| 好久久免费视频高清| 女同久久精品国产99国| 中文字幕av无码不卡免费 | 欧美另类视频一区二区三区| 欧美成a人片在线观看| 日本午夜影院| 久久精品一卡日本电影| 亚洲精品在线观看91| 青青草欧美| 无码 在线 在线| 婷婷亚洲天堂| 日本道综合一本久久久88| 天天色天天综合网| 亚洲精品麻豆| 久久久久青草线综合超碰| 亚洲综合香蕉| 色综合五月| 中文字幕在线观看日本| 99精品视频播放| 2022精品国偷自产免费观看|