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

頭部旋轉運動下自適應非接觸魯棒性心率檢測方法*

2022-03-18 10:14:32巴圖巴雅爾歐赟趙躍進3孔令琴3董立泉3劉明3惠梅
物理學報 2022年5期
關鍵詞:測量信號

巴圖巴雅爾·歐赟 趙躍進3)? 孔令琴3)? 董立泉3) 劉明3) 惠梅

1) (北京理工大學光電學院,北京 100081)

2) (北京理工大學,精密光電測試儀器及技術北京市重點實驗室,北京 100081)

3) (北京理工大學長三角研究院(嘉興),嘉興 314019)

基于人臉視頻的生理信號檢測面臨的主要挑戰是運動偽影噪聲.針對受試者頭部剛性旋轉運動引起的偽影噪聲,本文提出利用頭部運動信息構建自適應濾波器的非接觸式心率檢測方法.該方法利用人臉二維和三維的特征點計算受試者運動中頭部的偏航和俯仰歐拉角度,并將其作為調控過程噪聲協方差的信號質量指數,進而構建了自適應Kalman 濾波器,實現了穩健的心率估計.實驗結果表明:本文提出的方法可有效抑制頭部剛性旋轉運動引起的噪聲,平均絕對誤差為2.22 beat/min,均方根誤差為2.76 beat/min,與現有方法相比準確度分別提升9%與24.6%,具有統計顯著性.本文提出的頭部旋轉角度自適應非接觸魯棒性心率檢測方法在自發運動的真實場景下能有效提升檢測的準確性,擴大了成像式光電容積描記技術在視頻健康監測領域的使用場景.

1 引言

心腦血管疾病已經成為人類健康的主要威脅,根據心腦血管疾病早發現、早治療,積極干預的防治需求[1],尋求一種操作簡單、結果準確、可重復性好、并可應用于大規模人群日常檢測的測量技術具有重要意義.成像式光電容積描記(imaging photoplethysmography,IPPG)技術是在光電容積描記(photoplethysmography,PPG)技術基礎上發展的一種基于成像設備的非接觸生理參數檢測技術[2,3],克服了傳統PPG 技術存在的缺陷,正逐漸應用于遠程醫療監控等領域,成為近年來研究的熱點[4,5].

基于IPPG 的心率測量技術研究已經取得了重要進展,可以在受試者靜止及輕微非剛性運動狀態下實現心率測量[6,7],運動偽影是目前限制該技術的主要瓶頸.解決運動偽影的常用方法主要包括基于數據的顏色通道線性組合法、基于先驗知識的顏色通道組合法及基于圖像預處理的方法.基于數據的顏色通道線性組合法是依據血紅蛋白的光吸收在整個光譜范圍內的相對強度差異性,利用數據統計線性組合人體皮膚在三色通道的信號實現運動魯棒性的方法,通常綠色通道具有較強的脈搏信號強度[8].2008 年,Hülsbusch[9]構建了雙色通道線性組合法,首次利用該強度的差異性將PPG 信號和運動偽影分成兩個獨立的信號;2010 年,Poh 等[10,11]擴展了這項工作,提出了通過獨立分量分析(independent component analysis,ICA)方法,并將非高斯性作為獨立標準線性組合三通道信號;Lewandowska 等[12]在2011 年提出主成分分析(principal component analysis,PCA)法重新定義了三通道的獨立線性組合,使用盲源分離技術(blind source separation,BSS)分離了PPG 信號及噪聲.基于數據的顏色通道線性組合法認為攜帶脈搏信號的成分是先驗未知的,同時顯示出最強的周期性,但由于運動偽影同樣具有強周期性使其不適用于真實運動場景.由于基于數據的線性組合法需要較長觀測時間,有研究提出了基于先驗知識的顏色通道組合法,根據不同的經驗推理加權組合三色通道信號克服運動偽影.其中具有代表性的研究如下:2013 年Haan 和Jeanne[13]提出的基于皮膚反射物理模型組合三通道信號獲得色度(chrominance,CHROM)信號的方法,首次實現利用IPPG技術測量單人受試者在健身器材上運動時的心率;2017 年,Wang 等[14]提出了顏色失真濾波算法(color-distortion filtering,CDF),選取不同的顏色特征變化作為噪聲和信號的分離依據,從而實現對RGB 信號濾波;2017 年,Wang 等[15,16]提出的基于正交膚色平面(plane orthogonal to skin,POS)的方法,在RGB 空間中尋找一個與皮膚正交的平面以提取脈沖信號.上述研究均是針對信號處理部分,實際上對圖像進行跟蹤配準等預處理也可以實現在一定程度減小運動偽差影響,如圖像配準法[17]、感興趣區域(region of interest,ROI)空間平均法[6]、ROI 跟蹤法[18]、依據時頻分析法的抗運動頻譜峰值跟蹤[19],及適用于長距離移動的基于自適應變焦系統心率檢測法[20]等.

然而,現有方法均只適用于呼吸、頭部輕微移動等非剛性運動場景,當受試者進行旋轉頭部和行走等自發性運動時無法準確估計心率.針對上述問題,本文提出了一種利用角度信息構造自適應濾波器抑制頭部剛性旋轉運動偽差的方法,提供單一穩健的心率估計.該方法在目前基于IPPG 心率測量研究的基礎上擴大了在自發運動狀況下的使用場景,滿足了在頭部剛性旋轉運動場景下非接觸的實時準確測量心率,從而有效監控運動中受試者的健康狀態,進一步推動了IPPG 在生理信號檢測和視頻健康監測領域的發展.

基于人臉視頻[21]的IPPG 心率測量技術,利用成像設備對包含被測部位的信息進行視頻采集[22],將脈搏信號即由血液容積變化引起的光強變化用視頻圖像的方式記錄下來,再通過對視頻圖像處理提取出脈搏波信號[23-25],最后通過脈搏波特征量的分析實現心率信號的提取[4,26,27].事實上,在普通的生活場景中,如在健身、長期遠程連接醫療監控技術時,難以避免受試者相對于相機的運動.下面將針對受試部位(臉部)剛性旋轉運動引起的運動偽差,分析其產生原因,并對本文提出的噪聲抑制方法進行闡述.

2 方 法

基于人臉視頻的IPPG 心率測量技術,利用成像設備對包含被測部位的信息進行視頻采集,將脈搏信號即由血液容積變化引起的光強變化用視頻圖像的方式記錄下來,再通過對視頻圖像處理提取出脈搏波信號,最后通過脈搏波特征量的分析實現心率信號的提取.

2.1 剛性運動偽差

在運動的場景中,將人臉視作一個剛體,即其在三維空間中的旋轉運動可以分為3 種類型:俯仰(pitch)、偏航(yaw)、翻轉(roll),如圖1 所示.

圖1 人臉的3 個旋轉運動自由度Fig.1.Three rotational freedom degrees of human head.

當生物皮膚組織在各個方向具有均勻照度的理想環境(積分球)中,且觀察角度小于60°時,可近似視作擴展的朗伯輻射源[28].假設人臉在理想光照環境中運動,即將人臉視作擴展朗伯輻射源,其輻亮度與方向無關.

當人臉保持靜止時,成像系統像平面的輻照度E可表示為

其中,τ表示成像系統透過率;L表示ROI 輻亮度;D/f′表示光學系統數值孔徑.

CCD 相機接收到的總功率P可表示為

其中,AEP表示入射光瞳面積.

原始IPPG 信號正比于CCD 接受總功率與ROI 的像面積之比

人臉作為擴展朗伯輻射源,所以當人臉以角度α旋轉時輻亮度L不隨角度α變化,只有ROI 的像面積隨角度α變化,因此(3)式應改寫為

由于當人臉做翻轉運動時CCD 相機獲取的像面大小實際上沒有改變,并且在真實場景中頭部的上下、左右方向的旋轉運動更頻繁,故在本文中僅考慮俯仰和偏航運動對心率信號的影響.

人臉做俯仰和偏航運動會引起 cos2α的變化,從而對IPPG 信號引入誤差.然而真實場景中的環境與均勻照度的積分球的差異較大,人臉ROI 作為接收面,入射光強會隨著接收面法線和輻射方向的夾角改變而變化,導致CCD 相機采集視頻時人臉ROI 作為輻射面的輻亮度L也會隨著之改變.

綜上所述,當人臉做俯仰或偏航運動時,原始IPPG 信號會受旋轉角度α變化及ROI 入射光強變化的調制,繼而引入了運動偽影.

2.2 旋轉運動自適應濾波器

依據貝葉斯算法,在測量心率前對真值的預測為該時刻的先驗概率,在得到測量值后對真實值的估計為后驗概率.由于后驗概率中包含測量過程中產生的噪聲,例如運動偽影等,故利用先驗概率修正信號實現對錯誤心率估計的濾波.本文僅考慮測量過程中旋轉頭部運動的影響,如無特別說明下文中提到的測量噪聲均指旋轉頭部運動引入的偽影.

依據2.1 節的分析,本文利用貝葉斯框架,設計了一種頭部旋轉運動自適應濾波器,其實質是根據人臉旋轉的角度信息自適應地改變過程噪聲協方差估計,在新數據到達時準確估計心率動態變化.

2.2.1 自適應過程噪聲協方差

為從人臉視頻中獲取頭部旋轉角度信息,即頭部姿態歐拉角,首先利用 CLNF 算法定位面部特征點坐標,并將二維點反向投影到三維人臉模型,得到三維點坐標,再通過標定獲得相機標定參數,最后利用最小二乘法求解透視位姿得出頭部姿態歐拉角[29-32].本文僅考慮俯仰角度(pitch)和偏航角度(yaw):

對所獲角度的絕對值進行歸一化并取反操作,使數值范圍在0—1 之間.并滿足旋轉角度較大時數值趨近于0,角度較小時數值趨近于1:

令信號質量指數為

利用信號質量指數θSQI調節測量噪聲協方差估計,構成頭部旋轉自適應協方差RSQI:

其中,R0取經驗值0.1[33].由此將獲得的自適應測量噪聲協方差與Kalman 濾波器結合,使得濾波器在頭部旋轉歐拉角較大即數據置信度較低時,測量噪聲協方差較大.

2.2.2 基于RSQI調控Kalman 增益構造自適應濾波器

本文提出的貝葉斯決策是基于張玉燕等[34]和Nemati 等[35]提出的Kalman 濾波器基礎上實現的.Kalman 濾波器提供了一種有效的遞歸解決方案,它通過組合之前狀態值和當前觀測值,推斷或重構最符合系統狀態真實值的最優估計.

假設第k個的心率真實狀態值為xk,觀測值為yk,分別滿足如下的離散時間差分方程:

其中,A為狀態轉移矩陣;xk-1為前一時刻的心率狀態值;uk為該時刻系統的控制輸入變量;H為系統的測量參數;wk-1為過程激勵噪聲,vk為測量噪聲,二者均為高斯噪聲且相互獨立,他們的正態分布均值為0,協方差分別為.由于人體在穩定狀態下心率不會突變,因此假定A為單位矩陣;忽略控制輸入變量、將測量參數H設為1;忽略過程激勵噪聲的影響,將Q設置為經驗值1.(10)式和(11)式改寫為

假設第k-1個心率為均值為、協方差為的正態分布,即第k-1 個時刻的后驗估計的協方差為

預測步驟

根據貝葉斯規則,第k個時刻的先驗估計及協方差為

即預測步驟為:當前時刻的先驗概率由上一時刻的后驗概率決定.

更新步驟

可以由預測步驟得到第k個時刻的狀態更新,即后驗估計及協方差為

(16)式和(17)式中,Kalman 增益是Kalman 濾波過程中的核心步驟.利用(9)式得到的自適應測量噪聲協方差來調節Kalman 增益系數,定義如下:

當出現偽影時,估計測量噪聲協方差RSQI會增大,顯著降低Kalman 增益,此時降低測量數據的置信度,更信任該時刻的先驗概率分布.

圖2 為本文提出自適應濾波器的算法流程圖.首先,使用頭部姿勢歐拉角信息構建心率信號質量指數,用于估計自適應測量噪聲協方差RSQI.然后,利用RSQI計算Kalman 增益系數;根據頭部姿勢角度值的動態變化,結合Kalman 增益系數、心率測量值和當前時刻的先驗概率,實現對新到達數據的自適應濾波.最終,將當前狀態的后驗概率更新為后續狀態的先驗概率.當某一時刻后不再有新的心率測量值傳入時,退出該更新迭代循環.

圖2 頭部旋轉角度信息作為自適應濾波器質量控制的流程圖Fig.2.Flow chart of head rotation angles as a quality control of our proposed adaptive filter.

3 實驗系統

3.1 實驗方案

實驗方案如圖3 所示.實驗中采集的視頻幀率為30 frame/s,分辨率為1024×768,并保存為AVI格式.所有測量均采用環境光穩定的房間,避免了環境光強突變對實驗結果的干擾.在所有實驗中,CCD 相機采集人臉視頻,然后通過數據線傳輸到電腦端進行存儲和分析.同時,利用指夾式血氧儀(德州儀器公司,AFE4400SPO2EVM)進行心率檢測,并將其測量結果作為參考真值.

圖3 實驗方案圖 (a) 實驗裝置示意圖;(b) 采集到的視頻片段((i)頭部姿勢的歐拉角為0°;(ii) 頭部俯仰運動;(iii) 頭部偏航)Fig.3.Experimental plan diagram:(a) Experiment set-up;(b) video clips ((i) the Euler angle of the head pose is 0°;(ii) the head pitch segment;(iii) the head yaw segment).

參加實驗的志愿者為21 名年齡在23—28 周歲之間的在校研究生,志愿者在參加實驗前被告知了實驗目的且均是自愿參加實驗,所有參加實驗的志愿者均身體健康無心臟病等疾病.實驗中,要求受試者距離CCD 相機0.5 m 保持靜止,并正向面對相機(視線與相機法線平行)隨機做出頭部俯仰及偏航的動作.

3.2 感興趣區域提取

本文設計了如下由面部關鍵點構建的動態臉頰與鼻翼ROI 區域提取方法:選取面部68 個關鍵點中第1,17,42 和47 個關鍵點構建人臉ROI 區域,如圖4(a)所示.由p1和p2關鍵點決定的矩形區域作為ROI 區域,如圖4(b)所示.

圖4 中,p1和p2關鍵點的坐標由下式確定:

式中pxn和pyn(n=1,16,41,46) 代表 點n的x坐 標和y坐標;加(減)50 代表向右(左)偏移50 個像素值,如圖4(b)中藍色圓點所示;h代表當前幀的人臉檢測高度.采用兩側關鍵點分別向內偏移50 個像素值決定矩形框左右邊界,保證ROI 區域內不會隨著頭部的旋轉運動(如圖4(c))而采集到環境背景,從而很好地避免了引入唇色噪聲和背景噪聲.

圖4 矩形ROI 區域選取示意圖 (a)面部68 個關鍵點檢測;(b)臉頰矩形ROI 區域;(c)偏航運動時ROI 區域Fig.4.Flow chart of rectangular ROI area selection:(a) 68 feature points detection;(b) rectangular ROI area;(c) ROI area during yaw motion.

3.3 數據處理

本文提出的心率檢測方法算法流程如圖5 所示.首先采用開源工具OpenFace2.0 對采集到的視頻進行人臉檢測,基于特定人臉關鍵點選取毛細血管分布豐富的臉頰和鼻翼區域作為ROI,對每一幀的人臉圖像提取ROI 區域可以有效克服因受試者頭部移動等運動引入的噪聲偽影.然后將每一幀ROI 的像素均值按視頻時間序列連接生成原始IPPG 信號.然后將每一幀ROI 的像素均值按視頻時間序列連接生成原始IPPG 信號.再采用先驗平滑濾波(prior smoothing filter,PSF)[36]對IPPG信號進行去趨勢處理,去除信號的低頻漂移趨勢噪聲.通過小波濾波算法消除高頻噪聲,具體為首先選擇db12 小波基對原始信號進行6 層分解,由于人體脈搏信號頻率一般在0.7—4.0 Hz 之間,因此選擇第4,5,6 層的細節信號重構脈搏波信號(對應的頻率范圍為0.469—3.75 Hz).最后對去噪后的信號建立10 s 的滑動窗口進行快速傅里葉變換[37],設置滑動步長為0.2 s,并對窗口內的有效數據進行補零操作避免頻譜泄漏,獲得其對應功率譜,并利用功率譜峰峰值所在頻率估算心率.并對指夾式血氧儀獲取的一維時域脈搏信號同樣以10 s 的窗口大小及0.2 s 的滑動步長進行快速傅里葉變換,通過其功率譜計算得到參考真實心率值.

圖5 本文提出方法的流程圖Fig.5.Flow chart of method proposed in this paper.

在上述基礎上,結合人臉二維及三維特征點計算頭部姿勢歐拉角度,進而根據角度絕對值獲得信號質量指數θSQI,利用θSQI估計測量噪聲協方差RSQI從而調控Kalman 增益,以此構建一個頭部旋轉運動的自適應濾波器.該濾波器能夠依據頭部旋轉角度的變化動態濾除旋轉運動引入的噪聲,進而實現單一穩健的心率檢測.

4 結果與討論

為驗證本文提出算法的可行性,將所提出的濾波器與目前性能較佳的CHROM 方法[13]結合并進行對比實驗.

圖6 為某一受試者的測試結果.圖中淺灰色和深灰色虛線分別為頭部偏航和俯仰角度變化曲線,藍線為參考設備指夾式血氧儀測得的心率值.綠線為僅采用CHROM 算法測得的心率值,紅線為加入本文提出的自適應濾波器算法后測得的心率值.從圖6 可以看出,當局部出現運動偽影時,CHROM不能有效抑制運動偽差,而結合本文提出的算法后,能穩健地抑制運動偽差,實現心率準確測量.

圖6 本文提出的濾波器與CHROM 結合的對比結果Fig.6.Comparison results of the filter proposed in this paper combined with CHROM.

為定量描述上述兩種方法的21 組實驗結果與真實心率值之間的一致性,用相關性圖進行了結果分析,結果如圖7 所示.

圖7 給出了CHROM 方法和本文提出的濾波器與CHROM 結合的21 組實驗結果相關性圖.圖中任意一點的橫軸表示實驗的心率估計值,縱軸表示指夾式血氧儀獲得的心率參考真值.從圖7 的相關性結果中可以看出,與傳統的CHROM 方法相比,利用本文所提方法獲得的心率估計值與心率真實值更趨于一致,表明頭部旋轉自適應濾波器顯著提高了心率估計結果與真實心率結果的一致性.

為了進一步證實這一結論說法,將本文設計的濾波器與不同的傳統方法(GREEN,CDF,POS,CHROM)結合前后的性能(誤差)進行定量分析,結果如表1 所列.

表1 為本文設計的濾波器與不同的傳統方法(GREEN,CDF,POS,CHROM)結合前后測得心率值的MAE 和RMSE 誤差結果,從表1 可以看出,4 種方法增加了本文的自適應濾波器后MAE及RMSE 均有明顯提升效果;其中在CHROM 的基礎上組合本文的自適應濾波器后的效果最優,其MEA 及RMSE 值分別為2.22 beat/min 和2.74 beat/min,較CHROM 的結果提升了9%和25%,R2相關系數提升了3%.尤其應注意,CHROM疊加CDF 再與本文的自適應濾波器結合后,R2相關系數最優,高達0.8901,說明CDF 方法與本文的濾波器結合能有效提升心率信號與真實值的相關性.上述結果分析表明,本文提出的頭部剛性旋轉自適應濾波器能有效抑制旋轉運動引入的偽影并且顯著提升估計心率值的準確性.

圖721 組實驗結果相關性圖(圖中紅線表示y=x 的線性關系) (a) CHROM 方法的相關性圖;(b)本文提出的濾波器與CHROM 結合的相關性圖Fig.7.Correlation plots of 21 groups of experimental results (the red lines in the plots indicate linear relationship of y=x):(a) Correlation plot of CHROM method;(b) correlation plot of combining CHROM with our proposed adaptive filter.

表1 本文提出的自適應濾波器與不同的傳統方法結合前后實驗結果對比Table 1.Comparison of experimental results before and after combining the filter proposed in this paper with different traditional methods.

需要說明的是,本文在處理頭部姿態歐拉角度時應用的歸一化操作使用了視頻片段中的最大歐拉角度用于確定信號質量指數,在實時監控受試者的心率值時,會有一定的處理時間延遲,因此我們會在未來的工作中繼續進行相關的實驗研究.Kalman 濾波器參數的初始化值是影響該技術性能的重要因素,本文并沒有確定心率檢測的最優濾波器參數,相關工作將在日后的研究中進行探討.此外,本文提出的由頭部姿態信息估計自適應測量噪聲協方差,并以此構建自適應Kalman 濾波器抑制運動偽差的方法,適用于任何由運動偽影引入噪聲的時間序列信號,可以擴展到其他非接觸式生理參數檢測及跟蹤技術領域.

5 結論

本文針對非接觸式心率檢測時受試者剛性運動引入偽影的問題,提出了一種根據頭部姿態角度作為主要失真指標提高心率測量準確率的方法.該方法利用人臉檢測技術獲得受試者頭部姿態歐拉角度并以此計算信號質量指數θSQI,根據θSQI估計測量噪聲協方差來調控Kalman 增益,從而構建一個頭部剛性旋轉運動的自適應濾波器,利用該濾波器有效提高了心率檢測準確性.與傳統的適應于運動場景的心率檢測算法的對比分析表明,本文提出并設計的自適應濾波器有效解決了頭部剛性旋轉運動引起偽差的問題.該方法使基于人臉視頻的心率檢測技術在現實場景中有效監控受試者自發運動時的健康狀態成為可能,進一步推動了IPPG 技術在生理信號檢測和視頻健康監測領域的發展.

猜你喜歡
測量信號
信號
鴨綠江(2021年35期)2021-04-19 12:24:18
完形填空二則
把握四個“三” 測量變簡單
滑動摩擦力的測量和計算
孩子停止長個的信號
滑動摩擦力的測量與計算
測量的樂趣
測量
基于LabVIEW的力加載信號采集與PID控制
一種基于極大似然估計的信號盲抽取算法
主站蜘蛛池模板: 国产日韩欧美在线视频免费观看| 日本免费新一区视频| 国产亚洲视频播放9000| 露脸一二三区国语对白| 91色综合综合热五月激情| 理论片一区| 四虎影视库国产精品一区| 欧美丝袜高跟鞋一区二区| 欧美激情二区三区| 欧美在线导航| 精品人妻一区无码视频| 精品国产免费观看| 九九免费观看全部免费视频| 在线观看国产精品第一区免费| 国产自产视频一区二区三区| 国产视频自拍一区| 亚洲精品图区| 国产aⅴ无码专区亚洲av综合网 | 国产在线97| 九色最新网址| 国产精品浪潮Av| 亚洲成人精品久久| 国产经典免费播放视频| 国产另类视频| 欧美激情综合| 国产精品高清国产三级囯产AV| 精品福利国产| 亚瑟天堂久久一区二区影院| 人妻少妇乱子伦精品无码专区毛片| 亚洲中文字幕久久精品无码一区| 久久综合九色综合97网| 国产精品第页| 欧美国产在线看| 91亚洲精品第一| 日本草草视频在线观看| 大学生久久香蕉国产线观看| 免费国产在线精品一区| 精品国产成人av免费| 野花国产精品入口| 91亚洲精品国产自在现线| 97国产一区二区精品久久呦| 国产免费一级精品视频| 国产91成人| 国产精品亚洲一区二区三区在线观看| 亚洲妓女综合网995久久| 国产成人无码综合亚洲日韩不卡| 国产亚洲精品自在线| 99视频国产精品| 国产亚洲欧美在线人成aaaa| 久久这里只有精品国产99| 国产一在线观看| 精品久久香蕉国产线看观看gif| 亚洲v日韩v欧美在线观看| 亚洲一级毛片在线观播放| 欧洲在线免费视频| 天天做天天爱夜夜爽毛片毛片| 免费网站成人亚洲| 日韩成人在线视频| 国产欧美精品午夜在线播放| 亚洲国产成人无码AV在线影院L| 91啦中文字幕| 久草视频一区| 波多野结衣中文字幕一区二区| 全部免费特黄特色大片视频| 玩两个丰满老熟女久久网| 国产伦片中文免费观看| 国产亚洲精久久久久久久91| 91年精品国产福利线观看久久| 日韩美毛片| 久久精品这里只有精99品| 性69交片免费看| 国产精品专区第一页在线观看| 亚洲专区一区二区在线观看| 中文毛片无遮挡播放免费| 免费看av在线网站网址| 亚洲欧美日韩另类在线一| 91美女视频在线| 欧美成人免费一区在线播放| 亚洲大尺码专区影院| 久久semm亚洲国产| 大学生久久香蕉国产线观看| 亚洲无码高清免费视频亚洲|