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

基于多模態特征子集選擇性集成建模的磨機負荷參數預測方法

2021-09-28 07:20:40柴天佑
自動化學報 2021年8期
關鍵詞:筒體模態特征

劉 卓 湯 健 柴天佑 余 文

磨礦過程難以檢測的磨機負荷參數(Mill load parameter,MLP),即料球比(Material to ball volume ratio,MBVR)、磨礦濃度(Pulp density,PD)和充填率(Ball charge volume ratio,CVR),與磨礦過程的產品質量和生產效率密切相關,對其進行實時檢測是實現選礦過程運行優化控制的關鍵因素之一[1-2].數據驅動的軟測量建模技術廣泛用于類似難測參數的推理估計[3].旋轉運行的磨機在筒體、軸承、研磨區域等不同位置所產生的機械振動/振聲信號,在產生機理、靈敏度和蘊含信息等方面存在差異性、冗余性與互補性[4-5].文獻[6]指出事物被體驗或表達存在多種特定方式(模態),可融合多模態的有價值信息實現分類或回歸.針對球磨機系統,蘊含著差異化MLP 信息的多模態機械振動/振聲信號頻譜特征常用于構建MLPF模型[7-9].

為構建具有可解釋性和較強泛化能力的MLPF模型,進行多模態頻譜特征的有效選擇是較為有效的策略.特征選擇算法能夠有效去除 “無關特征”與“冗余特征”,并確保重要特征不丟失[10].針對機械振動頻譜進行特征子集的選擇更具有價值[11].此外,不同MLP 與高維特征間的映射關系也呈現出差異性.

基于單個輸入特征與MLP 等難測參數間的相關系數能夠選擇線性相關特征,如文獻[12]結合多目標優化算法和相關系數進行微陣列數據的特征選擇,文獻[13]提出基于相關系數的多目標半監督特征選擇方法,文獻[14]提出基于熵的相關系數的特征聚類方法對特征子集進行快速聚類.針對基于相關系數的線性方法難以描述復雜非線性映射關系的缺點,互信息法可有效選擇與難測參數相關的非線性特征[15-16],如文獻[17-18]提出了基于個體最佳互信息和條件互信息的特征選擇方法.針對高維頻譜數據,如何自適應確定特征選擇閾值并進行有效的線性和非線性特征子集的選擇是待解決的開放問題.

在選擇包含不同數量特征的線性和非線性頻譜特征子集后,還需解決MLPF 的構建問題.通常,上述線性和非線性機械頻譜特征子集間存在冗余性和互補性.基于這些特征子集所構建的線性或非線性模型針對不同MLP 參數的預測性能也存在差異性.集成建模通過組合多個異質或同質子模型的輸出提高預測模型的穩定性和魯棒性,其難點問題是如何提高子模型間的多樣性.文獻[19]指出子模型多樣性的構造策略包括樣本空間重采樣、特征空間特征子集劃分或特征變換等,其中基于特征空間的構造策略能夠獲得更佳的預測性能.

基于較高頻率分辨率獲得的單/多尺度頻譜間具有較強的共線性.潛結構映射或偏最小二乘(Partial least squares,PLS)算法能夠提取低維潛在變量(Latent variable,LV)構建回歸模型,適合對高維頻譜數據建模[20-21].為提高MLPF 模型泛化性能,文獻[22]提出基于 “操作輸入特征”與 “采樣訓練樣本”的雙重維度集成構造策略,用以構建SEN MLPF 模型;該方法構建的模型能夠有效融合多源有價值信息,與運行專家感知MLP 的機制相類似,但所構建的MLPF 模型結構較為復雜.此外,上述方法屬于同質子模型集成,并且未對頻譜特征進行線性或非線性子集的選擇.面對小樣本,基于單尺度頻譜構建的MLPF 模型的泛化性能較好,但在模型可解釋性和洞悉研磨過程機理等方面存在欠缺.隨機權神經網絡(Random weight neural network,RWNN)是隱含層輸入權重隨機產生、輸出權重采用Moore-Penrose 廣義逆方法計算的單隱含層神經網絡[23-24],具有較快的學習速度.為克服工況漂移帶來的泛化性能下降問題,文獻[25-27]構建了基于更新樣本和遷移學習的MLPF 模型.為選擇更有價值的多模態機械信號,文獻[28] 提出了面向MLPF 的信號分析評估與優化組合算法,但其并未進行多模態特征子集的選擇.

在已有研究中,關于多模態特征選擇與融合的研究多偏向基于圖像的應用領域,如文獻[29]利用多核學習方法融合對圖像所提取的顏色、紋理、輪廓、深度等多模態特征進行分類檢測;文獻[30]提出融合語音信號和腦電信號的多模態情感識別;文獻[31]綜述多模態功能神經影像的多元機器學習融合方法;文獻[32]綜述了深度多模態學習,指出采用正則化方法對多模態融合結構進行學習和優化是研究熱點;文獻[33]提出采用多種深度神經網絡的多模態融合策略手勢識別.上述研究的多模態主要是從數據來源的視角進行,存在的問題是深度融合方法難以對所選擇多模態特征進行合理的物理解釋.

因此,基于多模態機械頻譜的MLPF 模型構建,需要解決以下2 個問題:1)如何進行線性特征和非線性特征子集的選擇;2)如何基于多模態特征子集構建差異性的集成子模型并進行有效選擇與合并以構建SEN 模型.綜上,本文提出了一種基于多模態特征子集SEN 建模的MLPF 方法.首先,對多模態機械信號進行時頻域變換得到高維頻譜數據;然后,采用相關系數法和互信息法對多模態頻譜變量進行線性和非線性特征子集的自適應選擇;接著,構建基于多模態特征的候選子模型后進行自適應選擇與合并,進而獲得基于SEN 的MLPF 模型.采用磨礦過程實驗球磨機的多模態機械振動信號仿真驗證了所提方法的有效性.

1 磨機機械信號的多模態特性分析

球磨機是依靠鋼球和礦石間的沖擊和研磨進行破碎的重型旋轉設備.磨機系統不同位置機械信號的產生機理如圖1 所示[28].

如圖1 所示,分層排列的鋼球隨球磨機旋轉以不同的沖擊力和周期被拋落或滑落,這些沖擊力相互疊加導致筒體振動;同時,磨機自身質量的不平衡和安裝偏置也引起筒體振動,這些疊加后的沖擊力導致筒體振動.顯然,磨機旋轉整周期的不同階段所蘊含的MLP 信息具有差異性.筒體振動經機械傳動設備的多級傳遞會濾掉高頻部分,其殘余振動導致磨機前后軸承座在水平和垂直方向的振動,再耦合其他來源振動后導致軸承座振動和.球磨機內部噪聲經連續反射后形成混合聲場,經磨機筒體傳出產生空氣噪聲,筒體振動的聲發射機理產生結構噪聲,此兩種噪聲是筒體近表面振聲的主要組成部分;此外,磨機研磨區域下方的振聲還包含來自其他磨機和鄰近設備的背景噪聲.

圖1 磨機系統不同位置機械信號的產生機理示意圖Fig.1 Generation mechanism of mechanical signals in different position of mill system

因此,球磨機系統不同位置測量的多模態機械信號所蘊含的MLP 信息具有冗余性和互補性.理論上,不同區域的頻譜具有其特定的物理含義.因此,有必要對多模態機械頻譜特征進行選擇性融合.

2 建模策略

本文提出了由時頻域變換、多模態線性/非線性特征子集選擇、多模態SEN 模型構建共3 個模塊組成的建模策略,如圖2 所示.

圖2 建模策略Fig.2 The proposed modeling strategy

圖2 中,J表示多模態機械信號的測量數量,Xj表示基于設計方案進行N次實驗所采集的長度為M的第j個模態的機械信號,Zj表示將Xj采用FFT(Fast Fourier transform)技術變換至頻域獲得的機械頻譜,表示基于Zj所提取的線性和非線性特征子集,表示基于線性特征子集構建的線性PLS 和非線性RWNN 子模型的預測輸出,表示基于非線性子集構建的線性PLS 和非線性RWNN子模型的預測輸出,y和表示MLP 及其預測輸出,Jsel和kfeasel表示MLPF 軟測量模型的集成尺寸和特征選擇系數.

3 建模算法

3.1 時頻域變換

以基于設計方案進行N次實驗所采集的長度為M的第j個模態的機械信號Xj為例,可用如下矩陣表示,

3.2 多模態線性/非線性特征子集選擇

對全部模態執行上述過程后,采用如下公式獲得線性特征選擇系數的上限和下限,

接著,計算第j個模態中最大和最小互信息值與均值的比值,如下所示,

對全部模態執行上述過程后,采用如下公式獲得非線性特征選擇系數的上限和下限,

進一步,采用如下公式獲得頻譜特征選擇系數的上限和下限,

本文中,對全部模態采用相同的特征選擇系數kfeasel.

按如下公式結合預測性能自適應確定Jfeasel個特征選擇系數的取值,

其中,kfeasel取為1 時,表示閾值為均值;表示用于計算Jfeasel個特征選擇系數的步長,采用如下公式獲得,

不同模態的線性特征和非線性特征的選擇閾值θfeasel采用如下公式獲得,

以第j個模態的第p個輸入特征為例,線性特征按如下規則進行選擇,

3.3 多模態SEN 模型構建

偏最小二乘(PLS)算法的目標是通過最大化輸入輸出數據間的協方差,將原始輸入特征空間的信息投影到由少數潛在變量(LV)組成新空間[11].以第jth個線性特征子集構建線性PLS子模型為例,其構建過程可表示為,

此處,需要采用優化算法從 4J個候選子模型中選擇Jsel個集成子模型(即SEN模型的集成尺寸)并對Jsel個集成子模型的預測輸出進行組合,得到最終基于SEN的MLPF 預測模型的輸出,即存在如下關系,

其中,fSEN(·) 表示對Jsel個集成子模型的預測輸出進行合并的算法.

針對上述問題,此處采用的策略是:首先選定用于合并集成子模型預測輸出的算法,然后以最小化SEN 模型的均方根相對誤差(Root mean square relative error,RMSRE)為準則,采用優化算法尋優fSEN(·)Jsel個集成子模型并對其輸出進行合并,得到集成尺寸為Jsel的SEN 模型.

用于對Jsel個集成子模型的預測輸出進行合并的算法fSEN(·) 主要包括計算集成子模型的加權系數法和采用線性、非線性回歸建模方法構建集成子模型輸出與SEN 模型輸出間的映射關系.其中,加權系數法包括簡單平均方法、自適應加權融合方法和誤差信息墑加權方法.本文采用自適應加權融合算法:

采用如下公式獲得SEN 模型的輸出,

4 實驗驗證

4.1 實驗數據描述

以基于實驗球磨機的高維筒體振動頻譜對MLPF進行建模以驗證本文所提方法.本實驗在直徑為602 mm 和長度為715 mm 的小型實驗磨機上進行,其中磨機筒體的旋轉速度為42 r/min.全部8 個模態機械信號的采樣頻率均為51 200 Hz,傳感器安裝位置和類型:1)固定在磨機筒體表面的2 個加速度傳感器;2)與磨機筒體表面相距2 mm 的2 個聲傳感器;3)位于磨機軸承座左側測量垂直振動、右側測量垂直和水平振動的3 個加速度傳感器;4)位于磨機研磨區域下方10 mm 的1 個聲傳感器.這些模態依次被標記為Ch1~Ch8,如圖3 所示[28].

圖3 實驗球磨機傳感器布置示意圖Fig.3 Layout of sensors for experimental ball mill

共進行了5 種工況下的527 次實驗:第1 次(B=292 kg,W=35 kg,M=25.5~174 kg);第2 次(B=340.69 kg,W=40 kg,M=29.7~170.1);第3 次(B=389.36 kg,W=40 kg,M=34.2~157.5 kg);第4 次(B=438.03 kg,W=35 kg,M=23.4~151.2 kg);第5 次(B=486.7 kg,W=40 kg,M=15.3~144.9 kg),其中,B、M、W分別代表鋼球、物料和水負荷.上述實驗均是在固定鋼球和水負荷,逐漸增加礦石負荷的情況進行的.將全部樣本中的4/5 用做訓練和驗證數據集,其余的用于模型測試.

4.2 實驗結果

4.2.1 時頻域轉換結果

首先,對時域信號進行濾波處理;然后,采用FFT技術將磨機運行中穩定旋轉周期的數據轉換至頻域,得到每個模態的多個旋轉周期的單尺度頻譜;最后,將這些穩定旋轉周期的譜數據進行平均,獲得最終維數為12 800 的建模頻譜.部分實驗的頻譜曲線詳見文獻[28].

4.2.2 特征子集選擇結果

基于317 個訓練數據,計算全部8 個模態的頻譜變量分別與三個磨機負荷參數(MBVR、CVR 和PD)間的相關系數和互信息值,此處以MBVR 磨機負荷參數的Ch1、Ch6 和Ch8 通道為例,如圖4 所示.

由圖4 和其他通道圖可知,基于相關系數和基于互信息的特征度量結果對不同模態機械頻譜特征具有差異性,并且采用不同閾值選擇的特征子集也具有差異性.

圖4 模態Ch1,Ch6,Ch8 的頻譜變量與MBVR 間的相關系數和互信息值Fig.4 Correlation coefficient and mutual information value between spectrum variable of mode Ch1,Ch6,Ch8 and MBVR

依據本文所提方法,基于相關系數值和互信息值的最小/最大值與各自均值的比值確定特征選擇系數的上限和下限,兩者相差越大表明頻譜特征間的差異越大.

本文以三個磨機負荷中的PD 模型為例,面向PD 的不同模態頻譜特征的特征選擇系數統計結果如下表所示.

由表1 可得到,相關系數閾值的下限和上限分別為0.8741 和1.057,互信息閾值的下限和上限分別為0.9228 和1.0599,進而得到PD 多模態特征選擇的下限和上限分別為0.9228 和1.057.

表1 面向PD 的不同模態頻譜特征的特征選擇系數統計表Table 1 Coefficients statistical table of different modal spectrum feature for PD

通過三個磨機負荷參數的不同模態頻譜特征的特征選擇系數統計表可知,不同模態針對不同磨機負荷參數選擇的特征子集具有差異性,基于這些特征子集構建線性和非線性子模型是非常必要的.

4.2.3 SEN 模型構建結果

本文中的線性建模方法為PLS,非線性建模方法為RWNN,采用驗證數據集確定前者的潛變量個數和后者的隱層節點個數.基于上述2 個特征子集和2 種建模方法,組合得到的候選子模型類型為4 種.為便于后文統計,對候選子模型的編碼采用如表2 所示的形式.

表2 中:在 “子模型特點”列中,前項表示特征類型,后項表示子模型類型,相應的 “lin”和 “non-lin”分別表示線性和非線性;在 “子模型名稱”列中,“Corr”和 “Mi”分別表示相關系數和互信息.

表2 候選子模型編碼Table 2 Coding of candidate sub-models

針對不同MLP,取特征選擇系數的數量Jfeasel均為10.采用不同特征選擇系數所構建的SEN 模型的預測誤差和所選擇的子模型編號如表3 所示.

由表3 可知:

表3 不同特征選擇系數時所構建的SEN 模型的預測誤差和所選擇的子模型編號Table 3 Prediction error of SEN model with different feature selection coefficients and selected sub-model number

1) MBVR 模型的最小預測誤差為0.04404,集成尺寸為8,其集成子模型集合為{25 17 18 27 22 19 30 24},子模型的預測誤差按照由高到低排序;結合表2 可知,這些集成子模型的具體含義為Mi-RWNN-Ch1、Corr-RWNN-Ch1、Corr-RWNN-Ch2、Mi-RWNN-Ch3、Corr-RWNN-Ch6、Corr-RWNNCh3、Mi-RWNN-Ch6、Corr-RWNN-Ch8;可見,線性特征-非線性子模型為3 個,非線性特征-線性子模型為5 個,表明基于相關系數構建的RWNN子模型性能更佳;模態通道為Ch1、Ch2、Ch3、Ch6和Ch8,分別代表筒體振動、筒體振動、筒體振聲、軸承振動和研磨區域振聲,其中Corr-RWNN-Ch8具有最小的預測誤差,表明不同模態特征蘊含的MBVR 信息具有差異性.

2) PD 模型的最小預測誤差為0.01452,集成尺寸為7,其集成子模型集合為{22 14 24 32 26 19 30},子模型的預測誤差按照由高到低排序;結合表2 可知,這些集成子模型的具體含義為Corr-RWNN-Ch6、Mi-PLS-Ch6、Corr-RWNN-Ch8、Mi-RWNN-Ch8、Mi-RWNN-Ch2、Corr-RWNN-Ch3、Mi-RWNNCh6;可見,線性特征-非線性子模型為3 個,非線性特征-非線性子模型為3 個,線性特征-線性子模型為3 個,表明基于互信息提取特征子集構建的RWNN 子模型性能更佳;模態通道為Ch6、Ch8、Ch2 和Ch3,分別代表軸承振動、研磨區域振聲、筒體振動和筒體振聲,其中Mi-RWNN-Ch6 具有最小的預測誤差,表明不同模態特征蘊含的PD 信息具有差異性,并且主要存在機械振動模態頻譜特征中.

3) CVR 模型的最小預測誤差為0.009280,集成尺寸為7,其集成子模型集合為{28 18 26 27 19 22 30},子模型的預測誤差按照由高到低排序;結合表2可知,這些集成子模型的具體含義為Mi-RWNNCh4、Corr-RWNN-Ch2、Mi-RWNN-Ch2、Mi-RWNN-Ch3、Corr-RWNN-Ch3、Corr-RWNN-Ch6、Mi-RWNN-Ch6;可見,線性特征-非線性子模型為3 個,非線性特征-非線性子模型為4 個,表明基于相關系數和互信息提取特征構建的RWNN 子模型性能能夠互補;模態通道為Ch4、Ch2、Ch3 和Ch6,分別代表筒體表面振聲、筒體表面振動、筒體表面振聲和軸承振動,其中Mi-RWNN-Ch6 具有最小的預測誤差,表明不同模態特征蘊含的CVR信息具有差異性,主要存在機械振動和筒體表面振聲頻譜特征中.

4.2.4 實驗結果比較

本文采用的是8 個模態機械信號數據,數據不同于其他已有磨機負荷預測方法.因此,此處進行單通道和本文所提方法的實驗結果比較.

由表4 和表5 可知:

表4 磨機負荷參數各通道與多模態特征子集選擇性集成模型的測試誤差比較Table 4 Comparison of test errors between various channels of mill load parameters and multi-modal feature subset SEN model

表5 磨機負荷參數各通道與多模態特征子集選擇性集成模型的平均測試誤差比較Table 5 Average test errors comparison of the various channels of mill load parameters and multi-modal feature subset SEN model

1) 本文所提方法針對3 個磨機負荷參數的平均預測誤差為0.02260,小于基于單通道的平均預測誤差,其中:Ch4 和Ch5 的預測誤差最大,分別為0.1019 和0.1056,其他5 個通道的預測誤差在0.0501~0.0774 之間.因此,本文所提方法的預測性能比基于單通道的方法提升了1 倍,表明了所提方法的有效性.

2) 不同通道的平均測試誤差間具有差異性,其中:筒體振動Ch1 和Ch2 的預測誤差差別較小,筒體振聲Ch3 和Ch4 的預測誤差相差2 倍,筒體左側軸承振動信號Ch5 比右側軸承振動信號Ch6和Ch7 的預測誤差大2 倍,研磨區域振聲Ch8 的預測誤差較小,這表明磨機不同位置的振動和振聲信號中包含著不同的磨機負荷參數信息,采用不同通道的多模態信號有助于實現磨機負荷參數的準確測量;此外,Ch4 和Ch5 預測誤差遠高于其他通道,表明這兩個通道可能存在故障.

3) 同一MLP 在不同機械通道上的測試誤差具有差異性,其中:針對MBVR,Ch8 的測試誤差最小,僅為0.08090,而Ch4 的最大,為0.2001,兩者相差2 倍.針對PD,Ch6 的測試誤差最小,為0.02431,Ch5 的最大,為0.08122,相差近4 倍;針對CVR,Ch6 的測試誤差最小,為0.01630,Ch1 的最大為0.04930,相差3 倍.此外,針對PD 和CVR,Ch6 通道均具有最佳的測試誤差.上述結果表明,針對同一個MLP 而言,各模態信號特征之間的確存在冗余性和互補性,進行多模態信號頻譜特征的選擇與融合是合理的.

4) 同一模態(通道)針對不同MLP 的測試誤差具有差異性:針對MBVR 的預測誤差均高于PD和CVR,如Ch3 的MBVR 預測誤差為0.09020,PD 和CVR 預測誤差分別為0.03111 和0.02922,相差近3 倍,說明各個通道對不同磨機負荷參數的敏感度不同;但是,這也可能與所采用的實驗磨機的特性是相關的.顯然,這需要進行更為深入的實驗研究.

5 結論

針對多模態高維機械頻譜數據輸入特征與MLP 間的可解釋映射模型難以構建的難題,本文提出了基于多模態特征異質模型集成的MLPF 方法.主要貢獻表現在:能夠依據多模態頻譜數據特性進行線性特征子集和非線性特征子集的自適應選擇,提出構建線性特征線性子模型、線性特征非線性子模型、非線性特征線性子模型、非線性特征非線性子模型的策略增強集成子模型間的差異性,提出了多模態頻譜子模型的自適應選擇與合并方法.通過磨礦過程實驗球磨機的高維機械振動和振聲頻譜數據仿真驗證了所提方法的有效性.

猜你喜歡
筒體模態特征
b型管板與筒體溫差應力的分析計算和評定
化工管理(2021年7期)2021-05-13 00:46:04
回轉窯筒體對接操作方法
水泥技術(2021年2期)2021-04-20 12:37:26
一種臥式筒體糞污發酵裝置的筒體設計與分析
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
抓住特征巧觀察
國內多模態教學研究回顧與展望
球磨機筒體鑄鋼端蓋裂紋的處理方法
基于HHT和Prony算法的電力系統低頻振蕩模態識別
由單個模態構造對稱簡支梁的抗彎剛度
計算物理(2014年2期)2014-03-11 17:01:39
主站蜘蛛池模板: 黄色三级网站免费| 欧美专区在线观看| 精品伊人久久久久7777人| 国产精品成人啪精品视频| 玩两个丰满老熟女久久网| 青青草国产精品久久久久| 久久窝窝国产精品午夜看片| 日韩精品一区二区三区大桥未久| 久久无码高潮喷水| 国产第一页免费浮力影院| 久久精品日日躁夜夜躁欧美| 黄片一区二区三区| 欧亚日韩Av| 成年人视频一区二区| 久久香蕉国产线看观| 91精品啪在线观看国产| 免费在线国产一区二区三区精品 | 2022精品国偷自产免费观看| 老色鬼欧美精品| 波多野吉衣一区二区三区av| 国产午夜不卡| 美女国产在线| 欧美日韩另类在线| AV在线麻免费观看网站| 中文字幕亚洲乱码熟女1区2区| 中文成人在线| 为你提供最新久久精品久久综合| 婷婷色婷婷| 国产理论精品| 亚洲国产中文精品va在线播放| 中日韩欧亚无码视频| 国产人人射| 日韩高清欧美| 夜夜操国产| 国产欧美视频综合二区| 福利在线一区| 亚洲啪啪网| 亚洲aⅴ天堂| 午夜国产在线观看| 日韩最新中文字幕| 成人伊人色一区二区三区| 国产人在线成免费视频| 国产成人精品免费av| 亚洲Aⅴ无码专区在线观看q| 免费av一区二区三区在线| 欧美激情伊人| 精品久久久久无码| 91久久青青草原精品国产| 国产成人三级在线观看视频| 国产主播一区二区三区| 欧美在线视频a| 久久熟女AV| 日韩一级二级三级| 中文字幕日韩视频欧美一区| 日日拍夜夜操| 伊人无码视屏| 国产午夜精品一区二区三| 呦女亚洲一区精品| 国产丝袜无码一区二区视频| 亚洲天堂啪啪| 精品国产乱码久久久久久一区二区| 狠狠做深爱婷婷久久一区| 亚洲人成网站观看在线观看| 国产亚洲成AⅤ人片在线观看| 久久精品免费看一| 99r在线精品视频在线播放 | 国产一级视频久久| 成人夜夜嗨| 日韩精品无码免费专网站| 欧美精品在线免费| 免费人成黄页在线观看国产| 97精品久久久大香线焦| 日本尹人综合香蕉在线观看| 午夜不卡视频| 亚洲中文字幕久久精品无码一区| www精品久久| 91在线中文| 精品無碼一區在線觀看 | 国产亚卅精品无码| 亚欧美国产综合| 毛片在线看网站| 国产精品视频久|