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

基于近紅外光譜技術的腐竹脂肪定量分析

2014-02-27 08:39:50王加華王一方韓東海
食品科學 2014年18期
關鍵詞:信息模型

王加華,王 軍,王一方,韓東海*

(1.許昌學院食品與生物工程學院,河南 許昌 461000;2.中國農業大學食品科學與營養工程學院,北京 100083;3.許昌市食品藥品監督管理局,河南 許昌 461000)

基于近紅外光譜技術的腐竹脂肪定量分析

王加華1,王 軍1,王一方2,3,韓東海2,*

(1.許昌學院食品與生物工程學院,河南 許昌 461000;2.中國農業大學食品科學與營養工程學院,北京 100083;3.許昌市食品藥品監督管理局,河南 許昌 461000)

采用近紅外光譜技術結合化學計量學方法,建立腐竹脂肪含量的快速分析方法。收集不同生產線、不同時間的腐竹樣本180 份,利用積分球附件采集漫反射光譜(4 000~10 000 cm-1)。為消除顆粒散射影響和光譜基線漂移,二階導數和卷積平滑用于光譜預處理。采用反向區間偏最小二乘法、組合區間偏最小二乘法、搜索組合移動窗口偏最小二乘法和遺傳偏最小二乘法優化建模變量,最終構建了定量預測模型。結果顯示,4 種方法均可有效地提取信息變量、降低模型維度、提高預測性能;遺傳偏最小二乘法一次優選獲得143 個變量,構建的模型性能最佳,其校正相關系數、校正均方根誤差、預測相關系數、預測均方根誤差分別為0.96、0.95、0.92和1.17。研究表明,經過信息變量提取后所構建的近紅外模型簡單、預測精度高,可用于腐竹脂肪含量的日常監測。

近紅外光譜;腐竹;脂肪;變量提取;定量分析

腐竹是我國典型的傳統干制豆制品,生產歷史悠久,用傳統加工工藝,反映地方和民族特色,它是一個民族長期適應的自然選擇。腐竹是豆漿中的蛋白質分子在變性過程中與脂肪分子相聚合而形成的薄膜(蛋白質-脂質膜),具有良好的風味性、營養性、健康性和安全性,能為人體提供均衡能量,長期以來深受人們的喜愛。

國內外學者大多關注腐竹成膜理論[1]、工藝優化[2-4]、膜結構和特性[5-7]等方面,而對于成品腐竹的品質檢測研究較少。脂肪是腐竹重要的營養成分之一,且脂肪含量

越高,膜的透水率越低、阻水性越強,防腐效果越好,貨架期就越長[8],因此腐竹脂肪的測定對于質量評級、生產控制具有重要意義。目前,腐竹脂肪測定采用常規測試方法,效率低、費時費力,不能滿足腐竹規模化生產的需求。因此,對于腐竹脂肪檢測,迫切需求一種快速、可實施在線監控的分析技術。

近紅外光譜記錄的是含氫基團X—H(X:C、N、O)單個化學鍵的基頻振動的倍頻和合頻信息,動植物性食品的成分大多由這些基團構成,基團的吸收頻譜表征了這些成分的化學結構和含量,因此近紅外光譜分析技術已廣泛用于果蔬[9-10]、畜產食品[11-12]、水產品[13-14]、糧油[15-16]、茶葉[17-18]、藥食材料[19-20]等品質指標檢測。腐竹的紅外光譜受含氫基團的重疊主導,信息豐富,為近紅外分析提供了理論基礎,但對于腐竹品質近紅外檢測應用尚未見報道。

本實驗以傳統食品腐竹為研究對象,比較研究反向區間偏最小二乘法(backward interval partial least squares,BiPLS)、組合區間偏最小二乘法(synergy interval partial least squares,SiPLS)、搜索組合移動窗口偏最小二乘法(searching combination moving window partial least squares,SCMWPLS)和遺傳偏最小二乘法(genetic algorithms partial least squares,GAPLS)優化腐竹脂肪信息變量的效果,并構建定量分析模型,以期為實現腐竹脂肪含量的綠色、快速檢測提供一定的參考依據。

1 材料與方法

1.1 材料與試劑

腐竹樣品 許昌某豆制品加工企業。

為獲取代表性樣品,分別在a、b、c、d四條不同生產線上,連續15 d收集生產的成品,在每條線生產的初、中和末3 個時期各收集3 個樣品,共有180 份獨立樣品。樣品密封于封口袋,包裝后運回實驗室,并分別標號后備用。

1.2 儀器與設備

ANTARISⅡ型傅里葉變換近紅外光譜儀(配備積分球采集附件、InGaAs檢測器) 美國Thermo Fisher Nicolet公司;R254S型索氏提取器 德國Behr實驗儀器設備公司;AUY220型電子分析天平(精度0.1 mg)日本島津公司;HB-DK-S26型電熱恒溫水浴鍋 北京恒奧德儀器儀表有限公司;DHG-9023A型鼓風烘箱 沈陽林頻實驗設備有限公司;JYL-C012型料理機 九陽股份有限公司。

1.3 方法

1.3.1 光譜采集

腐竹樣品經粉碎后,過18 目篩后加入到樣品杯,采集腐竹的積分球漫反射光譜。為獲取代表性的腐竹光譜,采用自動旋轉樣品杯附件采集,采集光譜范圍4 000~10 000 cm-1,分辨率8 cm-1,旋轉杯偏心距為8 mm,掃描32 次取平均。

1.3.2 脂肪化學值測定

采用索氏抽提法測定腐竹脂肪含量,方法參照GB/T 5009.6—2003《食品中脂肪的測定》。

1.3.3 化學計量學方法

BiPLS和SiPLS是N?rgaard[21]提出的區間偏最小二乘(interval partial least squares,iPLS)法的改進和演化。BiPLS是將整個光譜分割成k 個等寬子區間,分別計算各子區間的交互驗證均方根誤差(root mean square error of cross-validation,RMSECV)值,然后逐步去除RMSECV最大的區間i(i=0、1、2…k),在剩余的k-i區間上建立偏最小二乘法(partial least squares,PLS)模型,并給出相應的RMSECV值。當RMSECV值最小時所對應的多個區間即為所優化的組合區間。SiPLS是在各子區間上,計算所有可能的j(2≤j≤h)個子區間組合模型,依RMSECV值大小給出各個組合區間,當RMSECV值最小時該組合區間為最優區間組合。

SCMWPLS[22]是基于移動窗口偏最小二乘法(moving window partial least squares,MWPLS)基本原理的演化,計算步驟簡述如下:1)在給定最大窗口p下,計算1~p所有窗口下的MWPLS,獲取不同寬度的信息區間,在所得信息區間里選擇殘差最小的區間作為基礎信息區間A;2)以基礎信息區間A為基準,在剩余信息區間(除基礎信息區間A)里依次選擇單個區間與信息區間A進行組合,并計算殘差,殘差最小時,得到最優組合信息區間A和B;3)再以信息區間A和B為基準,重復步驟2),得到最佳信息區間組合信息區間A+B+C;4)按上述方法重復計算,至組合計算次數結束。結果輸出每步組合下的信息區間和殘差。最小殘差信息組合即為最優結果。

GAPLS[23]是引用生物界物種競爭選擇的進化機制,以適應度函數為依據,通過對群體中個體施加遺傳操作,如選擇、交叉、變異,來實現群體內個體結構重組的迭代優化。采用遺傳算法選擇特征變量,由于每次運行初始參數不同,如初始群體選擇,交叉變異位置等,輸出的0~1二進制編碼各異,因此,獨立運行100 次,以RMSECV為目標函數選取出現頻率較高的變量。

BiPLS、SiPLS、SCMWPLS和GAPLS程序均在MATLAB V7.0中實現。

2 結果與分析

2.1 光譜特征及預處理

腐竹是蛋白質分子在變性過程中與多糖和脂肪通

過分子間的相互作用而形成的可食性膜,具有多孔網絡結構,光譜信息豐富,且相互重疊。腐竹原始積分球漫反射光譜如圖1a所示,光譜形狀相似,在5 155 cm-1和6 890 cm-1處為水分子O—H伸縮和HOH彎曲振動的組合頻,8 310 cm-1處為水分子O—H比較弱的組合頻吸收。受顆粒散射影響,光譜基線漂移嚴重,二階導數可以消除基線漂移,且可以放大光譜信號,圖1b為二階導數(2D)和卷積平滑(7點3階多項式S-G平滑)處理后的光譜,處理后消除了基線漂移,4 230~6 080 cm-1區間信號更加豐富,而且也使6 3 2 9~6 8 2 8、8 060~8 650 cm-1區間的光譜差異凸顯。

圖1 腐竹原始光譜(a)和二階導數平滑處理后光譜(b)Fig.1 Spectra of yuba obtained from original data (a) and 2-Der with S-G smothering pretreatment (b)

在建立校正模型前,采用TQ軟件程序V8.0,在95%置信度下進行Chauvenet檢驗(置信度95%),剔除異常光譜。剩下的樣本依據化學值大小排序,依3∶1左右分為校正集和預測集,化學值的最大和最小樣本歸為校正集,統計參數如表1所示。

表1 腐竹校正集和預測集樣品化學值統計表Table 1 Statistics of fat content in yuba for calibration and prediction data sets

2.2 BiPLS模型

考察區間分割數對選擇結果及模型的影響,將整條光譜(1 557 個數據點)分為k個子區間(k=10~50,間隔5),在每種分割情況下運行BiPLS程序。采用留一法計算RMSECV,最大因子數設為10。

當分割數為40時,BiPLS所得RMSECV值最小為1.45,入選子區間為[8 31 6 21 34 4 13 20 26],所對應的信息區間分別是4 451~4 597、4 752~4 898、5 053~5 199、5 805~5 951、6 858~7 155、7 760~7 907、8 512~ 8 659 cm-1和8 964~9 110 cm-1,如圖2上方所示。BiPLS共選擇了351 個數據點,占全譜的25.5%。區間變量表征了腐竹脂肪分子振動信息,如4 587、4 545 cm-1附近是C=O伸縮振動合頻吸收,5 865、5 900 cm-1附件主要是甲基C—H伸縮振動一倍頻吸收,7 092 cm-1為油脂中O—H一倍頻吸收,5 128 cm-1和8 622 cm-1是C=O伸縮振動的二倍頻和四倍頻吸收,8 547 cm-1為不飽和脂肪酸中烯烴的C—H二倍頻吸收。

圖2 BiPLS和SiPLS信息區間選擇結果Fig.2 Informative regions obtained by BiPLS and SiPLS procedures

在上述信息區間內建立BiPLS模型,校正和預測結果如表3所示,其校正相關系數(R)和校正均方根誤差(root mean square error of calibration,RMSEC)分別為0.90、1.48;其預測相關系數(r)、相對預測均方根誤差(relative prediction mean square error,RMSEP)分別為0.89和1.53,優于全譜PLS模型(r=0.86,RMSEP=1.68),但數據點顯著減少。

2.3 SiPLS模型

在使用SiPLS程序時,不同的子區間數和組合數對輸出結果均有影響。在本實驗中,也將光譜分為k個子區間(k=10~50,間隔5),在不同分割數下,分別計算不同組合數(1~4)下的結果。

當分割數為40時,入選子區間為[8 20 28 31],所得RMSECV值最小為1.56,此時所選擇的信息區間是5 053~5 199、6 858~7 004、8 061~8 208 cm-1和8 512~8 659 cm-1,如圖2下方所示。SiPLS共選擇了156 個數據點,占全譜的10%。SiPLS選擇波段中5 053~5 199、6 858~7 004 cm-1和8 512~8 659 cm-1與BiPLS選擇結果相同,反映了腐竹脂肪特征吸收;而在8 208 cm-1附近主要反映了乙基C—H二倍頻吸收。在信息區間建立模型SiPLS,預測結果如表3所示,其r、RMSEP分別為0.87和1.60,其模型性能接近BiPLS模型,但是建模數據進一步減少到156 個。

2.4 SCMWPLS模型

SCMWPLS運行參數如下:基礎區間計算窗口大小為k(k=4~300,間隔1),次級區間窗寬為2,組合次數為100,最大因子數為10。

根據上述方法運行SCMWPLS程序,當k=201時,所得RMSECV值最小為1.54,對應基礎信息區間為7 328~8 100 cm-1,主要表征的是C—H的二倍頻和合頻吸收。在此基礎區間上,SCMWPLS進行最優搜索組合,最終RMSECV最小達到1.48,所獲得的次級信息區間分別為5 053~5 151、8 613~8 624 cm-1,共有29 個數據點,5 128 cm-1附近是—CO2R官能團的C=O伸縮振動二倍頻吸收,而在8 622附近是C=O伸縮振動四倍頻吸收。SCMWPLS共選擇了230 個數據點,占全譜的14.8%,在上述信息變量下建立SCMWPLS模型,預測結果如表3所示,其r、RMSEP分別為0.90和1.49,其模型預測性能略優于BiPLS和SiPLS模型。

2.5 GAPLS模型

當光譜變量數超過200時,運行GAPLS可能導致過擬合的風險[21]。本實驗在此作兩步處理:1)將整條光譜數據分割成8 個子區間,在每個子區間下運行GAPLS,將每個子區間下的最小RMSECV值對應的變量選出,重新構造一個光譜矩陣M,在新矩陣M下建立GAPLS1模型;2)在光譜矩陣M下,再運行GAPLS,二次選擇信息變量,并建立GAPLS2模型。

GAPLS運行參數分別為:種群大小30,入選變量數最大值30,變異概率0.01,交叉概率0.5,最大因子數10。獨立運行100 次,計算每個數據點標識為“1”的概率。入選變量以RMSECV為目標函數,當滿足F檢驗(P<0.1)時為最佳結果。

首先將1 557個數據點分割成8 個子區間,在每個子區間下運行GAPLS,輸出最小RMSECV值和對應的變量數,如表2所示,在第6個子區間(7 837~8 601 cm-1)內,RMSECV值最低,且模型具有較低的維度,表明在此區間變量表征了腐竹脂肪分子振動信息。

表2 不同子區間下GAPLS選擇結果Table 2 Results of GAPLS in different intervals

GAPLS一次優化共選擇了143 個數據點,占全譜的9.2%,采用上述143 個變量重新構造一個光譜矩陣M(圖3),新矩陣光譜差異明顯,排除了較多的非目標信息變量。在矩陣M下建立GAPLS1模型,其R、RMSECV分別為0.96和0.95,采用外部44 個樣品檢驗模型預測性能,散點圖如圖4所示,預測結果r、RMSEP分別為0.92和1.17(表3),顯著優于全譜PLS、BiPLS和SiPLS模型,且變量數最小。

圖3 經分段處理后GAPLS重構光譜矩陣Fig.3 Reconstruction spectrum matrix obtained by GAPLS after segmented treatment

圖4 GAPLS1模型預測結果散點圖Fig.4 Plots of measured firmness value vs. GAPLS1 predicted value

圖5 GA-PLS程序運行結果(a)變量選擇概率(b)不同變量數下的RMSSEECCVV值Fig.5 Variables selections accomplished by GA-PLS

在新矩陣M下再次運行GAPLS選擇變量,入選變量數與RMSECV值關系如圖5b所示,當有25 個變量入選時,即圖5a中水平線為所選擇的變量,RMSECV值達到最小(1.564)。由表3可得,GAPLS二次所選的25 個變量中有15 個變量來源于表2中的第6個子區間,進一步說

明了該區間包含了較多的腐竹脂肪信息變量,主要是脂肪烴的C—H二倍頻吸收,與表2所示該區間RMSECV值最小相一致。采用所選25 個變量構建模型GAPLS2,預測結果如表3所示,其r、EMSEP分別為0.87和1.60,與BiPLS和SiPLS相當,但數據使用量僅為全譜的1.6%。

表3 腐竹脂肪的不同PLS模型及性能評價結果Table 3 Results of PLS modeling for fat content of yuba and evaluation of the performance of their corresponding models

3 3 結 論

應用積分球附件采集腐竹樣品近紅外漫反射光譜,經過導數和平滑處理,以消除顆粒散射影響和基線漂移,探討BiPLS、SiPLS、SCMWPLS和GAPLS四種方法優化腐竹脂肪信息變量的有效性,并構建預測模型。結果表明經變量優化后,BiPLS、SiPLS、SCMWPLS和GAPLS四種方法的入選變量數分別為351、156、230 個和143 個,分別占全譜的變量數(1 557 個變量數)的22.5%、10%、14.8%、9.2%。用GAPLS一次優選變量(143 個數據點)所構建的GAPLS1模型最優,其R、RMSECV、r、RMSEP分別為0.96、0.95、0.92和1.17;在重構光譜矩陣基礎上,進行GAPLS二次優化獲得25 個變量,所構建的GAPLS2模型性能與全譜相當。結果表明:近紅外光譜技術結合化學計量學方法可用于構建腐竹脂肪定量模型,信息變量提取可簡化模型,提高預測性能,同時顯示了近紅外光譜技術在腐竹營養指標快速評價方面的應用潛力。

[1] CHEN Y, YAMAGUCHI S, ONO T. Mechanism of the chemical composition changes of yuba prepared by a laboratory processing method[J]. Journal of Agricultural and Food Chemistry, 2009, 57(9): 3831-3836.

[2] LONG Lei, HAN Zhi, ZHANG Xiujin, et al. Effects of different heating methods on the production of protein-lipid film[J]. Journal of Food Engineering, 2007, 82(3): 292-297.

[3] MARIA B P, JOHN M K. Drying temperature effect on water vapor permeability and mechanical properties of whey protein-lipid emulsion films[J]. Journal of Agricultural and Food Chemistry, 2000, 48(7): 2687-2692.

[4] 謝麗燕, 林瑩, 譚瑤瑤, 等. 正交試驗優化傳統腐竹制作工藝[J]. 食品科學, 2014, 35(2): 36-40.

[5] ZHELUDEVA S, NOVIKOVA N, STEPINA N, et al. Molecular organization in protein-lipid film on the water surface studied by X-ray standing wave measurements under total external reflection[J]. Spectrochimica Acta Part B, 2008, 63(12): 1399-1403.

[6] KOKOSZKA S, DEBEAUFORT F, LENARTA A, et al. Liquid and vapour water transfer through whey protein/lipid emulsion films[J]. Journal of the Science of Food and Agriculture, 2010, 90(5): 1673-1680.

[7] CHEN Y, ONO T. The mechanisms for yuba formation and its stable lipid[J]. Journal of Agricultural and Food Chemistry, 2010, 58(10): 6485-6489.

[8] 歐錦強, 王興國, 金青哲. 大豆組分對腐竹性能的影響[J]. 中國油脂, 2005, 30(2): 37-40.

[9] NICOLA? B M, BEULLENSA K, BOBELYNA E, et al. Nondestructive measurement of fruit and vegetable quality by means of NIR spectroscopy: a review[J]. Postharvest Biology and Technology, 2007, 46(2): 99-118.

[10] 劉燕德, 高榮杰, 孫旭東. 便攜式水果內部品質近紅外檢測儀研究進展[J]. 光譜學與光譜分析, 2010, 30(10): 2874-2878.

[11] 樊玉霞. 豬肉肉糜品質與安全可見/近紅外光譜快速檢測方法的實驗研究[D]. 杭州: 浙江大學, 2011.

[12] 王田子, 鄭麗敏, 田立軍, 等. 近紅外在乳及乳制品質量檢測中的研究進展[J]. 光譜學與光譜分析, 2010, 30(12): 3208-3212.

[13] 徐文杰, 李俊杰, 賈丹, 等. 近紅外光譜技術分析草魚營養成分[J].食品科學, 2013, 34(20): 161-164.

[14] 徐文杰, 劉茹, 洪響聲, 等. 基于近紅外光譜技術的淡水魚品種快速鑒別[J]. 農業工程學報, 2014, 30(1): 253-258.

[15] 王加華, 王一方, 屈凌波. 糧食品質近紅外光譜無損檢測研究進展[J].河南工業大學學報: 自然科學版, 2011, 32(6): 80-87.

[16] 劉建學, 劉珊珊. 芝麻油真偽快速檢測方法研究[J]. 中國糧油學報, 2012, 27(12): 116-121.

[17] 張龍, 潘家榮, 朱誠. 基于近紅外光譜的不同發酵類型茶葉判別[J].食品科學, 2012, 33(20): 149-152.

[18] 寧井銘, 宛曉春, 張正竹, 等. 近紅外光譜技術結合人工神經網絡判別普洱茶發酵程度[J]. 農業工程學報, 2013, 29(6): 255-260.

[19] 湯麗華, 劉敦華. 基于近紅外光譜的枸杞化學成分定量分析[J]. 現代食品科技, 2013, 29(9): 2306-2310.

[20] 趙羚志. 短波近紅外光譜技術結合人工神經網絡用于藥物無損定量分析的研究[D]. 長春: 吉林大學, 2009.

[21] N?RGAARD L, SAUDLAND A, WAGNER J, et al. Interval partial least squares regression (iPLS): a comparative chemometric study with an example from near-infrared spectroscopy[J]. Applied Spectroscopy, 2000, 54(3): 413-419.

[22] DU Yiping, LIANG Yizeng, JIANG Jianhui, et al. Spectral regions selection to improve prediction ability of PLS models by changeable size moving window partial least squares and searching combination moving window partial least squares[J]. Analytica Chimica Acta, 2004, 501(2): 183-191.

[23] LEARDI R. Application of genetic algorithm-PLS for feature selection in spectral data sets[J]. Journal of Chemometrics, 2000, 14(5): 643-655.

Determination of Fat Content in Yuba by Near Infrared Spectroscopy and Chemometrics

WANG Jia-hua1, WANG Jun1, WANG Yi-fang2,3, HAN Dong-hai2,*

(1. College of Food and Biological Engineering, Xuchang University, Xuchang 461000, China; 2. College of Food Science and Nutritional Engineering, China Agricultural University, Beijing 100083, China; 3. Xuchang Food and Drug Administration, Xuchang 461000, China)

The objective of this study was to develop a method to determine the fat content in yuba by nea r infrared (NIR) spectroscopy combined with chemometrics. A total of 180 yuba samples collected at different occasions from different production lines were tested by NIR spectroscopy. The diffuse ref l ectance spectra (4 000 10 000 cm-1) were collected using an integrating sphere attachment. In order to eliminate the particle scattering and baseline drift, the NIR ref l ectance spectra were preprocessed by 2ndorder derivative with Savitzky-Golay. Backward interval partial least squares (BiPLS), synergy interval partial least squares (SiPLS), searching combination moving window partial least squares (SCMWPLS) and genetic algorithms partial least squares (GAPLS) were employed to extract informative variables and construct quantitative models for the fat content in yuba. After comparison, the best model was obtained by GAPLS method with 143 data points. The correlation coeff i cient (r) was 0.96 and the root mean square error of cross-validation (RMSECV) was 0.95 in calibration set, and the r was 0.92 and the root mean square error of prediction (RMSEP) was 1.17 in prediction set. This work demonstrates that variables extraction methods not only allow selection of the NIR informative variables for the fat content of yuba and simplify the models, but also highlight the potential of NIR technique for assessing the quality of yuba on-line.

near infrared spectroscopy (NIRS); yuba; fat; variables extraction; quantitative detection

O657.3;TS214.2

A

1002-6630(2014)18-0136-05

10.7506/spkx1002-6630-201418027

2013-12-03

國家自然科學基金青年科學基金項目(31401579);河南省科技攻關計劃項目(122102210247)

王加華(1979—),男,副教授,博士,研究方向為食品質量控制與檢測技術。E-mail:w.jiahua@163.com

*通信作者:韓東海(1958—),男,教授,博士,研究方向為食品質量無損檢測技術。E-mail:handh@cau.edu.cn

猜你喜歡
信息模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
訂閱信息
中華手工(2017年2期)2017-06-06 23:00:31
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
展會信息
中外會展(2014年4期)2014-11-27 07:46:46
一個相似模型的應用
信息
建筑創作(2001年3期)2001-08-22 18:48:14
健康信息
祝您健康(1987年3期)1987-12-30 09:52:32
主站蜘蛛池模板: 国产精品专区第1页| 又爽又黄又无遮挡网站| 欧美a在线看| 亚欧成人无码AV在线播放| 亚亚洲乱码一二三四区| 72种姿势欧美久久久久大黄蕉| 久久久久国产精品嫩草影院| 亚洲国产AV无码综合原创| 香蕉国产精品视频| 国产性爱网站| 亚洲丝袜第一页| 国产日韩精品一区在线不卡| 国产第一页屁屁影院| 精品免费在线视频| 看国产毛片| 国产91丝袜在线播放动漫 | 亚洲天堂成人在线观看| 99精品欧美一区| 久久精品人人做人人爽| 爱爱影院18禁免费| 天堂网亚洲系列亚洲系列| 毛片免费高清免费| 成人一级黄色毛片| 99在线视频精品| 亚洲中文字幕97久久精品少妇| 色综合婷婷| 高清无码不卡视频| 国产美女在线观看| 国产91线观看| 精品久久久久无码| 第一页亚洲| 亚洲首页在线观看| 精品无码日韩国产不卡av| 91福利在线观看视频| 日韩AV手机在线观看蜜芽| 久久久久88色偷偷| 国产午夜不卡| 色综合久久88| 国产精品第一区在线观看| 国产成人精品亚洲日本对白优播| 午夜综合网| 亚洲国产天堂在线观看| 国产网站一区二区三区| 永久免费av网站可以直接看的| 国产成人精品午夜视频'| 9久久伊人精品综合| 在线亚洲小视频| 在线精品亚洲一区二区古装| www欧美在线观看| 思思99思思久久最新精品| 国产亚洲精久久久久久久91| 精品福利国产| 暴力调教一区二区三区| 91精品国产91久无码网站| 天天躁狠狠躁| 国产欧美日韩va| 国产欧美成人不卡视频| 国产成人综合网| 国产区免费| 久久这里只有精品23| 五月六月伊人狠狠丁香网| 国产地址二永久伊甸园| 无码免费视频| 在线观看欧美国产| 亚洲91精品视频| 幺女国产一级毛片| 青青青国产视频手机| 久操线在视频在线观看| av在线人妻熟妇| 国产精品性| 亚卅精品无码久久毛片乌克兰| 经典三级久久| 天天做天天爱夜夜爽毛片毛片| 久久无码免费束人妻| 亚洲 欧美 偷自乱 图片| 午夜一级做a爰片久久毛片| 欧洲欧美人成免费全部视频 | 免费人成黄页在线观看国产| 亚洲国产天堂久久九九九| 中文字幕1区2区| 真实国产乱子伦高清| 一级毛片免费观看不卡视频|