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

現實場景中非接觸式心率檢測方法研究

2022-05-15 06:35:48楊學志吳克偉
計算機工程與應用 2022年9期
關鍵詞:信號實驗檢測

牟 睿,陳 鯨,楊學志,吳克偉,方 帥

1.合肥工業大學 計算機與信息學院,合肥230009

2.工業安全與應急技術安徽省重點實驗室,合肥230009

3.合肥工業大學 科研院,合肥230009

隨著社會的發展,人們的生活水平在不斷提高,生活方式也發生了很大變化。人們越來越注重自己的身體健康情況,生命科學已然成為了人類研究的重點。心率是人體最基本也是最重要的生命特征,它反映人體健康水平,更是與人的心臟疾病密切相關,因此心率檢測成為了解自身健康的重要手段[1]。傳統的接觸式檢測法中傳感器設備需要與人體皮膚接觸,無法滿足特定環境和特定人群的心率檢測。為了避免接觸,尋找一種基于視頻的心率檢測技術具有重要意義?;诔上袷焦鈱W體積描記技術(imaging photoplethysmography,IPPG)來進行心率檢測是最多的一類不依賴于硬件的心率檢測方法[2]。IPPG 是一種利用血液中血紅蛋白對不同波段的光線吸收不同,通過測量傳輸或反射回來的光線變換來提取血容量脈沖,進而檢測心率的一類代表性、非接觸和低成本的心率測量方法。但這類方法需要被測試人員處于光照穩定的環境中。若干學者將心沖擊描記技術應用于非接觸式的心率測量中[3],BCG克服了IPPG需要穩定光照的缺點,并且具有非接觸式的特點,僅需成像設備采集一段視頻就能進行心率檢測。

文獻[4]創新性地提出了一種基于BCG原理的非接觸式心率測量方法。該方法從一段記錄頭部運動的視頻中,通過追蹤大量特征點后使用主成分分析來講軌跡分解成一系列的主成分,進而提取心率。然而,由于文獻[4]只使用了Y軸方向上的軌跡成分信號進行分析,因此往往需要跟蹤多達500~1 000個特征點的運動。識別這些特征點,然后對大量的特征點進行跟蹤將會需要很長的求解運行時間,難以滿足實際心率測量的需要[5]。文獻[6]通過框選臉部區域,定位出Shi-Tomas 角點[7]之后,采用KLT[8]跟蹤算法來進行目標點跟蹤,獲取人體頭部的運動軌跡信號,完成了基于單角點的心率檢測。雖然該方法克服了文獻[4]需要跟蹤大量特征點,但是在復雜的光照條件下,會增加角點計算的不可靠性導致算法魯棒性較差。文獻[9]提出了一種基于預決策金字塔層數的視頻微弱運動放大技術,對原始輸入視頻中的微弱運動(包括頭部運動)進行放大,再對放大后的視頻進行處理,從而提取心率信號,但該方法在放大頭部運動的同時也放大了視頻中其他噪聲信號,心率測量精度較低[10]。文獻[11]利用激光作為主動光源,結合附著在頭部的平面鏡與人體身后的平面鏡,實現頭部運動的放大;同時利用加權質心跟蹤算法提取頭部運動軌跡得到BCG信號。該方法需要較為復雜的實驗環境和需要提前進行距離的矯正,不便于實際心率測量的應用。文獻[12]同樣通過追蹤多個特征點,取所有特征點垂直方向的平均位移值作為BCG 信號,最后通過改進的經驗模態分解[13](improved complete ensemble empirical modes decomposition,ICEEMD)對BCG 信號進行處理獲取心率值。雖然該方法在靜態場景中利用ICEEMD獲得的脈搏波信噪比較小,但是在具有較多干擾的場景下,該方法的魯棒性較差。

針對上述方法的缺陷,現提出一種方法能夠克服光照干擾并能快速實現心率檢測SD-BCG,提高現實環境中視頻心率檢測的魯棒性。該方法通過方向可調濾波器過濾包含頭部運動信息的視頻信號獲得頭部在垂直方向上的振動軌跡,計算振動軌跡的相位差分信號,以此弱化光照變化及運動干擾對振動檢測的影響。在此基礎上,通過帶通濾波和主要成分分析等方法從相位差分信號中分離出脈搏信號,并利用傅里葉變換從脈搏信號的頻率中提取心率。

1 本文算法

本文提出一種基于方向可調濾波器的頭部微小振動檢測方法(SD-BCG)。該算法以頭部視頻為基礎,利用設計的空間方向可調濾波器提取豎直方向上的頭部微小振動獲得頭部BCG信號。檢測心率的方法具體框圖如圖1所示。首先對輸入視頻進行人臉檢測,選擇臉頰或者額頭ROI區域,將該區域所有像素點輸入空間方向可調濾波器中得到垂直方向的相位差信號,對其進行時域帶通濾波得到包含心率信息的相位差信號,再使用主成分分析[14](principal component analysis,PCA)對信號進行降維處理,選擇最大頻率和其一次諧波占總功率譜百分比最大的信號,將該信號進行心率計算就能得到心率值。

1.1 特征區域選擇

實驗首先需要獲取人臉部區域的視頻,即選擇感興趣的區域(region of interest,ROI)??梢允謩踊蚓幊虒崿F自動地選取較小的區域,當編程實現ROI 選擇時:實驗使用OpenCV中的ViolaJones臉部檢測器[15]找到包含人臉的矩形區,如圖1 人臉檢測與特征提取提取所示,將眼球之間的距離定義為4d。那么ROI區域為距兩眼連線上方水平距離d,寬度為眼距中間3d,高度為3/2d的紅色矩形區域??紤]到實驗過程中被測人員會有不自覺的眨眼動作,眼睛的微弱運動會影響實驗的魯棒性,因此,實驗把矩形區域的人眼綠色區域去掉,移除的人眼區域占整個矩形區域的20%到50%。

圖1 SD-BCG算法流程圖Fig.1 SD-BCG algorithm workflow

1.2 構造空間方向可調濾波器與相位差信號提取

1.2.1 空間方向可調濾波器

實驗方法需要提取豎直方向上頭部微小振動,根據文獻[16]可知,可以通過讓圖像與可調濾波器卷積能夠獲取任意方向上圖像的細節信息。高斯的二階導數是光滑的并且與邊界非常相似,它可以寫成成一個圓對稱的窗函數和一個多項式的乘積。因此本文使用高斯的二階導數作為基函數,插值函數為高斯的二階導數在90°方向上對應的常數??臻g方向可調濾波器的一般表達式:

式中,q(x,y,θ)表示θ方向的濾波圖像,hi(x,y)表示基濾波器,ωi(θ)表示插值函數。

然后定義三個基濾波器hi(x,y)以高斯二階導數作為基函數,θ=0,π/3,2π/3 作為基濾波器的方向,則有表達式:

根據文獻[17]選插值函數ωi(θ)去掉所有非零項在去除實部和虛部得到等式:

對公式(5)進行計算,得到:

由公式(2)~(4)和(6),對θ取90°,得到能提取豎直方向信息的空間濾波器,濾波器公式為:

將圖像與θ=π/2 方向上的濾波器做卷積運算,可以有效地提取圖像在與θ=π/2 方向上的細節特征,這里就可以提取到每幅圖像豎直方向的相位信息。

1.2.2 相位差信號提取

上文設計的空間方向可調濾波器可以提取圖像豎直方向的相位信息,將視頻每幀卷積濾波器得到每幀的相位信號然后進行算數處理得到相位差信號。假設視頻輸入信號在時刻t,以視頻中ROI 區域內的1 個像素點為例,位置x處的像素,該像素點強度表示為I(x,t),隨著時間的變換,像素豎直方向的位移變化導致像素點的強度發生變化,可以表示成如下公式:

式中,δ(t)表示豎直方向的運動幅度。首先,利用傅里葉級數分解,將位移圖像輪廓f(x-δ(t))表示成多個復雜的正弦曲線之和,具體公式如公式(9)所示:

其中,每一個子帶都對應著1個頻率ω,?ω則代表著相位信號。根據FFT(fast Fourier transform,FFT)變換的平移性質,可以計算出I(x,t)的傅里葉變換結果為:

對于上式中的每一個頻率ω,其復指數形式均如式(11)所示:

在t=0 的時候,相位是?ω,其他時刻的相位為?ω-ωδ(t)。將t≠0 時刻的相位與t=0 時刻的相位做減法,得相位差表示如下:

則點x處的相位差是ROI 區域所有像素點的1 維信號,500 是視頻幀的個數。然后計算ROI 中所有像素點的相位差信號,就完成了頭部相位信號的提取。

1.3 時域濾波

通過方向可調濾波器得到的頭部相位差信號包含多種組成成分,需要進行時域濾波選擇出包含盡可能多心率信息的信號。通常在靜息條件下,一個正常成人心率在[0.75,2]Hz,也就是心跳每分鐘跳[45,120]次[18]。頻率在0.75 Hz 以下的信號主要由低頻運動導致,如呼吸和姿勢變化,需要舍棄這部分信息,頻率在[2,5]Hz的信號在后續實驗過程峰值檢測中選擇中提供一次諧波信號,該部分信號需要保留。考慮保留相位差信號的頻帶范圍,結合巴特沃斯濾波器具有最大平坦幅度的特性,選擇通頻帶為[0.75,5]Hz的5階巴特沃斯濾波器對頭部相位差信號進行濾波,去除不包含心率信息的信號。頭部其中一個橡樹點的相位信號差如圖2(a)所示,經過帶通濾波后相位信號差如圖2(b)所示。

圖2 相位差時域圖Fig.2 Wave in time domain of phase difference

1.4 主成分分析

經時域濾波后的頭部相位差信號包含了由心血脈沖引起的頭部運動以及呼吸、前庭運動和其他面部表情變換等運動信號。必須將這種混合運動分解成子信號來隔離脈沖。將ROI 區域的所有像素點在每個幀上的多維位置作為一個單獨的數據點,并使用PCA 找到一組位置變化的主要維度。選擇一個維度,將數據點序列投影到該維度獲得心率的脈沖信號。

將每一個點的相位差信號作為1個樣本數據,組成1個輸入數據矩陣:

其中,t表示視頻幀數,N表示ROI分塊中的像素個數,N=60×60-20×20×2,共計2 800個像素點。

求輸入矩陣的協方差矩陣如式(14)、式(15)所示:

協方差矩陣的特征向量和特征值如式(16)所示:

式中,Λm為對角矩陣,對應求出的特征為λ1,λ2,…,λN;Φm為對應的特征向量,?1,?2,…,?m組成特征向量矩陣。

將特征向量按特征值的大小進行從上到下的排列,取前T個特征值所對應的特征向量,并組成矩陣QL。這里,T個特征值對應著T個主要成分,其中參數T是通過計算主成分的貢獻率而得到的。T個主成分一一對應的信號都可以通過對特征向量矩陣進行反投影得到[18],如式(17)所示:

視頻中有頭部在其中異常移動,如吞咽、姿勢調整等動作。這種移動增加了位置向量的方差,從而影響PCA 分解。為了處理這個問題,可以在執行PCA 之前丟棄具有最大L2 范數的向量的百分比α。然而,所有最大L2范數仍然必須在投影步驟中使用用來產生完整的信號。實驗表明,當α等于25%時,能夠有效地消除頭部異常移動帶來的影響,實驗效果最佳。因此本實驗α值設定為25%。

1.5 信號選擇

BCG頭部相位差信號通過帶通濾波和PCA主成分分析后得到需要從中提取出心率信息的特征向量。雖然排序靠前的主成分方差較小能夠解釋大部分心率信息,但是方差最小的主成分信號可能無法清晰地表示為心率信號,需要選擇周期性最強的主成分。實驗證明原始相位差信號的98%以上的信息在前5個主成分,故選擇前5個主成分其中一個解釋為心率信號,然后計算主成分最大頻率和其一次諧波頻率相對于功率譜的占比,選擇占比最大的主成分,即為提取的脈搏波波形。所提取的脈搏波圖如圖3所示。

圖3 脈搏波示意圖Fig.3 Schematic diagram of pulse wave

1.6 心率計算

上文提取到代表心率信號的主成分信號。進行呼吸率估計最常用的方法是用快速傅里葉變換(FFT)分析功率譜,最大功率對應的頻率即為心率[19]。估算心率流程如下:

輸入視頻總幀數為N=500 幀,及一維信號長度為500的脈搏信號y和功率譜p(f),計算功率譜中峰值所對應的頻率fmax。

根據公式(18)求得的最大頻率,計算心率值:

式中,HRy為心率值,實現心率檢測。

2 實驗裝置和結果分析

2.1 實驗裝置

實驗裝置如圖4 所示。使用一個普通的網絡攝像頭連接筆記本電腦,攝像頭對準人體頭部,距離被測人員1 m 左右距離。所有被測試人員需要靜止坐在固定椅子上,保持自然狀態。

圖4 實驗場景Fig.4 Experimental scene

電腦通過Matlab 控制攝像頭錄制視頻。設置攝像設備采集圖像為RGB 彩色空間,幀率為30 frame/s,視頻分辨率為1 280×720,總計500 幀圖像。本文對40 位測試者進行了對比實驗,其中包括28 名男性和12 名女性,志愿者的年齡分布于18~60歲之間。在每次測試的時候,采集受測者的頭部視頻信息,同時使受測者佩戴基于指尖的PWS-20D 脈搏波檢測設備,該設備采用光學體積描記術原理,基于指尖進行脈搏波的采集,可以精確地采集受測者的脈搏波信息,在拍攝視頻時,使用該設備測量脈搏波作為參考值。

為了全面的估計本文提出SD-BCG,本文復現了文獻[4]、文獻[7]心率檢測算法作為對比分析,采用了5 個相關參數對其定量描述。

如果測得N組心率值和真實值數據,定義心率值與脈搏波檢測設備之間的誤差為eN=PBCG-Preal,其中PBCG為通過SD-BCG 測得的心率值,Preal為通過脈搏波檢測設備測得的心率值。

平均誤差Me表示對在N組實驗數據,對誤差的絕對值求平均,用來表示其期望值,其計算公式:

標準方差Se表示N組實驗數據的離散程度,其計算公式:

均方根誤差RMSE表示N組實驗數據中檢測值與真實值的偏差,其計算公式為:

平均準確率HRac表示N組實驗數據中檢測值相對于真實值的準確性,其計算公式為:

皮爾森相關系數ρ表示N組實驗數據中檢測值與真實值的線性關系,其計算公式為:

2.2 實驗裝置

2.2.1 抗光照變換能力的心率檢測

由于本文方法是基于頭部運動提取心率,所以在光照變化對實驗結果影響較小。為了充分驗證所提出的方法抗光照干擾的能力,在自然光、日光燈和室內環境3個不同光照度環境下,按照上述實驗步驟,對10位志愿者進行非接觸式心率檢測,其中,每一種光照環境采集5組視頻,視頻長度為20 s,幀率為30 frame/s,共計50 組數據,同時使用PWS-20S 脈搏波檢測儀同步檢測心率信息。結果如表1所示。

表1 不同光照條件下的心率檢測結果Table 1 Heart rate detection results under different light conditions

如表1 可知,在三個不同的光照強度下的Me低于1.6 beat/min,Se和RMSE低于2.3 beat/min,HRac和ρ都在94%以上。數據說明:不同光照強度下,SD-BCG都可以準確地檢測出心率,具有良好的準確性和穩定性。經實驗分析,本算法具有較好的對光照變化的抗干擾能力是因為提取BCG 運動未利用亮度信息,場景變換導致性能略微變化其原因可能是:光照變換可能導致采集的視頻質量較差,讓算法性能有輕微降低。

2.2.2 抗運動干擾能力的心率檢測

上述實驗驗證了本文方法的抗光照干擾能力,下面進一步驗證本文算法的準確性和抗運動干擾的能力。在本次實驗中,每位志愿者可以自由地左右小幅度搖晃頭部或者上下小幅度晃動頭部,按照本文算法的實驗步驟,在自然光條件下,對10名志愿者進行非接觸式心率檢測,其中包括靜止狀態的SD-BCG檢測和頭部輕微搖晃的SD-BCG檢測2組數據,同時使用PWS-20S脈搏波檢測儀同步檢測心率信息,每名受測者包括非接觸式和接觸時心率測量各3次,共60對數據。具體的實驗結果如表2所示,為了進一步直觀展示本文算法的準確性和抗運動干擾的能力,使用Bland-Altman法和一元線性回歸法對實驗結果進行一致性分析,圖5展示了實驗結果的散點圖。

表2 不同運動場景下的心率檢測結果Table 2 Heart rate detection results under different exercise scenarios

由圖5(a)可知,在靜態場景和動態場景,通過SDBCG 檢測的心率數據集中在線性回歸線附近,動態場景比靜態場景略微分散,線性回歸線的斜率接近1,表明SD-BCG 方法的檢測值與心率參考值相關性很強。由圖5(b)的Bland-Altman法對檢測結果進行分析,靜態場景置信區間為[-2.348 6,2.456 1]動態場景的置信區間為[-3.198 4,3.305 5],說明在靜態場景和動態場景,SDBCG方法檢測心率與接觸式脈搏波測量設備具有很好的一致性。

圖5 結果散點圖Fig.5 Result scatter plot

由表2可知,靜態場景和左右晃動場景下,Me小于3 beat/min,Se和RMSE均小于3.3 beat/min,HRac和ρ均大于0.95。實驗數據表明:SD-BCG 的性能受左右運動相較于上下運動影響較小,更能抗左右方向上的運動干擾,能夠在運動干擾的條件下較為準確完成心率檢測。

2.2.3 逐步實驗

上述實驗驗證了SD-BCG 檢測心率的正確性和抗干擾的能力,下一步研究每個步驟的單獨貢獻。本次實驗中,光照環境為自然光,志愿者保持靜止狀態坐在攝像機前,對10位志愿者進行非接觸式心率檢測,實驗步驟按照表所示分別進行,每種實驗步驟進行一次,采集3組視頻,同時使用PWS-20S 脈搏波檢測儀同步檢測心率信息。如表3所示。

表3 可知,實驗步驟會影響本文算法的整體性能,缺少算法的步驟會對心率檢測結果產生嚴重影響。算法步驟中對實驗結果Me的準確率分別上升了3.38%、7.29%和1.52%,Mmax的準確率分別上升了5.05%、7.28%和2.19%,其中對算法性能影響最大的步驟是PCA,影響最小的步驟是SS。經實驗分析:未經PCA 處理的信號包含很多與心率信號無關的信息,帶通濾波并不能過濾這部分信息,經PCA 處理后能夠較為準確的心率信號,SS該步驟主要的作用是完成信號的篩選,對實驗結果影響程度較低。

表3 對所提方法各步驟的檢測結果Table 3 Detection results for each step of proposed method

2.2.4 干擾環境下的對比實驗

抗光照變換能力的心率檢測驗證了SD-BCG 進行心率檢測具有抗光照干擾的能力,抗運動干擾能力的心率檢測驗證了SD-BCG 進行心率檢測具有抗運動干擾的能力,為了進一步驗證SD-BCG

抗運動干擾和抗光照干擾的具體性能。設置了兩組實驗,第一組實驗為志愿者靜止坐在攝像機前,選擇不同的光照環境,第二組實驗為穩定的光照環境,志愿者坐在攝像機前方,允許自由地左右小幅度搖晃頭部。按上述實驗步驟,完成SD-BCG 心率檢測,同時使用PWS-20S 脈搏波檢測儀同步檢測心率信息。為了充分驗證所提方法抗運動干擾和抗光照干擾能力,在上述兩組實驗場景中分別與文獻[4]、文獻[6]和文獻[12]進行對比。其中文獻[4]利用LK 光流法追蹤頭部多個特征點獲取BCG 信號完成心率檢測,文獻[6]利用KLT 追蹤單個取Shi-Tomas 角點后通過KLT 光流追蹤角點Y軸坐標獲得BCG 信號,文獻[12]利用LK 光流追蹤頭部多個特征點后對特征點信號進行ICEEMD 處理獲得心率值。將實驗結果與文獻[4,6,12]比較,對比結果如表4所示。

表4 干擾場景下的心率檢測結果Table 4 Heart rate detection results in disturbed scenes

根據表中數據分析,在有光照干擾或者運動干擾的實驗場景下,本文所提出的方法Me均小于3 beat/min,Se和RMSE均小于3.5 beat/min,HRac和ρ均大于0.95。本文所提出的方法相較于獻[4,6,12]的方法在有光照干擾下,準確率提高分別提高了8.28%、5.97%和4.44%在運動干擾下準確率分別提高了11.86%、5.09%和4.76%。且各方面指標均優于相比于文獻[4,6,12]的方法,具有明顯優勢,主要原因如下:

(1)本文方法是利用方向可調濾波器提取豎直方向上的相位信息,然后進行帶通濾波過濾與心率信息無關的頻率,在通過主成分分析得到噪聲較小的心率信號。因為提取的信號是豎直方向的相位信號所以本文方法受光照影響程度低具有較強的抗光照干擾的能力,同時具有一定的抗運動干擾能力。

(2)文獻[4,6,12]都是利用運動追蹤獲取特征點的運動軌跡從而得到關于心率的BCG 信號,這類方法對光照和運動干擾較為敏感,檢測精度隨之降低。雖然文獻[12]對信號進行了進一步的ICEEMD信號處理,去處了部分同頻噪聲干擾,可抗干擾能力與本文算法相比較弱。

3 結束語

非接觸式BCG心率檢測技術檢測微弱頭部運動聚焦在通過運動追蹤獲取運動信息。但是,由于運動追蹤本身存在的缺點,存在運動干擾、算法復雜度高和易受光照影響等問題。針對以上問題,提出一種基于頭部運動的非接觸式心率測量方法。該算法利用方向可調濾波器提取頭部運動的相位差信息,通過對相位差信息進行濾波和主成分分析后計算出心率值。在多組實驗數據集上進行對比實驗首先進行了抗光照干擾實驗和抗運動干擾實驗驗證了本文算法的抗干擾能力,接著分析了本文算法每步驟的貢獻度,實驗證明主成分分析能夠有效的提高檢測正確率。最后與其他三種主流BCG檢測心率方法對比,有光照干擾下,準確率提高分別提高了8.28%、5.97%和4.44%,在左右方向上的運動干擾下準確率分別提高了11.86%、5.09%和4.76%。實驗結果表明,本文在存在光照干擾和左右方向上的運動干擾下能夠有效提高心率檢測精度。在實驗中發現,當頭部存在豎直方向上的運動干擾時,本文方法效果沒有特別理想,在未來工作中會提高算法魯棒性。

猜你喜歡
信號實驗檢測
記一次有趣的實驗
“不等式”檢測題
“一元一次不等式”檢測題
“一元一次不等式組”檢測題
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
做個怪怪長實驗
基于FPGA的多功能信號發生器的設計
電子制作(2018年11期)2018-08-04 03:25:42
NO與NO2相互轉化實驗的改進
實踐十號上的19項實驗
太空探索(2016年5期)2016-07-12 15:17:55
主站蜘蛛池模板: 麻豆精品在线视频| 丁香六月综合网| 国内精品视频| 少妇露出福利视频| 亚洲国产无码有码| 2021亚洲精品不卡a| 永久成人无码激情视频免费| 国产国产人成免费视频77777| 欧美精品1区2区| 亚洲第一中文字幕| 性色一区| 亚洲成网站| 免费A级毛片无码无遮挡| 亚洲无码视频一区二区三区| 亚洲精品免费网站| 久青草国产高清在线视频| 午夜视频日本| 91网在线| 天堂成人在线视频| 8090成人午夜精品| 国产午夜人做人免费视频中文| 久久窝窝国产精品午夜看片| 国产成人一区免费观看| 欧美中文字幕一区| 国产精品思思热在线| 国产拍揄自揄精品视频网站| 国产99欧美精品久久精品久久| 无遮挡国产高潮视频免费观看| 免费人成在线观看成人片| 香蕉视频在线观看www| 国产第一页屁屁影院| 亚洲美女高潮久久久久久久| 日本欧美午夜| 亚洲清纯自偷自拍另类专区| 成人福利在线观看| 人人爱天天做夜夜爽| 亚洲性日韩精品一区二区| 欧美成人综合在线| 国产精品性| 久久亚洲中文字幕精品一区| 综合网天天| 国产精品亚洲一区二区三区z| 欧美亚洲综合免费精品高清在线观看| 性欧美久久| 波多野结衣久久高清免费| 一级全黄毛片| 午夜精品福利影院| 国产欧美日韩另类精彩视频| 成人在线天堂| 亚洲色图欧美一区| 九九热这里只有国产精品| 久久久噜噜噜| 片在线无码观看| 色婷婷色丁香| 色天堂无毒不卡| 精品国产美女福到在线直播| 亚洲精品成人片在线观看| 日韩乱码免费一区二区三区| 乱人伦视频中文字幕在线| 精品国产电影久久九九| 美女一区二区在线观看| 国产成人夜色91| 91精品免费高清在线| 欧美亚洲中文精品三区| 亚洲国产精品一区二区高清无码久久| 欧美午夜视频| 久热精品免费| 最新国产在线| 久久无码免费束人妻| 极品av一区二区| 中文字幕亚洲另类天堂| 九九久久99精品| 国产综合色在线视频播放线视 | 久久精品人人做人人爽电影蜜月 | 精品91视频| swag国产精品| 国产91全国探花系列在线播放 | 亚洲人成影院在线观看| 欧美激情成人网| 国产幂在线无码精品| 中文字幕免费视频| 最近最新中文字幕在线第一页|