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

嵌入式電容層析成像系統的成像加速方法研究

2024-04-25 00:00:00陸程程胡紅利唐凱豪單宇杰黃敬軒董海健
西安交通大學學報 2024年3期

摘要:為了提高嵌入式系統設計中電容層析成像(ECT)的圖像重構速度,研究了一種針對進階精簡指令集機器加上現場可編程門陣列(ARM+FPGA)硬件架構的圖像重構算法加速技術。針對廣泛應用且魯棒的Landweber迭代算法(ILA),首先分析算法結構,然后基于FPGA的流水線特點,改進ILA涉及的循環結構,從而達到加速的效果。同時,針對ARM+FPGA架構的特點,討論了ARM核與FPGA核各自的任務分配方式,進一步優化了算法速度。為了驗證算法的有效性,分別在使用MATLAB編程和使用提出的加速方法搭建的ZYNQ平臺進行了圖像重構實驗,從圖像重構耗時、圖像相對誤差和圖像相關系數3個指標論證提出方法的有效性。實驗結果顯示,使用搭建的ZYNQ平臺進行Landweber算法成像時,每個圖像的運行時間比使用MATLAB編程的運行時間減少了30%~40%。該研究在保持重構精度的同時有效提升了迭代算法的速度,對于ECT系統的硬件加速具有一定適用性。

關鍵詞:電容層析成像;圖像重構算法;硬件加速;嵌入式系統

中圖分類號:TN911.7 文獻標志碼:A

DOI:10.7652/xjtuxb202403017 文章編號:0253-987X(2024)03-0183-10

Research on an Imaging Acceleration Method for Embedded Systems of Electrical Capacitance Tomography

Abstract:To increase the image reconstruction speed of electrical capacitance tomography (ECT) in embedded system design, this paper proposes an acceleration method for image reconstruction algorithm in the context of the ARM+FPGA hardware architecture. The structure of the widely-used and robust iterative Landweber algorithm (ILA) is analyzed and then the loop structure within the ILA is modified taking advantage of the FPGA’s pipeline features to achieve acceleration. Furthermore, according to the characteristics of the ARM+FPGA architecture, the task allocation of ARM and FPGA cores is discussed to further optimize the speed of the algorithm. To validate the effectiveness of the proposed algorithm, image reconstruction experiments are conducted on both a desktop computer (using MATLAB programming) and the ZYNQ platform (using the acceleration method proposed in this paper). The effectiveness of the proposed method is demonstrated using three parameters: image reconstruction time, relative error of the reconstructed image, and image correlation coefficient. The experimental results show that when using the platform built in this paper for imaging with ILA, the runtime for each image is reduced by 30%—40% compared with the time when using MATLAB programming on a PC. This study effectively increases the speed of iterative algorithm while maintaining reconstruction accuracy, demonstrating its applicability to hardware acceleration of ECT systems.

Keywords:electrical capacitance tomography; image reconstruction algorithm; hardware acceleration; embedded system

為了滿足經濟社會發展對清潔能源的需求、推動天然氣對傳統高碳化石能源進行替代,天然氣行業得到了大力發展。在天然氣運輸過程中,氣-水兩相流是管道內常見的流動狀態,監測并推斷出管道內兩相流流動特性不僅能保障輸送過程中的安全,還能隨時根據監測結果進行調控。因此,長久以來,工業生產中的氣-水兩相流檢測問題一直都是國內外專家學者的研究重點與難點[1-5]。

傳統的電容層析成像(ECT)傳感器由多個敏感電極組成,均勻貼附在被測管道上,通常電極數量為8、12或16個[6],在敏感電極外部,采用屏蔽罩以抑制外界對傳感器的電磁干擾。傳感器工作時,敏感電極陣列接收激勵信號,并由開關陣列電路控制成對工作,直至所有電極完成工作。如此,系統采集到的每一組電信號均包含了極板間電容的所有信息。隨后,通過信號調制系統對采集到的電信號進行調制,以獲得電極板間的電容。調制后的信號傳遞到信號處理系統進行圖像重構,生成介質分布圖像數據。最后,數據傳輸至顯示系統中并顯示。

這一重構過程又被稱為ECT技術的逆問題,領域內對此進行過大量研究,其中使用的算法稱為圖像重構算法,目前的算法主要有迭代法與非迭代法兩種。非迭代法即單步法,將圖像重構近似為線性過程,具有運行速度快、易于實現、計算量小等優點。由于介電常數分布和電容本質上是非線性關系,幾乎不可能通過任何簡化的線性模型獲得準確的解決方案[6]。

應用于ECT的圖像重構迭代算法包含代數重構法(ART)[7]、同步迭代重構法(SIRT)[8]、Landweber迭代法[9]、共軛梯度法(CG)[10]等。目前,由于Landweber迭代法(ILA)相對于其他方法兼顧了時間成本與圖像質量,所以基于Landweber的改進算法[11-12]一直是研究的熱點問題并被廣泛應用。

[HJ1.8mm]在傳統的ECT系統中,圖像重構算法一般通過計算機在通用處理器上運行,但這種方式存在性能限制。與通用處理器相比,可編程門陣列(FPGA)具有更高的計算效率和更低的功耗[14-17]。此外,FPGA還支持深度可變的流水線結構,具有大量的并行計算資源,在一個時鐘周期之內就可以完成多種復雜的功能,具有很大的運算優勢[18-21],并增加了傳感器的靈活性以及實時性。

為提升基于嵌入式硬件平臺的ECT系統的成像速度,本文在FPGA與ARM架構的硬件平臺上搭載了Landweber圖像重構算法,并研究了算法加速策略,提高了嵌入式ECT系統的實時成像性能。

1 電容層析成像的原理與算法

1.1 電容層析成像原理

電容層析成像是一種電學層析成像技術,產生于20世紀90年代。ECT傳感器的結構及電極位置對測量結果精確性有很大影響,結合實際工業的嵌入式優化需要,本文研究的傳感器模型結構為圓形12電極配置[6],如圖1所示。

圖1中ECT傳感器的數學模型可以描述為時諧形式、帶有Dirichlet邊界條件的方程

目前,ECT的主要原理是使用敏感電極板測量電容,再利用圖像重構算法顯示敏感域內與電容相關的介電常數的分布情況,最終得到不同介電常數的介質分布圖像。依據ECT的靈敏度理論,其測量值和介電常數離散分布存在一種線性關系

z=Sg(2)

式中:z為歸一化的電容向量;S為系統的靈敏度矩陣;g為圖像的灰度向量。式(2)描述了ECT的測量過程。靈敏度矩陣S即可通過電磁場數值計算的方法求解,電容向量z可以通過傳感器電路測量,于是ECT的成像的原理即已知S和z,根據式(2)逆推灰度向量g,這一過程被稱為ECT的逆問題。

圖1結構的成像區域被均勻剖分成1961個矩形單元,因此靈敏度矩陣S的大小為66×1961,待求的灰度向量g的大小為1961×1。

1.2 Landweber迭代重構算法原理

本文擬對ECT領域內廣泛使用的Landweber迭代算法(ILA)進行改進,以實現嵌入式平臺上的圖像重構[9]。ILA是一種梯度下降算法,其思想是用迭代的方式逼近式(2)中S的逆矩陣S-1。ILA常用于求解ECT逆問題中L2范數優化模型,基于此模型,由式(2)逆推g的問題可以描述為

令式(3)的目標函數為f,那么

gf=ST(Sg-z) (5)

于是,用梯度下降法求解式(3)的迭代格式可以表達為

gk+1=gk-αgf(gk)=gk+αST(z-Sgk) (6)

式(6)即為ILA的迭代格式。其中,α為松弛因子,通常取為

式中:σ(S)為靈敏度矩陣S的奇異值。按式(7)的方式確定α可以保證迭代格式(6)最終收斂[22]。

2 嵌入式ECT系統的硬件設計

本文設計的硬件系統面向油氣開采規劃和生產等過程的應用場景[16]。在該場景中,為了實現對對象的可靠監測,結合實驗室采集的具體需求,管道內流速約為0.9m/s,因此本文采用的高速攝像機拍攝速度為240幀/s,分辨率為1080像素,系統中算力可以滿足每幀拍攝耗時小于等于0.005s的工業需求。嵌入式ECT主要包含傳感器、測量單元、控制與資料解釋單元。其中,傳感器為團隊自主設計制造的12電極ECT傳感器[22];測量單元為基于時諧電流法的微小阻抗測量電路[23];控制與資料解釋單元是基于ZYNQ平臺實現的嵌入式上位機[16],用于控制ECT的輪巡測量、與測量單元通信和圖像重構。

為了保證ECT傳感器的幾何參數具有較高的均一性,本文采用柔性電極表貼工藝制造了12電極ECT傳感器[24]:管壁材料為聚醚醚銅,采用3D打印工藝制造;電極是基于柔性印制電路(FPC)技術制造的。將FPC電極表貼在聚醚醚銅管壁上,并焊接同軸線纜,最終形成傳感器。傳感器的核心部件與傳感器實物照片分布如圖2和圖3所示。

為滿足高速輪巡測量的要求,本文設計了多路復用開關陣列電路作為傳感器的接口。具體采用模擬多路復用選擇器ADG1606,開關時間為ns級,導通電阻為4.5Ω,其平坦度為1.1Ω,可以滿足數據采集的低失真要求。受控于控制與資料解釋單元,該接口可配置每個電極處于激勵狀態或檢測狀態。接口電路的實物如圖4所示。

為實現可靠的微小電容測量,本文設計了基于時諧電流法原理[23]的測量單元,原理如圖5所示,其中Cx為待測電容,Cs1、Cs2為引線對地的雜散電容。圖5中,采用MAX038芯片設計了激勵信號源,以產生高頻精密的正弦信號,輸出頻率范圍為0.1Hz~20MHz。運算放大器型號為AD8066,具有145MHz工作帶寬。

圖5中Cf、Rf以及運算放大器將經由待測電容流過的電流轉換為電壓輸出Vo1,設激勵信號為Vi=Asin(ωt+θ),根據電路理論,該電路的電容-電壓轉換輸出為

3 FPGA+ARM硬件的加速策略

為進一步保證足夠短的圖像重構耗時,本文研究了針對硬件結構和ILA算法特點的硬件加速方法。圖像重構算法加速本質是對大數據量矩陣運算的加速,在C/C++語言中,矩陣運算是通過for語句的循環來完成的,因此對相關算法的加速關鍵在于對嵌套循環的優化。本文采取了展開與流水線優化組合的優化方法。流水線優化是在硬件資源充足時,在一次循環還未完成時就開始下一次的循環,使得吞吐量得到提升;展開則是把循環復制成幾個電路使其同時運行,兩種方法結合可以進一步提升算法加速的效率。

如圖6所示,在設置了流水線之后,每個循環周期消耗的時鐘周期并沒有發生改變,但是在前一次循環執行結束之前下一次循環已經開始,增加了程序的并行性。

3.1 考慮依賴性的流水線加速

如果只將ILA圖像重構算法中的一部分進行移植并使用硬件加速,就需要保證這一部分的程序所需的交互數據量少、加速前運行時間較長,這樣才能夠保證獲得足夠的程序運行在時間方面的優化。如前文所述,ILA可以分成兩部分嵌套循環,其偽代碼如表1所示。策略1的硬件加速部分偽代碼如表2所示。

實驗中注意到在循環loop1_2中對于浮點數的乘法和累加的運算分別需要3、4個時鐘周期,而循環的流水線優化也被這一過程影響,無法實現并行化。實際上,造成這樣一種現象的原因是由于浮點累加計算的程序中存在一種“反依賴關系”,具體解釋為對變量進行讀操作之后馬上對其進行寫操作。

上述將部分算法進行硬件加速運算的策略可以將特定的算法最大程度的并行化,進而獲得快速的運算。但是,這種加速策略需要處理器系統(PS)端與可編程邏輯(PL)端聯合參與圖像重構的運算,這就使得硬件加速的效果受制于數據交互的傳輸速率以及PS端程序的執行時間。因此,本文將ILA整體全部遷移至PL端并進行加速優化,本文這種方法稱為策略2。

一般情況下,需要加以關注的依賴關系可以分為循環內部的依賴關系和循環不同迭代間的依賴關系。針對于這種情況有一個有效的解決方案:使用一種分組的方式消除變量之間的反依賴關系,將同一個組內所用的寄存器通過存儲器分割的方式完全展開,使得浮點數運算可以并行運行,在循環完畢后將所有組的元素累加,最終得到運算結果。本文稱這種改進寄生關系的優化策略為策略3。策略3的偽代碼如表3所示。

這種方式進一步結合了數組優化與循環優化方式,在總的運算量不變、流水線優化延遲時間不變的情況下,單次循環所需運行周期的增加與計算量增加意味著并行度的進一步增加。因此,總的延遲時間減少。

3.2 結果分析

這些策略的流程優化首先在Vivado軟件上進行測試,實際優化用時將在實驗中實際測量。各優化策略所需的運算時間即時間復雜度如表4所示。

表4中,時間是指系統得到輸入后與計算完畢后得到輸出結果之間的時間差。結果表明,改進寄生關系的策略3優化效果在3種策略中所需的運算時間大幅下降。

4 實 驗

基于ZYNQ平臺的ECT系統原理和實驗平臺如圖7和圖8所示。

在傳統的電學層析成像系統中,信息處理部分一般采用計算機通用處理器進行運算。使用通用處理器開發周期短,但它具有更大的體積、更高的功耗和成本。如果能夠加入設計好的專用硬件則可以使用并行架構加速運算,減少擁有大計算量的圖像重構運算的運算時間,提高系統的吞吐量,使系統能夠在高動態要求的情況下得以應用。現場可編程門陣列(FPGA)具有極好的邏輯容量以及高度的使用靈活特性,是硬件并行計算器的理想平臺。

在ZYNQ中,硬件部分在ECT系統中主要作用包含數據傳輸以及圖像重構硬件加速算法的處理。使用Vivado來進行編譯、綜合以及實現,該系統的整體框架如圖9所示。

ECT系統采集到的數據通過通用異步收發器(UART)串口傳輸到ZYNQ系統中,而后使用其中PS部分的ARM處理器,將數據通過PL端的直接內存訪問(DMA)控制器傳輸到FPGA加速模塊中進行處理。其中FPGA模塊使用高層綜合(HLS)搭建實現,該部分主要對ILA中的大型矩陣運算進行循環優化。當加速模塊運算結束后,將處理結果傳回PS端。在此過程中,由PS部分的ARM處理器對整個流程進行控制,除此之外,PS端還負責圖像重構計算過程的流程控制,包括數據輸入、DMA與FPGA的配置以及數據的輸出,其流程如圖10所示。

通過對氣-水兩相流進行相關參數的采集,并將數據采集的結果交給ZYNQ進行硬件加速,處理后的結果通過上位機顯示出來,并最終得到重構圖像。本文測量的介質分布模型如圖11所示。

實驗中,綜合考慮圖像質量,ILA的迭代次數為300次;進行圖像重構的總用時在230~360ms,而ZYNQ平臺的總用時為150~225ms。以計算機的重構時間減去ZYNQ重構時間Δt在計算機重構用時中所占比例來衡量硬件加速的加速效果。在ZYNQ以及在計算機進行數據處理后的各重構圖像如表5所示。

從表5可以看到,ZYNQ和計算機重構出的圖像有著一定的差距,這主要是因為:在ZYNQ中所參與運算的數據是16位的,而計算機上參與運算的數據是32位的,這就導致了圖像重構算法運行過程中會產生截斷誤差,進而影響了圖像重構的效果;其次在ZYNQ上固定的迭代因子α以及閾值的選取也對圖像重構的結果造成影響。

為了評價生成圖像的質量,引入圖像相關系數cimg、圖像相對誤差eimg[23]兩種圖像重構質量指標。圖像相關系數主要用于評價兩個圖像之間的相似程度,表達式為

圖像相對誤差eimg可以衡量重構介質分布與原始的介質分布之間的差異,表達式為

式中:greco和greal分別為重構和真實灰度向量。在實際使用中,式(10)中的范數為向量的二范數。

通過這些評價指標對ILA重構圖像和改進后重構圖像的質量進行評價,結果如表6所示。可以看到,在數據評價上,除層流以外其他流型的圖像相似度比較高,圖像相對誤差較低,證明了ZYNQ平臺運行圖像重構運算的可靠性。

由于數據位寬和迭代系數等問題,使用ZYNQ進行圖像重構仍然與計算機圖像重構有一定差別,但是根據結果可以看出誤差尚在可接受范圍內。在流型識別以及其他方面的應用中會大大縮減時間成本,提高效率。

需要說明的是,MATLAB軟件的求解器內置對矩陣求解的算法優化流程,相反嵌入式平臺受限于芯片資源,只能通過自定義硬件電路來實現算法加速。本文FPGA+ARM的硬件加速方案經過實驗對照,比使用MATLAB編程的優化效果更加顯著。所以除了便攜性外,在相似資源的使用下,FPGA+ARM方案仍可以提供更高的計算性能和效率,這是由其深度定制的并行架構帶來的優勢,仍具有在特定應用場景下替代甚至超過目前通用方案的潛力。

5 結 論

本文設計了一種在ARM+FPGA構架上的圖像重構加速技術和嵌入式ECT系統,并對該系統性能做出了評價,得到主要研究成果如下。

(1)在圖像重構算法方面,本文沿用了經典的Landweber算法,在FPGA的協同下,對其循環結構進行了考慮依賴性的流水線加速,提出3種優化策略。結果表明,改進反依賴關系后的加速策略運算時間減少50%以上, 優化效果比改進前得到很大提升。

(2)在硬件設計部分,本文針對ECT系統的圖像重構算法運行時間較長的問題進行了改進,開發了一種基于ZYNQ的ECT硬件加速系統。搭建了ECT氣-水兩相流流型顯示系統實驗平臺,并進行實驗驗證。實驗結果表明,該硬件加速系統在使用ILA算法進行成像時每個圖像所需運行時間比使用MATLAB編程運行的用時減少了30%~40%。由此證明了利用FPGA進行大規模矩陣運算的高效性和快速性。

參考文獻:

[1]COSTA M G, LEITE J M, BECKEDORFF L, et al. Static pressure behavior of gas-liquid flows along a Venturi [J]. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 2021, 43(11): 498.

[2]ABRAR U, YOUSAF A, JAFFRI N R, et al. Analysis of complex solid-gas flow under the influence of gravity through inclined channel and comparison with real-time dual-sensor system [J]. Electronics, 2021, 10(22): 2849.

[3]JEANMEURE L F C, DYAKOWSKI T, ZIMMERMAN W B J, et al. Direct flow-pattern identification using electrical capacitance tomography [J]. Experimental Thermal and Fluid Science, 2002, 26(6/7): 763-773.

[4]WANG Shengnan, YE Jiamin, YANG Yunjie. Quantitative measurement of two-phase flow by electrical capacitance tomography based on 3D coupling field simulation [J]. IEEE Sensors Journal, 2021, 21(18): 20136-20144.

[5]TIAN Wenbin, SUO Peng, LIU Dong, et al. Simultaneous shape and permittivity reconstruction in ECT with sparse representation: two-phase distribution imaging [J]. IEEE Transactions on Instrumentation and Measurement, 2021, 70: 1-14.

[6]PENG Lihui, YE Jiamin, LU Geng, et al. Evaluation of effect of number of electrodes in ECT sensors on image quality [J]. IEEE Sensors Journal, 2012, 12(5): 1554-1565.

[7]張順利, 張定華, 李明君, 等. 基于SIMD技術的錐束ART算法快速并行圖像重建 [J]. 儀器儀表學報, 2010, 31(3): 630-634.

ZHANG Shunli, ZHANG Dinghua, LI Mingjun, et al. Fast parallel image reconstruction with cone-beam ART algorithm based on SIMD technology [J]. Chinese Journal of Scientific Instrument, 2010, 31(3): 630-634.

[8]顏華, 孫燦烽, 王伊凡. 基于改進SIRT算法的聲學CT溫度場重建 [J]. 沈陽工業大學學報, 2021, 43(5): 550-556.

YAN Hua, SUN Canfeng, WANG Yifan. Temperature field reconstruction of acoustic CT based on improved SIRT algorithm [J]. Journal of Shenyang University of Technology, 2021, 43(5): 550-556.

[9]王化祥, 王超, 陳磊. 基于Landweber迭代的圖像重建算法 [J]. 信號處理, 2000, 16(4): 354-356.

WANG Huaxiang, WANG Chao, CHEN Lei. An image reconstruction algorithm based on the Landweber iteration method [J]. Journal of Signal Processing, 2000, 16(4): 354-356.

[10]李秀艷, 韓倩, 汪劍鳴, 等. 基于改進共軛梯度法的ERT圖像重建 [J]. 儀器儀表學報, 2016, 37(7): 1673-1679.

LI Xiuyan, HAN Qian, WANG Jianming, et al. ERT image reconstruction based on improved CG method [J]. Chinese Journal of Scientific Instrument, 2016, 37(7): 1673-1679.

[11]REAL R, JIN Qinian. A revisit on Landweber iteration [J]. Inverse Problems, 2020, 36(7): 075011.

[12]嚴春滿, 穆哲, 張道亮, 等. 基于改進Landweber算法的ECT圖像重建 [J]. 傳感技術學報, 2019, 32(10): 1522-1526.

YAN Chunman, MU Zhe, ZHANG Daoliang, et al. ECT image reconstruction based on improved Landweber algorithm [J]. Chinese Journal of Sensors and Actuators, 2019, 32(10): 1522-1526.

[13]ZHANG Xuehui, WANG Huaxiang, CUI Ziqiang, et al. A novel ECT system based on FPGA and DSP [C]//Second International Conference on Innovative Computing, Information and Control (ICICIC 2007). Piscataway, NJ, USA: IEEE, 2007: 510.

[14]CASTILLO E, MORALES D P, MARTINEZ-OLMOS A, et al. Parametrized ECT processing over FPGA for a reconfigurable application [C]//2015 Conference on Design of Circuits and Integrated Systems (DCIS). Piscataway, NJ, USA: IEEE, 2015: 1-6.

[15]SHEHAZ F. An electrical capacitance tomography system for real-time process imaging [C]//Applications of Digital Image Processing XLIII. Bellingham, WA, USA: SPIE, 2020: 1151036.

[16]ALLAM A, DEABES W. Model-based hardware-software codesign of ECT digital processing unit [J]. Modelling and Simulation in Engineering, 2021, 2021: 4757464.

[17]DEABES W. FPGA implementation of ECT digital system for imaging conductive materials [J]. Algorithms, 2019, 12(2): 28.

[18]MERIBOUT M, TENIOU S. A pipelined parallel hardware architecture for 2-D real-time electrical capacitance tomography imaging using interframe correlation [J]. IEEE Transactions on Very Large Scale Integration (VLSI) Systems, 2017, 25(4): 1320-1328.

[19]MERIBOUT M, SAIED I M. Real-time two-dimensional imaging of solid contaminants in gas pipelines using an electrical capacitance tomography system [J]. IEEE Transactions on Industrial Electronics, 2017, 64(5): 3989-3996.

[20]MA Min, YA Chaoqi. Analysis of gas/solid two-phase flow imaging system based on ECT technology [C]//2017 IEEE International Conference on Imaging Systems and Techniques (IST). Piscataway, NJ, USA: IEEE, 2017: 1-6.

[21][JP+3]CHEN Dixiang, YANG [JP+2]Wuqiang, PAN Mengchun. [LL]Design of impedance measuring circuits based on phase-sensitive demodulation technique [J]. IEEE Transactions on Instrumentation and Measurement, 2011, 60(4): 1276-1282.

[22]唐凱豪. 用于氣水兩相流可視化監測的電容層析成像技術關鍵問題研究 [D]. 西安: 西安交通大學, 2022.

[23]王化祥. 電學層析成像 [M]. 北京: 科學出版社, 2013.

[24]唐凱豪, 胡紅利, 李林, 等. 利用場量旋轉變換的電容層析成像靈敏度系數簡易計算方法 [J]. 西安交通大學學報, 2019, 53(3): 75-80.

TANG Kaihao, HU Hongli, LI Lin, et al. Simplified method for electrical capacitance tomography sensitivity coefficient computation with specific-electrode-excited field quantities rotation transformation [J]. Journal of Xi’an Jiaotong University, 2019, 53(3): 75-80.

主站蜘蛛池模板: 一级毛片高清| 欧美激情视频在线观看一区| 全部毛片免费看| 爱色欧美亚洲综合图区| 全部免费特黄特色大片视频| 国产激爽爽爽大片在线观看| 无码人中文字幕| 毛片免费高清免费| 2021精品国产自在现线看| 国产在线八区| 久久国产精品娇妻素人| 区国产精品搜索视频| 国产av色站网站| 久久天天躁狠狠躁夜夜2020一| 亚洲经典在线中文字幕| 国产亚洲高清视频| 免费毛片全部不收费的| 久久久成年黄色视频| 亚洲欧美一区二区三区麻豆| 国产无码在线调教| 天天综合亚洲| 国产精品无码久久久久AV| 日韩乱码免费一区二区三区| 91精品国产自产91精品资源| 欧美成人亚洲综合精品欧美激情| 99在线视频免费观看| 国产成人91精品免费网址在线| 好久久免费视频高清| 日本在线欧美在线| 国产一在线观看| 91www在线观看| 午夜啪啪福利| 青青草原国产av福利网站| 久久久亚洲国产美女国产盗摄| 国产小视频a在线观看| 欧美一级在线看| 亚洲Av激情网五月天| 99在线观看国产| 国产精品亚洲а∨天堂免下载| 久久精品一品道久久精品| 欧美日韩亚洲综合在线观看| аv天堂最新中文在线| 久久精品人妻中文视频| 99视频精品全国免费品| 久久黄色视频影| 亚洲AV电影不卡在线观看| 日韩视频福利| 精品人妻系列无码专区久久| 亚洲成人一区二区三区| 欧美.成人.综合在线| 国产一级妓女av网站| 色婷婷在线播放| 亚洲中文无码h在线观看| 欧美一区二区三区香蕉视| 亚洲一道AV无码午夜福利| 亚洲精品动漫在线观看| 少妇人妻无码首页| 四虎永久在线视频| 中字无码精油按摩中出视频| 中日韩一区二区三区中文免费视频| 高潮毛片无遮挡高清视频播放 | 亚洲欧美国产高清va在线播放| 日韩在线播放中文字幕| 中日韩一区二区三区中文免费视频 | 国产成人久久综合一区| 色综合a怡红院怡红院首页| 国产亚洲精久久久久久无码AV| 国产内射一区亚洲| 亚洲午夜久久久精品电影院| 亚洲人成人无码www| 人人看人人鲁狠狠高清| 久久精品人人做人人爽电影蜜月| 久久综合亚洲色一区二区三区| 国产精品成人啪精品视频| 久久国产拍爱| 国产亚洲一区二区三区在线| 欧美成人精品高清在线下载| 亚洲一级无毛片无码在线免费视频| 欧美不卡在线视频| 成人日韩视频| 天堂岛国av无码免费无禁网站| 日韩精品高清自在线|