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

基于熱紅外成像的坡面薄層水流流速測(cè)量方法

2022-01-27 03:22:48史海靜郭明航展小云丁成琴
關(guān)鍵詞:測(cè)量

張 艷,史海靜,郭明航,趙 軍,展小云,丁成琴

基于熱紅外成像的坡面薄層水流流速測(cè)量方法

張 艷1,3,史海靜1,2,3※,郭明航1,2,趙 軍1,2,展小云2,丁成琴2

(1. 中國科學(xué)院水利部水土保持研究所黃土高原土壤侵蝕與旱地農(nóng)業(yè)國家重點(diǎn)實(shí)驗(yàn)室,楊凌 712100;2. 西北農(nóng)林科技大學(xué)水土保持研究所,楊凌 712100; 3. 中國科學(xué)院大學(xué),北京 100049)

流速是表征水流水力學(xué)特性的重要物理量。為了準(zhǔn)確獲取薄層水流流速,基于熱紅外成像技術(shù)、計(jì)算機(jī)視覺識(shí)別技術(shù),設(shè)計(jì)了一種薄層水流流速測(cè)量系統(tǒng)。該系統(tǒng)通過對(duì)熱示蹤劑的自動(dòng)控制以及對(duì)其熱成像圖的瞬時(shí)采集、影像校正、噪點(diǎn)去除、質(zhì)心確定等手段,獲取薄層水流的流速等參數(shù),從而實(shí)現(xiàn)對(duì)坡面薄層水流流速的動(dòng)態(tài)觀測(cè)。該系統(tǒng)精確度和準(zhǔn)確度高,可從不同時(shí)間和空間尺度上更加準(zhǔn)確地觀測(cè)熱示蹤劑動(dòng)態(tài)運(yùn)移過程。系統(tǒng)的測(cè)量標(biāo)準(zhǔn)差為0.020 m/s,觀測(cè)精度可達(dá)到98.33%,觀測(cè)的時(shí)間分辨率為1/9 s,空間分辨率為2 mm。為了驗(yàn)證該系統(tǒng)的準(zhǔn)確度,以流量法為基準(zhǔn),與傳統(tǒng)示蹤法(染料示蹤法、鹽示蹤法)比較,結(jié)果表明,該系統(tǒng)的準(zhǔn)確度高于染料示蹤法和鹽示蹤法。其中熱紅外成像觀測(cè)系統(tǒng)的相對(duì)誤差均在±10%以內(nèi);染料示蹤法的相對(duì)誤差均大于10%;鹽示蹤法52%的相對(duì)誤差在±10%以內(nèi)。利用熱紅外成像系統(tǒng)獲取的影像,可對(duì)不同時(shí)刻示蹤段水流的發(fā)生、發(fā)展的過程進(jìn)行追溯,也可以測(cè)量示蹤劑沿水流方向的運(yùn)移速度及垂直于水流方向的擴(kuò)散速度,計(jì)算熱示蹤劑的彌散系數(shù)等。該技術(shù)可應(yīng)用于降雨侵蝕、徑流沖刷等方面的研究,對(duì)于進(jìn)一步深化土壤侵蝕過程與機(jī)理研究具有重要的意義。

流速;坡面;熱紅外成像觀測(cè)系統(tǒng);薄層水流;熱紅外標(biāo)靶;精準(zhǔn)度評(píng)價(jià)

0 引 言

流速作為水動(dòng)力學(xué)重要的參數(shù)之一,是水土流失預(yù)測(cè)模型建立的關(guān)鍵因素,且坡面流速的準(zhǔn)確測(cè)量對(duì)于提高模型精度和定量分析徑流攜沙能力等具有重要意義[1]。而坡面水流不同于明渠水流,水深較淺一般為毫米量級(jí),稱為薄層水流。它通常在降雨或集中徑流沖刷條件下形成,且容易受到坡度、流量、下墊面狀況等的影響[2-4]。因此,對(duì)于坡面薄層水流流速的準(zhǔn)確測(cè)量是目前土壤侵蝕定量研究中的難點(diǎn)之一。

常用的坡面水流流速測(cè)量方法主要包括:示蹤法[5-8]、流量法[9-10]、旋漿式微流速測(cè)定法[11-13]、粒子圖像法[14-16]以及20世紀(jì)60年代發(fā)展起來的相關(guān)法[17-19]等。傳統(tǒng)的示蹤法包括染料示蹤法和電解質(zhì)示蹤法,示蹤法在流速和泥沙含量較小時(shí)能很好地測(cè)量水流速度,但因示蹤劑自身擴(kuò)散和觀察困難的問題使得測(cè)量的流速不準(zhǔn)確,且在計(jì)算平均速度時(shí)都需要引用經(jīng)驗(yàn)參數(shù)[6]。流量法基于流量水深比測(cè)量水流流速,通常適用于明渠水流,對(duì)于非常規(guī)則的渠道,是目前測(cè)量流速較為準(zhǔn)確的方法之一,但只能測(cè)量斷面的平均流速[10]。旋漿式微流速測(cè)定法一般也用于明渠水流,但漿的旋轉(zhuǎn)對(duì)水流有干擾作用,因此也不能準(zhǔn)確測(cè)量薄層水流流速[13]。粒子圖像法可以準(zhǔn)確測(cè)量清水流某一點(diǎn)的流速,但當(dāng)水流中含有其他物質(zhì)時(shí)會(huì)影響粒子的傳播,降低測(cè)量的準(zhǔn)確度[14]。相關(guān)法一般用于兩相流、多相流的測(cè)定,但在測(cè)量過程中難以區(qū)分清水速度和懸浮物速度,且在水流中放置傳感器會(huì)干擾水流流速的測(cè)定[17]。以上的觀測(cè)技術(shù)和方法,對(duì)薄層水流流速的測(cè)量存在干擾大、缺乏動(dòng)態(tài)監(jiān)測(cè)等問題。要么只能獲取某一個(gè)時(shí)間段的平均流速,要么只能獲取某一個(gè)空間點(diǎn)的瞬時(shí)流速,無法準(zhǔn)確估算示蹤段水流的整體情況。近年來,熱紅外成像技術(shù)[20-22]逐漸被應(yīng)用于薄層水流流速的研究中,為連續(xù)觀測(cè)薄層水流形態(tài)以及流速空間分布的數(shù)字化表達(dá)提供了可能。但對(duì)于熱示蹤劑的定量化控制、熱成像影像解譯坐標(biāo)系的建立,以及影像解算及處理等方面并不完善,并沒有形成一個(gè)完整的應(yīng)用觀測(cè)系統(tǒng)。

綜上,本文擬耦合熱紅外成像和計(jì)算機(jī)視覺識(shí)別技術(shù),研制一套對(duì)薄層水流流速進(jìn)行實(shí)時(shí)、動(dòng)態(tài)、連續(xù)觀測(cè)的系統(tǒng)。利用該系統(tǒng)對(duì)熱示蹤劑進(jìn)行自動(dòng)的控制,并對(duì)其熱成像圖進(jìn)行同步采集,獲取具有高時(shí)空分辨率的影像,通過解算得到不同時(shí)間點(diǎn)和距離點(diǎn)的薄層水流流速,從而為研究坡面土壤侵蝕過程、機(jī)理、建模等提供理論依據(jù)和研究基礎(chǔ)。

1 熱紅外成像觀測(cè)系統(tǒng)的工作原理與系統(tǒng)組成

1.1 工作原理

熱紅外成像觀測(cè)系統(tǒng)的工作原理主要是基于熱紅外成像和計(jì)算機(jī)視覺識(shí)別技術(shù)完成對(duì)熱示蹤劑瞬態(tài)遷移過程的監(jiān)測(cè),利用影像的表面溫度分布圖以及溫度差來區(qū)分背景水流和熱示蹤劑,從而對(duì)熱示蹤劑的運(yùn)移過程進(jìn)行解算。首先利用熱紅外攝像機(jī)獲取坡面薄層水流熱示蹤劑的運(yùn)移視頻,通過無線網(wǎng)絡(luò)對(duì)視頻數(shù)據(jù)進(jìn)行傳輸并按幀率提取出影像;然后利用熱紅外標(biāo)靶建立坐標(biāo)系,通過對(duì)影像的校正、噪點(diǎn)去除、質(zhì)心確定等手段,獲取薄層水流的流速等參數(shù),從而實(shí)現(xiàn)對(duì)坡面薄層水流流速的動(dòng)態(tài)觀測(cè)。

1.2 系統(tǒng)組成

熱紅外成像觀測(cè)系統(tǒng)由熱示蹤劑控制子系統(tǒng)、影像采集與傳輸子系統(tǒng)以及影像解算子系統(tǒng)組成,且各子系統(tǒng)由不同的軟硬件單元組成(圖1)。其中熱示蹤劑控制子系統(tǒng)用于制備和添加熱示蹤劑;影像采集與傳輸子系統(tǒng)包括熱紅外攝像機(jī)、標(biāo)靶以及無線路由器,主要用于熱示蹤劑瞬態(tài)遷移過程的監(jiān)測(cè)以及數(shù)據(jù)的傳輸;影像解算子系統(tǒng)用于識(shí)別和確定熱示蹤劑的前緣點(diǎn)和質(zhì)心點(diǎn)從而計(jì)算水流流速。

圖1 熱紅外成像觀測(cè)系統(tǒng)結(jié)構(gòu)

1.2.1 熱示蹤劑控制系統(tǒng)

熱示蹤劑控制系統(tǒng)主要負(fù)責(zé)提供持續(xù)恒溫的熱水源,該系統(tǒng)由電加熱器、溫度傳感器和熱水泵組成。電加熱器用以加熱水體,溫度傳感器用以感應(yīng)熱示蹤劑的溫度并調(diào)節(jié)水溫至設(shè)定溫度,熱水泵用以調(diào)節(jié)熱示蹤劑滴入量和時(shí)間間隔,從而定時(shí)定量地滴入熱示蹤劑。對(duì)加入示蹤劑的溫度和體積多次試驗(yàn),設(shè)定熱示蹤劑初始溫度為100℃,出水時(shí)間為0.2 s(10 mL),間隔時(shí)間為9 s。

1.2.2 影像采集與傳輸系統(tǒng)

影像采集與傳輸系統(tǒng)主要負(fù)責(zé)熱示蹤劑運(yùn)移影像的采集和傳輸,包括1臺(tái)FLIR ONE 3.1.0 熱紅外攝像機(jī)、4個(gè)熱紅外標(biāo)靶和1個(gè)無線路由器。其中,F(xiàn)LIR ONE 3.1.0熱紅外相機(jī)的分辨率為 4 032′3 016像素,相機(jī)幀率為 9 幀/s,固定于水槽上方1.5 m處進(jìn)行影像采集,拍攝水槽的范圍為1.2 m(長)′0.8 m(寬);熱紅外標(biāo)靶主要用于坐標(biāo)系的建立;無線路由器負(fù)責(zé)將拍攝的視頻數(shù)據(jù)傳入存儲(chǔ)系統(tǒng)。

1.2.3 影像解算系統(tǒng)

影像解算系統(tǒng)主要負(fù)責(zé)提取高幀率熱成像影像并進(jìn)行解算。主要由數(shù)據(jù)存儲(chǔ)單元和計(jì)算匹配單元組成。影像解算系統(tǒng)的軟件有2個(gè)模塊:配置模塊和計(jì)算模塊。配置模塊主要實(shí)現(xiàn)輸入?yún)?shù)、輸出參數(shù)、計(jì)算參數(shù)、作者和版權(quán)信息參數(shù)等的配置。計(jì)算模塊主要實(shí)現(xiàn)影像分解、圖片校正、物體識(shí)別、流速計(jì)算、結(jié)果輸出等功能。在模塊的開發(fā)過程中采用JAVA、GDAL、OPENCV軟件。

2 熱紅外成像觀測(cè)系統(tǒng)測(cè)量流程與關(guān)鍵技術(shù)

2.1 測(cè)量流程

使用熱紅外成像觀測(cè)系統(tǒng)測(cè)量水流流速的具體流程如下:1)利用熱示蹤劑控制系統(tǒng)將熱示蹤劑添加到水流中;2)采用熱紅外相機(jī)對(duì)熱示蹤劑進(jìn)行拍攝;3)對(duì)獲取的熱紅外影像解算得到熱示蹤劑的流速等參數(shù)?;谝陨喜襟E生成熱示蹤劑運(yùn)移的每一幀影像,都可以解算出相應(yīng)的前緣速度、質(zhì)心速度等物理量。其流程圖如圖2所示。

圖2 影像解算流程

2.2 關(guān)鍵技術(shù)

2.2.1 坐標(biāo)系的構(gòu)建

坐標(biāo)系建立的目的是將所有影像放在統(tǒng)一坐標(biāo)系下進(jìn)行解算,建立的關(guān)鍵是要有可識(shí)別的原點(diǎn)和標(biāo)靶點(diǎn)。由于熱紅外成像只能通過溫度來感應(yīng)和識(shí)別物體,因此,在熱紅外相機(jī)拍攝中普通光學(xué)的標(biāo)靶未能成像或成像模糊。本文根據(jù)鞏稼民等[23]提出的熱紅外標(biāo)靶,利用電熱絲制作了“田”字格熱紅外標(biāo)靶,但該方法因電熱絲溫度控制困難且比熱容小,易導(dǎo)致熱紅外標(biāo)靶出現(xiàn)成像不一致、輪廓模糊等問題。如圖3所示,圖3a為開始加熱的“田”字格熱紅外標(biāo)靶、圖3b為溫度恒定的“田”字格熱紅外標(biāo)靶。從圖3a和圖3b中可以明顯看出,標(biāo)靶成像不一致且隨著溫度的升高,電熱絲發(fā)熱成像逐漸模糊,出現(xiàn)熱擴(kuò)散現(xiàn)象。為解決該問題,本研究選用比熱容較大的水,研制了熱紅外標(biāo)靶,即圓形熱水標(biāo)靶。將熱水加入形狀規(guī)則的圓柱體鋁盒中(圖3d),由于水的比熱容較大,在成像時(shí)能夠形成清晰的邊界,因此,用熱水制作的圓形標(biāo)靶具有輪廓清晰、成像一致的優(yōu)點(diǎn)(圖3c)。該標(biāo)靶的直徑為5.5 cm,高為3 cm。為了保證標(biāo)靶里熱水的溫度恒定,與影像獲取時(shí)間同步,每2 min更換一次熱水,確保圓形標(biāo)靶圖像清晰。熱紅外標(biāo)靶確定后,在水槽的兩側(cè)選定測(cè)量區(qū)域放置4個(gè)標(biāo)靶,根據(jù)標(biāo)靶的已知距離來確定像素尺寸、對(duì)圖像進(jìn)行校正、確定水槽區(qū)域、確定拍攝視頻方向、標(biāo)定選定區(qū)域的坐標(biāo),從而統(tǒng)一所有影像的坐標(biāo)便于計(jì)算和分析。

2.2.2 前緣點(diǎn)的識(shí)別與確定

在薄層水流中添加熱示蹤劑,其隨水流方向運(yùn)移且運(yùn)移的最前緣為前緣點(diǎn)。前緣點(diǎn)的確定是計(jì)算前緣速度的基礎(chǔ)。而在前緣點(diǎn)確定之前,首先要在獲取的影像中確定標(biāo)靶位置以及熱區(qū)域,其過程如下:將熱紅外相機(jī)拍攝的視頻按幀率導(dǎo)出彩色圖,圖像中有紅綠藍(lán)3個(gè)波段。因標(biāo)靶和熱示蹤劑的溫度較高,所以彩色圖中紅色和綠色波段值較高。利用彩色轉(zhuǎn)灰度算法,將彩色圖轉(zhuǎn)換為灰度圖。在灰度圖上,利用霍夫圓算法查找圓,找到標(biāo)靶圓心和半徑,確定標(biāo)靶位置并建立興趣區(qū)(region of interest),對(duì)興趣區(qū)按一定算法去除噪點(diǎn)。然后將背景水流溫度的最高溫度設(shè)定為閾值溫度,溫度高于閾值溫度的像素區(qū)域?yàn)闊崾聚檮﹨^(qū)域。由于背景溫度隨周圍溫度的變化而變化,因此與影像獲取時(shí)間同步,每獲取一次影像,計(jì)算一次閾值溫度。閾值溫度的計(jì)算根據(jù)Abrantes等[22]的方法,將熱示蹤劑加入前10 s的水流中出現(xiàn)的最高溫度設(shè)置為閾值溫度,用熱示蹤劑的溫度減去閾值溫度就可確定熱示蹤劑區(qū)域。熱示蹤劑區(qū)域確定后,將其二值化保存,最后根據(jù)二值化圖識(shí)別并確定前緣點(diǎn)(圖4)。

a. 彩色圖a. Color figureb. 灰度圖b. Grayscale figurec. 興趣區(qū)c. Region of interestd. 二值化圖d. Binarization figure

2.2.3 質(zhì)心點(diǎn)的識(shí)別與確定

由于熱示蹤劑在薄層水流中的運(yùn)移形狀不規(guī)則,常規(guī)的方法難以尋找及計(jì)算質(zhì)心。本研究基于灰度或二值化輪廓圖像[24]快速計(jì)算質(zhì)心。將彩色圖轉(zhuǎn)為灰度圖后,根據(jù)溫度閾值,查找滿足要求的邊沿像素所在的相應(yīng)行列位置,然后進(jìn)行算術(shù)平均求得質(zhì)心坐標(biāo),從而確定質(zhì)心點(diǎn)(圖5)。具體的過程如下:首先從每行的左端和右端查找滿足設(shè)定閾值的第1個(gè)像素,并記錄列坐標(biāo)值y及個(gè)數(shù)m,然后從每列的上端和下端查找滿足設(shè)定閾值的第1個(gè)像素并記錄行坐標(biāo)值x及個(gè)數(shù)n,最后通過公式(1)算術(shù)平均求取熱成像圖的質(zhì)心坐標(biāo)。

式中x為行坐標(biāo)值;y為列坐標(biāo)值;n為行坐標(biāo)個(gè)數(shù);m為列坐標(biāo)個(gè)數(shù);、、、為分別為目標(biāo)圖像的行數(shù)和列數(shù)。

2.2.4 前緣速度和質(zhì)心速度的計(jì)算

熱紅外成像技術(shù)的原理通過計(jì)算熱示蹤劑在單位時(shí)間內(nèi)的移動(dòng)距離計(jì)算水流流速,前緣速度和質(zhì)心速度的計(jì)算見式(2)和式(3)。熱紅外成像技術(shù)獲取的影像為高幀率影像,本研究中使用的熱紅外相機(jī)的幀率為9幀/s,即每1/9 s就可獲取一張熱紅外影像,并解算得到熱示蹤劑運(yùn)移的區(qū)域和距離。根據(jù)每張影像解算的熱示蹤劑區(qū)域和距離可計(jì)算得到熱示蹤劑在每幀影像運(yùn)移的速度即為瞬時(shí)速度Vx。前緣瞬時(shí)速度為熱示蹤劑在單位時(shí)間內(nèi)最前緣點(diǎn)的運(yùn)移距離,質(zhì)心瞬時(shí)速度為熱示蹤劑在單位時(shí)間內(nèi)質(zhì)心點(diǎn)的運(yùn)移距離(圖6)。由于流量法計(jì)算得到的是平均流速,為了基于流量法驗(yàn)證本研究測(cè)量流速的準(zhǔn)確性,將前緣瞬時(shí)速度和質(zhì)心瞬時(shí)速度分別算術(shù)平均得到水流平均流速,見式(4),進(jìn)而對(duì)流量法和熱紅外成像技術(shù)測(cè)量的流速進(jìn)行對(duì)比研究。

前緣速度:

質(zhì)心速度:

平均流速:

式中D2為從1時(shí)刻到2時(shí)刻示蹤劑前緣運(yùn)移的路程,m;D1為從1時(shí)刻到2時(shí)刻示蹤劑質(zhì)心運(yùn)移的路程,m;Vx為示蹤劑的瞬時(shí)速度,為一次測(cè)量所解譯的瞬時(shí)速度的個(gè)數(shù)。

注:1為熱示蹤劑質(zhì)心點(diǎn)從1到2時(shí)刻的運(yùn)移距離,m;2為熱示蹤劑前緣點(diǎn)從1到2時(shí)刻的運(yùn)移距離,m;1和2為時(shí)間,s

Note:1is the migration distance of thermal tracer centroid point from1to2, m;2is the migration distance of the leading edge point of the thermal tracer from1to2, m;1and2are the times, s.

圖6 熱示蹤劑運(yùn)移圖

Fig.6 Thermal tracer migration diagram

3 熱紅外成像觀測(cè)系統(tǒng)的應(yīng)用

3.1 試驗(yàn)過程

3.1.1 試驗(yàn)設(shè)計(jì)

試驗(yàn)在黃土高原土壤侵蝕與旱地農(nóng)業(yè)國家重點(diǎn)實(shí)驗(yàn)室人工模擬徑流大廳進(jìn)行。試驗(yàn)平臺(tái)規(guī)格為10 m(長)′3 m(寬),試驗(yàn)水槽為可調(diào)坡度水槽(圖7),規(guī)格為 4.6 m(長)′0.196 m(寬)′0.1 m(深),設(shè)置其下墊面為光滑的玻璃。光滑玻璃的表面粗糙度為0.03m[25],可近似看作粗糙度為0,在形狀規(guī)則的水槽中可通過流量法精準(zhǔn)計(jì)算水流流速。在水槽進(jìn)水口處設(shè)置穩(wěn)流槽,使水流從穩(wěn)流槽中溢流而出,并且均勻地分布在坡面上。在水槽出水口處放置徑流箱,利用水泵將徑流箱中的水流傳輸?shù)焦┧b置中,以達(dá)到水流循環(huán)使用的效果。在試驗(yàn)開始前,調(diào)節(jié)供水裝置的流量,多次試驗(yàn)取樣證明流量設(shè)計(jì)值與實(shí)際流量值一致。在放水沖刷過程中,利用熱紅外攝像機(jī)對(duì)熱示蹤劑進(jìn)行拍攝,同時(shí)對(duì)獲取的影像數(shù)據(jù)進(jìn)行傳輸存儲(chǔ),防止數(shù)據(jù)的丟失和混亂,最后利用影像解譯系統(tǒng)對(duì)所獲取的影像進(jìn)行解譯。

圖7 試驗(yàn)水槽

3.1.2 精確度檢驗(yàn)

精確度是指對(duì)某一對(duì)象多次重復(fù)測(cè)量數(shù)值的離散表達(dá),離散小則精度高。為檢驗(yàn)熱紅外成像觀測(cè)系統(tǒng)的測(cè)量精度,在水槽坡度為5°,流量為0.2 L/s條件下重復(fù)30次,總時(shí)長為5 min,進(jìn)行水流流速觀測(cè)。以熱示蹤劑的質(zhì)心速度為觀測(cè)對(duì)象,因?yàn)樵诜潘疀_刷條件下,示蹤劑峰值溫度的位置較難確定,因此將熱示蹤劑的質(zhì)心速度作為水流的平均流速[26],利用SPSS 18.0軟件對(duì)30次的重復(fù)測(cè)量結(jié)果進(jìn)行統(tǒng)計(jì)分析。

3.1.3 準(zhǔn)確度檢驗(yàn)

準(zhǔn)確度是表征測(cè)量值與真值之間的關(guān)系,用相對(duì)誤差來表述。在放水沖刷條件下,通過3個(gè)坡度(0、5、10°)、7種流量(0.2、0.3、0.4、0.5、0.6、1.0、1.5 L/s)共21種水流狀況下,每種水流狀況下拍攝時(shí)長為2 min。以流量法為基準(zhǔn),分別計(jì)算熱紅外成像技術(shù)、染料示蹤法、鹽示蹤法測(cè)量的相對(duì)誤差,從而評(píng)估熱紅外成像技術(shù)、染料示蹤法、鹽示蹤法對(duì)水流流速觀測(cè)的準(zhǔn)確度。流量法是目前測(cè)量平均流速較準(zhǔn)確的方法之一,因此采用流量法測(cè)量的流速值作為真值,其計(jì)算方法為式(5)。染料示蹤法是將高錳酸鉀溶液加入被測(cè)薄層水流中,記錄其通過測(cè)量區(qū)域的時(shí)間,然后根據(jù)手工測(cè)量的距離和記錄時(shí)間計(jì)算薄層水流流速,其中高錳酸鉀溶液配制為每500 mL去離子水中加入5 g高錳酸鉀,該方法測(cè)量得到的是水流的前緣速度。鹽示蹤法是在薄層水流中滴入鹽溶液,并在其下游另一個(gè)確定位置對(duì)水流的電導(dǎo)率進(jìn)行測(cè)量,通過計(jì)算鹽溶液滴入薄層水流至電導(dǎo)率最大值出現(xiàn)之間的時(shí)間、水流流過的距離計(jì)算薄層水流流速,其中鹽溶液的濃度為每1 000 mL去離子水中加入5 g的食用鹽,時(shí)間用秒表計(jì)時(shí)。為保證試驗(yàn)的可靠性,對(duì)坡度、流量的每種組合條件下進(jìn)行3次重復(fù)。該方法測(cè)量得到的是水流的質(zhì)心速度。

式中為過流斷面流體體積,L/s;為水流水深,本文用測(cè)針儀測(cè)量精度為0.1 mm;為水槽寬度,本文為0.196 m。

3.2 結(jié)果與分析

3.2.1 精確度分析

精確度檢驗(yàn)結(jié)果表明,測(cè)量的熱示蹤劑質(zhì)心速度平均值為0.517 m/s,測(cè)量標(biāo)準(zhǔn)差為0.020 m/s,觀測(cè)精度可達(dá)到98.33%。利用單一樣本K-S(Kolmogorov-Smirnov)檢驗(yàn)該質(zhì)心速度的分布,統(tǒng)計(jì)分析得到K-S的值為0.505,值為0.961(>0.05),由此可知熱紅外成像技術(shù)的觀測(cè)結(jié)果服從(0.517, 0.020)正態(tài)分布(如圖8所示)。

注:m為平均值;d為標(biāo)準(zhǔn)差;N為樣本數(shù)。

3.2.2 準(zhǔn)確度分析

準(zhǔn)確度檢驗(yàn)結(jié)果表明(表1、圖9),熱紅外成像觀測(cè)系統(tǒng)的準(zhǔn)確度高于染料示蹤法和鹽示蹤法。熱紅外成像技術(shù)測(cè)量的質(zhì)心速度與實(shí)際值之間的最大相對(duì)誤差、最小相對(duì)誤差分別為–9.61%和0.16%,且相對(duì)誤差范圍均在±10%以內(nèi)(表1、圖9a)。染料示蹤法測(cè)量的流速值與實(shí)際值之間的最大相對(duì)誤差為75.91%,最小相對(duì)誤差為12.03%,且染料示蹤法測(cè)量的流速值均大于流量法測(cè)量的流速值(表1、圖9c)。鹽示蹤法流速測(cè)量值與實(shí)際值之間的最大相對(duì)誤差為26.67%,最小相對(duì)誤差為0.11%,相對(duì)誤差在±10%以內(nèi)的樣本數(shù)有52%(表1、圖9d)。

以流量法為基準(zhǔn),熱紅外成像技術(shù)測(cè)得的前緣速度的最大相對(duì)誤差、最小相對(duì)誤差分別為28.05%、–4.81%,且只有33%的樣本相對(duì)誤差在±10%以內(nèi)(表1、圖9b)。這表明熱紅外成像技術(shù)的質(zhì)心速度的準(zhǔn)確度高于前緣速度的準(zhǔn)確度。其中質(zhì)心速度的相對(duì)誤差均小于前緣速度的相對(duì)誤差,這與de Lima等[21-22]的研究一致,因?yàn)榱髁糠y(cè)量的流速為水流的平均流速,所以在熱紅外成像技術(shù)測(cè)速時(shí)用熱示蹤劑的質(zhì)心速度來表征水流平均流速。而染料示蹤法測(cè)量的也是前緣速度,其觀測(cè)樣本有95%的相對(duì)誤差大于熱紅外成像技術(shù)前緣速度的相對(duì)誤差,這表明熱示蹤劑的擴(kuò)散小于染料示蹤劑的擴(kuò)散。且由于示蹤段內(nèi)人腦的反應(yīng)等主觀原因,時(shí)間記錄有偏差,導(dǎo)致染料示蹤法測(cè)得的流速值誤差較大。

表1 不同方法流速測(cè)量結(jié)果

圖9 水流流速相對(duì)誤差頻率分布

無論是染料示蹤法還是鹽示蹤法,在測(cè)量過程中由于人為記錄,只能讀取一個(gè)時(shí)間值,無法測(cè)算示蹤段水流的整體情況。而熱紅外成像技術(shù)可以準(zhǔn)確地記錄示蹤段水流的整體情況,并且通過影像解算,可在不同的時(shí)間節(jié)點(diǎn)和空間點(diǎn)上識(shí)別示蹤劑所形成圖斑的前緣點(diǎn)和質(zhì)心點(diǎn),從而準(zhǔn)確地計(jì)算水流的前緣速度以及質(zhì)心速度。本研究中,熱紅外成像技術(shù)觀測(cè)的時(shí)間分辨率為1/9 s,空間分辨率為2 mm,實(shí)現(xiàn)了從不同時(shí)空尺度上對(duì)薄層水流流速的觀測(cè)(圖10)。

注:圖中紅點(diǎn)為質(zhì)心點(diǎn),紅線為前緣線。坡度為0、流量為0.2 L·s-1。

相比于染料示蹤法和鹽示蹤法,熱紅外成像觀測(cè)技術(shù)的主要優(yōu)點(diǎn)在于,熱圖像中示蹤劑具有較高的可視化且前緣的能見度高,可以準(zhǔn)確尋找前緣點(diǎn)和示蹤劑區(qū)域。選定示蹤劑區(qū)域后,可以快速確定前質(zhì)心,計(jì)算示蹤劑的前緣速度和質(zhì)心速度。在染料示蹤法中,示蹤劑的邊緣難以確定,即使利用光學(xué)相機(jī)對(duì)染料示蹤劑進(jìn)行拍攝,也很難估算示蹤劑的質(zhì)心,尤其在野外工作中,因?yàn)楣鈱W(xué)相機(jī)對(duì)光照條件嚴(yán)格且光學(xué)圖像處理復(fù)雜。在鹽示蹤法中,盡管鹽示蹤劑的保守性高,但要想實(shí)現(xiàn)對(duì)水流流速的動(dòng)態(tài)觀測(cè),只能通過安裝多個(gè)傳感器來實(shí)現(xiàn),但傳感器對(duì)薄層水流的干擾作用大,不利于水流流速的準(zhǔn)確測(cè)量。此外,安裝傳感器測(cè)量鹽溶液的濃度需要一個(gè)最小水深,導(dǎo)致小于該水深的水流流速不能測(cè)量或測(cè)量誤差大。而熱紅外成像技術(shù)可以實(shí)現(xiàn)對(duì)水流流速的無接觸式的監(jiān)測(cè)且不需要最小水深。在示蹤劑材料的選取上,染料示蹤劑通常為高錳酸鉀等帶顏色的溶液,鹽示蹤劑通常為鹽溶液,它們對(duì)環(huán)境污染大,而熱示蹤劑通常為熱水或冰塊,是一種清潔、環(huán)保的材料。

同時(shí),利用熱紅外成像技術(shù)獲取的影像,可對(duì)不同時(shí)刻示蹤段水流的發(fā)生、發(fā)展的過程進(jìn)行追溯,也可以測(cè)量示蹤劑沿水流方向的運(yùn)移速度及垂直于水流方向的擴(kuò)散速度,計(jì)算熱示蹤劑的彌散系數(shù)等。本研究中試驗(yàn)設(shè)置的下墊面為光滑的玻璃,后續(xù)將進(jìn)一步開展基于不同下墊面包括人工草皮以及土壤表面等的放水沖刷試驗(yàn),通過熱紅外成像技術(shù)來研究不同下墊面情況對(duì)薄層水流流速的影響,尤其是發(fā)生土壤侵蝕的下墊面。未來,該技術(shù)可應(yīng)用于降雨侵蝕、徑流沖刷等方面的研究,對(duì)于進(jìn)一步深化土壤侵蝕過程與機(jī)理研究具有重要的意義。

4 結(jié) 論

1)提出了一種能夠?qū)ζ旅姹铀髁魉龠M(jìn)行動(dòng)態(tài)觀測(cè)的熱紅外成像觀測(cè)系統(tǒng)。通過統(tǒng)一坐標(biāo)系的建立,高幀率影像的采集,及影像的校正、噪點(diǎn)去除、質(zhì)心確定等關(guān)鍵技術(shù),提取坡面薄層水流流態(tài)動(dòng)態(tài)變化及空間分布信息,從而實(shí)現(xiàn)對(duì)坡面薄層水流流速的動(dòng)態(tài)觀測(cè)。

2)該系統(tǒng)精確度和準(zhǔn)確度高,可從不同時(shí)間和空間尺度上更加準(zhǔn)確地觀測(cè)熱示蹤劑動(dòng)態(tài)運(yùn)移過程。系統(tǒng)的測(cè)量標(biāo)準(zhǔn)差為0.020 m/s,熱示蹤劑的質(zhì)心速度觀測(cè)精度可達(dá)到98.33%,觀測(cè)的時(shí)間分辨率為1/9 s,空間分辨率為2 mm。與傳統(tǒng)示蹤法(染料示蹤法、鹽示蹤法)比較,該系統(tǒng)的準(zhǔn)確度高于染料示蹤法和鹽示蹤法。其中熱紅外成像觀測(cè)系統(tǒng)的相對(duì)誤差均在±10%以內(nèi);染料示蹤法的相對(duì)誤差均大于10%;鹽示蹤法52%樣本的相對(duì)誤差在±10%以內(nèi)。

利用熱紅外成像技術(shù)獲取的影像,可對(duì)不同時(shí)刻示蹤段水流的發(fā)生、發(fā)展的過程進(jìn)行追溯,也可以測(cè)量示蹤劑沿水流方向的運(yùn)移速度及垂直于水流方向的擴(kuò)散速度,計(jì)算熱示蹤劑的彌散系數(shù)等。該技術(shù)可應(yīng)用于降雨侵蝕、徑流沖刷等方面的研究,對(duì)于進(jìn)一步深化土壤侵蝕過程與機(jī)理研究具有重要的意義。

[1] 劉泉. 徑流流速虛擬儀器測(cè)試系統(tǒng)的研究[D]. 武漢:華中農(nóng)業(yè)大學(xué),2005.

Liu Quan. Research on Runoff Current Velocity Measuring System of Virtual Instrument[D]. Wuhan: Huazhong Agricultural University, 2005. (in Chinese with English abstract)

[2] 白玉潔,張風(fēng)寶,楊明義,等. 急陡黃土坡面薄層水流水力學(xué)參數(shù)變化特征[J]. 土壤學(xué)報(bào),2018,55(3):641-649.

Bai Yujie, Zhang Fengbao, Yang Mingyi, et al. Variation of hydraulic parameters of shallow flow on steep loess slope[J]. Acta Pedologica Sinica, 2018, 55(3): 641-649. ( in Chinesewith English abstract)

[3] 張琪琳,王占禮,張慶瑋,等. 雨強(qiáng)及坡度對(duì)黃土區(qū)草地坡面水流流速的影響[J]. 人民黃河,2018,40(4):96-99.

Zhang Qilin, Wang Zhanli, Zhang Qingwei, et al. Experiment on influence of rainfall intensity and slope on flow velocity in loess area grassland[J]. Yellow River, 2018, 40(4): 96-99. ( in Chinese with English abstract)

[4] 陳麗燕.復(fù)雜下墊面水流流速測(cè)量與估算[D].北京:中國農(nóng)業(yè)大學(xué),2018.

Chen Yanli. Measuring and Estimating Shallow Water Flow Velocity under Complex Underlying Soil Surface Conditions[D]. Beijing: China Agricultural University, 2018. ( in Chinese with English abstract)

[5] 張光輝.坡面薄層流水動(dòng)力學(xué)特性的實(shí)驗(yàn)研究[J].水科學(xué)進(jìn)展,2002,13(2):159-165.

Zhang Guanghui. Study on hydraulic properties of shallow flow[J]. Advances in Water Science, 2002, 13(2): 159-165. (in Chinesewith English abstract )

[6] Lei T W, Xia W S, Zhao J, et al. Method for measuring velocity of shallow water flow for soil erosion with an electrolyte tracer[J]. Journal of Hydrology,2005, 301(1): 139-145.

[7] 邢行,陳曉燕,韓珍,等. 飽和與非飽和黃綿土細(xì)溝徑流水動(dòng)力學(xué)特征及侵蝕阻力對(duì)比[J]. 水土保持學(xué)報(bào),2018,32(3):92-97.

Xing Hang, Chen Xiaoyan, Han Zhen, et al. Comparation of hydrodynamic characteristics and flow resistance under rill erosion between saturated and unsaturated loess soil[J]. Journal of Soil and Water Conservation, 2018, 32(3): 92-97. (in Chinese with English abstract)

[8] Chen C, Ban Y Y, Wang X F, et al. Measuring flow velocity on frozen and non-frozen slopes of black soil through leading edge method[J]. International Soil and Water Conservation Research, 2017, 5(3): 180-189.

[9] 陳麗燕,雷廷武,啜瑞媛. 坡面薄層水流優(yōu)勢(shì)流速研究[J].中國水土保持科學(xué),2016,14(5):130-137.

Chen Liyan, Lei Tingwu, Chuo Ruiyuan. Studies on peak velocity of shallow water flow on slopes[J]. Science of Soil and Water Conservation, 2016, 14(5): 130-137. (in Chinese with English abstract)

[10] 夏衛(wèi)生,雷廷武,趙軍. 泥沙含量對(duì)鹽液示蹤法經(jīng)驗(yàn)系數(shù)影響的研究[J].農(nóng)業(yè)工程學(xué)報(bào),2003,19(4):97-100.

Xia Weisheng, Lei Tingwu, Zhao Jun. Research of empirical coefficient of salt tracer method aff ected by sediment[J].Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2003, 19(4): 97-100. (in Chinese with English abstract)

[11] 李剛,黃同愿,邵桂芳,等. 基于復(fù)雜可編程邏輯器件的流速快速測(cè)量系統(tǒng)[J]. 傳感器技術(shù),2004,23(6):23-26.

Li Gang, Huang Tongyuan, Shao Guifang, et al. High speed system to measure flow velocity based on CPLD[J]. Journal of Transducer Technology, 2004, 23(6): 23-26. (in Chinese with English abstract)

[12] 張科利,秋吉康弘,張興奇. 坡面徑流沖刷及泥沙輸移特征的試驗(yàn)研究[J]. 地理研究,1998,17 (2):3-5.

Zhang Keli, Akiyoshi Yasuhiro, Zhang Xinqi. A laboratory study on rill erosion and sediment delivery on the slope[J]. Geographical Research, 1998, 17(2): 3-5. (in Chinese with English abstract)

[13] 范建成,邢小寧. 光纖通孔式流速旋漿傳感器研制[J]. 泥沙研究1997(3):73-76.

Fan Jiancheng, Xing Xiaoning. Optical fiber miniature propeller probe with passing hole[J]. Journal of Sediment Research, 1997(3): 73-76.(in Chinese with English abstract)

[14] Akbarpour F, Fathimoghadam M, Schneider J. Application of LSPIV to measure supercritical flow in steep channels with low relative submergence[J]. Flow Measurement and Instrumentation, 2020, 72: 101718.

[15] Shi S B, Wang D W, Qian Y L, et al. Liquid-phase turbulence measurements in air-water two-phase flows using particle image velocimetry[J]. Progress in Nuclear Energy, 2020, 124: 103334.

[16] Mazumder Q, Nallamothu V T, Mazumder F. Comparison of characteristic particle velocities in solid-liquid multiphase flow in elbow[J]. International Journal of Thermofluids, 2020, 5-6: 100032

[17] 王一鳴,于海波,何建宇,等. 用于測(cè)量氣固兩相流質(zhì)量流量的梯度相關(guān)方法[J]. 中國農(nóng)業(yè)大學(xué)學(xué)報(bào),1997,2(2):105-108.

Wang Yiming, Yu Haibo, He Jianyu, et al. Correlation gradient method for measuring the mass flow in two phase flow cases[J]. Journal of China Agricultural University, 1997, 2(2): 105-108. ( in Chinese with English abstract)

[18] 宿成基,陸志明,夏虹. 利用流體中溫度噪聲信號(hào)的相關(guān)分析進(jìn)行流速測(cè)量[J]. 哈爾濱船舶工程學(xué)院學(xué)報(bào),1993,14(2):49-55.

Su Chenji, Lu Zhiming, Xia Hong. Fluid veloeity measurement by using the correlation analysis of the temperature noise signals in fluid[J]. Journal of Harbi Shipbuilding Engineering Institute, 1993, 14(2): 49-55. (in Chinesewith English abstract )

[19] 周紅仙,周有平,王毅. 用自相關(guān)法測(cè)量橫向流速[J]. 物理實(shí)驗(yàn),2012,32(5):6-8.

Zhou Hongxian, Zhou Youping, Wang Yi. Measuring transverse velocity by autocorrelation[J]. Physics Experimentation, 2012, 32(5): 6-8. (in Chinese with English abstract)

[20] de Lima J L M P, Abrantes J R C B. Can infrared thermography be used to estimate soil surface microrelief and rill morphology?[J]. Catena, 2014, 113(1): 314-322.

[21] de Lima J L M P, Abrantes J R C B. Using a thermal tracer to estimate overland and rill flow velocities[J]. Earth Surface Processes and Landforms, 2014, 39(10): 1293-1300.

[22] Abrantes J R C B, Moruzzi R B, Alexanres S, et al. Comparison of thermal, salt and dye tracing to estimate shallow flow velocities: Novel triple-tracer approach[J]. Journal of Hydrology, 2018, 557: 362-377.

[23] 鞏稼民,郭濤,孫宇航,等. 一種新型紅外靶標(biāo)的設(shè)計(jì)[J].半導(dǎo)體光電,2015,36(5):804-807.

Gong Jiamin, Guo Tao, Sun Yuhang, et al. Design and application of a new type IR drone[J]. Semiconductor Optoelectronics, 2015, 36(5): 804-807. (in Chinese with English abstract)

[24] 袁仁宗,李建勛. 基于灰度圖像的快速質(zhì)心算法[J]. 中國信息化,2013(8):466-467.

Yuan Renzong, Li Jianxun. Computation of centroid of gray level image based on judgement and arithmetic mean[J].Informatization in China, 2013(8): 466-467.(in Chinese with English abstract)

[25] 胡志亮.直流下有機(jī)玻璃表面粗糙度對(duì)表面電荷影響的研究[D]. 保定:華北電力大學(xué),2017.

Hu Zhiliang. Study on the Influence of Surface Roughness on the Surface Charge of PMMA Under DC Voltage[D]. Baoding: North China Electric Power University, 2017.(in Chinese with English abstract)

[26] Zhuang X H, Wang W,Ma Y Y, et al. Spatial distribution of sheet flow velocity along slope under simulated rainfall conditions[J]. Geoderma,2018, 321: 1-7.

Thermal infrared imaging measurement method for shallow flow velocity

Zhang Yan1,3, Shi Haijing1,2,3※, Guo Minghang1,2, Zhao Jun1,2, Zhan Xiaoyun1,2, Ding Chengqin2

(1.712100;2.712100;3.100049,)

A flow velocity is one of the most important physical parameters to quantify the hydraulic characteristics of water flow. The accurate measurement of shallow flow velocities can greatly contribute to understanding and simulating the sediment transport and soil erosion. In this study, a novel observation system of thermal infrared imaging and computer vision was established to measure the velocities of shallow flow. The observation system consisted of three subsystems, such as the thermal tracer control, image acquisition, and transmission, as well as image calculation. Specifically, the control subsystem of the thermal tracer was mainly responsible for the constant temperature of hot water, including the electric heater, temperature sensor, and hot water pump. The subsystem of image acquisition and transmission consisted of a FLIR ONE 3.1.0 thermal infrared camera, 4 thermal infrared targets, and a wireless router, particularly for the thermal infrared images of thermal tracer migration. The subsystem of image calculation was used to extract the high frame rate from the thermal infrared images, including the data storage and computing matching. As such, the velocity of shallow flow was dynamically monitored using the automatic control of thermal tracer, instantaneous image acquisition, image correction, noise removal, and centroid determination. Furthermore, a series of experiments were conducted to verify the system in the Simulated Runoff Hall of the State Key Laboratory of Soil Erosion and Dryland Farming on the Loess Plateau. The experimental glass tank was in the size of 4.6 × 0.196 × 0.1 m3at a gradient of 15° to the horizontal, together with the different flows (0.2-1.5 L/s). The results showed that the observation system presented a higher accuracy than before, suitable for the dynamic transport of thermal tracer in the different temporal and spatial scales. Specifically, an excellent performance was achieved, where the measurement standard deviation of the system was 0.0201 m/s, the observation accuracy reached 98.33%, the observation time resolution was 1/9 s, and the spatial resolution was up to 2 mm. Moreover, the accuracy of the observation system with the thermal infrared imaging was much higher than that with the traditional tracer techniques (dye and salt tracer). The maximum and minimum relative errors of the observation system were –9.61% and 0.16%, respectively, and the range of the relative errors was within±10%. The maximum relative error of the dye tracer was 75.92%, and the minimum relative error was 12.03%. The maximum and minimum relative errors of salt tracer were –26.67% and 0.11%, respectively, where 52% of the samples presented the relative errors within±10%. Correspondingly, the observation system with thermal infrared imaging was provided a reliable way to measure the shallow flow velocity. By contrast, either the dye tracer or salt tracer cannot calculate the overall situation of water flow, due to the manually recording one-time value during measurement. Fortunately, the observation system with thermal infrared imaging can be widely expected to accurately record the overall situation of the water flow. More importantly, the leading edge points and centroid points of the thermal tracer area can be identified at different time intervals, thereby accurately calculating the leading edge velocity and centroid velocity of the water flow. The images can also be used to trace the development of water flow in various tracer sections at different times. As such, it can be possible to measure the tracer migration velocity along the direction of water flow and the diffusion velocity perpendicular to the direction of water flow in the future. This technique can also be applied to monitor the rainfall erosion and runoff scour in the process of soil erosion.

flow velocity; slope; thermal infrared imaging observation system; shallow flow; thermal infrared target; accuracy evaluation

10.11975/j.issn.1002-6819.2021.21.013

P335;S157

A

1002-6819(2021)-21-0108-08

張艷,史海靜,郭明航,等. 基于熱紅外成像的坡面薄層水流流速測(cè)量方法[J]. 農(nóng)業(yè)工程學(xué)報(bào),2021,37(21):108-115.doi:10.11975/j.issn.1002-6819.2021.21.013 http://www.tcsae.org

Zhang Yan, Shi Haijing, Guo Minghang, et al. Thermal infrared imaging measurement method for shallow flow velocity[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2021, 37(21): 108-115. (in Chinese with English abstract) doi:10.11975/j.issn.1002-6819.2021.21.013 http://www.tcsae.org

2021-05-12

2021-10-10

中國科學(xué)院戰(zhàn)略性先導(dǎo)科技專項(xiàng)(XDA20040202);黃土丘陵區(qū)地形微生境分類評(píng)價(jià)與潛在植物群落分布模擬(XAB2020YN04);國家科技基礎(chǔ)條件平臺(tái)建設(shè)項(xiàng)目(2005DKA32300)

張艷,研究方向?yàn)槠旅嫠鞅O(jiān)測(cè)方法。Email:zhangyan199@mails.ucas.ac.cn

史海靜,博士,副研究員,研究方向?yàn)樗帘3峙c科研信息化。Email:shihaijingcn@ nwafu. edu.cn

猜你喜歡
測(cè)量
測(cè)量重量,測(cè)量長度……
把握四個(gè)“三” 測(cè)量變簡單
滑動(dòng)摩擦力的測(cè)量和計(jì)算
滑動(dòng)摩擦力的測(cè)量與計(jì)算
測(cè)量的樂趣
二十四節(jié)氣簡易測(cè)量
日出日落的觀察與測(cè)量
滑動(dòng)摩擦力的測(cè)量與計(jì)算
測(cè)量
測(cè)量水的多少……
主站蜘蛛池模板: 成年免费在线观看| 欧美成人免费一区在线播放| 亚洲人成网站色7777| 日韩高清在线观看不卡一区二区| 久久中文字幕不卡一二区| 亚洲人成日本在线观看| 亚洲欧美人成人让影院| 国产系列在线| 国产成人无码播放| 国产91高清视频| 欧美成人手机在线观看网址| 亚洲永久精品ww47国产| 国产精品福利导航| 伊人久久大香线蕉aⅴ色| 久久精品电影| 在线播放真实国产乱子伦| 国产成人精品男人的天堂| 久久天天躁夜夜躁狠狠| 色综合热无码热国产| 亚洲av日韩综合一区尤物| 久久婷婷五月综合97色| av色爱 天堂网| 欧美高清三区| 欧美三級片黃色三級片黃色1| 国产精品一区二区在线播放| 国产超碰在线观看| 极品尤物av美乳在线观看| 中文字幕第4页| 久久久久久高潮白浆| 国国产a国产片免费麻豆| 日韩一区二区三免费高清| 天天干天天色综合网| 亚洲熟女中文字幕男人总站| 国内精品免费| 国内a级毛片| 一区二区三区成人| 中文字幕永久视频| 秋霞国产在线| 午夜啪啪福利| 99无码中文字幕视频| 波多野结衣二区| 亚洲区一区| 亚洲第一极品精品无码| 亚洲 日韩 激情 无码 中出| 亚洲视频色图| 在线观看亚洲人成网站| 小13箩利洗澡无码视频免费网站| 色偷偷综合网| 亚洲中文字幕在线观看| 国产白浆一区二区三区视频在线| 五月婷婷综合在线视频| 伊人久久精品亚洲午夜| 亚洲欧美国产五月天综合| 57pao国产成视频免费播放| 欧美成人精品一级在线观看| 91免费在线看| 日韩av高清无码一区二区三区| 成人av专区精品无码国产| 亚洲成人www| 亚洲无码37.| 成人综合久久综合| 97国产精品视频自在拍| 老司机午夜精品网站在线观看| 色妞永久免费视频| av在线手机播放| 麻豆a级片| 中文字幕一区二区视频| 欧美三级日韩三级| 国产99欧美精品久久精品久久| 久久国产V一级毛多内射| 99热这里只有精品在线播放| av天堂最新版在线| 亚洲视频无码| 一级毛片免费高清视频| 国产女人18水真多毛片18精品| 91精品国产情侣高潮露脸| 欧美天堂在线| 污网站在线观看视频| 毛片a级毛片免费观看免下载| 永久天堂网Av| 成人免费午间影院在线观看| 精品国产一区91在线|