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

地表細小可燃物含水率實測值與自動測量儀器值的比較和校正

2023-02-11 11:19:32范佳樂郭妍于宏洲孫龍胡海清胡同欣
森林工程 2023年1期
關鍵詞:測量模型

范佳樂,郭妍,于宏洲,孫龍,胡海清,胡同欣

(東北林業大學 林學院, 哈爾濱 150040)

0 引言

近年,森林火災頻發,林火預測預報研究備受關注[1]。森林可燃物含水率直接影響林火發生的難易程度[2],并間接影響著火強度、火蔓延速度等火行為[3]。因此,可燃物含水率被廣泛用作火險等級評估和火行為模型的輸入數據[4-5]。特別是地表細小死可燃物含水率,是指示森林火險等級的重要指標[6]。美國和加拿大等國的森林火險等級系統都將可燃物含水率作為最基礎的輸入數據[7-8]。需要注意的,當可燃物含水率的誤差為±50%(即平均為4%的可燃物含水率的±2%)時,會導致輸出變量(如火災蔓延速率等)中產生高達80%的誤差[9]。可燃物含水率的預測精度要達到平均絕對誤差在1%附近才能滿足火災行為預測和火災風險預測的精度要求[10-11,29,37]??梢姡瑴蚀_連續實時地獲取森林可燃物含水率是提高森林火險預報精度的核心任務[12]。而準確連續地實測可燃物含水率的難度很大,因此需要一種簡單準確的方法對可燃物含水率進行自動連續的監測。

目前常用的可燃物含水率測定方法主要有測量法和預測法2種。測量法主要包括稱重法、探針法和濕度棒法。其中稱重法的結果最為準確,但需要耗費大量的人力、物力和時間,且不能實時地得到可燃物含水率,不能滿足林火預測預報的要求[13];探針法可以快速地測量可燃物含水率,但由于可燃物的結構和水分異質性以及有機質含量高等原因,探針法的準確性和穩定性較差[14];濕度棒法相對于稱重法較為快速,但由于放置可燃物棒的位置沒有與地面接觸,測得的含水率偏低[14],準確度較差[15]。預測法是當前的主要研究內容,指利用可燃物含水率的影響因子通過一系列經驗或機械模型來對可燃物含水率進行預測[16-17]。主要有遙感估測法、氣象要素回歸法、平衡含水率法和過程模型法4類方法[18]。當前最完整和使用最廣泛的方法是基于Catchpole等[19]提出的基于時滯和平衡含水率法,通過Nelson[20]模型,利用氣象因子直接預測可燃物含水率[21],準確性較高,但仍有一定誤差[22-23]。而且由于森林生態系統的異質性,典型可燃物類型的基礎數據薄弱以及可燃物含水率數據觀測困難等原因,導致模型的外推性不是很好[24]。綜上,稱重法的準確性最高,但存在的問題是獲取數據的時間最長、難度最大。本研究中所應用的可燃物含水率自動測量儀器正是基于稱重法以固定的步長測量可燃物含水率并同步測量空氣溫度和空氣相對濕度等氣象數據然后儲存記錄,實現了可燃物含水率的自動實時監測工作,解決了對地表細小可燃物含水率長期高精度自動化測量的難題。但在儀器的前期使用過程中發現其與人工實測數據仍存在一定的誤差,需要分析誤差原因,找尋校正方法[25]。提高儀器的準確性,使其滿足林火預測預報對可燃物含水率的精度要求。

鑒于此,本研究以哈爾濱市城市林業示范基地內的蒙古櫟(Quercusmongolica)林、樟子松(Pinussylvestrisvar.mongolica)林、蒙古櫟-樟子松混交林為對象,通過儀器測量和人工實測2種方法獲得含水率數據和氣象因子數據,將二者進行比較,檢驗儀器測量含水率數據的準確性,結合環境因子和儀器的工作原理分析誤差原因,找尋校正方法。本研究結果可以為校正儀器監測數據提供理論依據,為今后進行快速和準確的火險預測預報提供重要的技術支撐。

1 研究區域概況

研究地點位于東北林業大學校內的哈爾濱市城市林業示范基地(126 °37 ′15 ″E,45 °43 ′11 ″N)。地形平緩起伏較小,海拔136~140 m。屬中溫帶大陸性季風氣候,四季分明,1月平均氣溫-15~-30 ℃,7月份平均氣溫19~23 ℃,年平均氣溫4.2 ℃。年均降水量569.1 mm,集中在6—9月份。土壤類型為地帶性黑鈣土。林分主要為人工林,以小面積純林和塊狀混交林為主。主要樹種有樟子松、胡桃楸(Juglansmandshurica)、蒙古櫟、紅松(Pinuskoraiensis)等。

2 材料和方法

2.1 儀器介紹

可燃物含水率自動測量儀器是一套基于稱重法能夠在野外長期自動連續測量地表細小可燃物含水率的儀器??傮w結構由采集主機、拉力傳感器、電動推桿、網兜、大氣溫濕度傳感器、風速傳感器、光照傳感器、雨量傳感器(選配)、安卓手機軟件、蓄電池、太陽能電池板、測量支架、北斗短報文模塊(選配)和北斗遠程數據接收服務器(選配)等組成,組裝完成后的主體大小約為55 cm×70 cm×120 cm,結構組成如圖1所示。工作原理為將可燃物裝于尼龍網兜內,下部與土壤表面接觸、上部無覆蓋物,保證其與周圍環境的水汽交換,通過高精度的拉力傳感器周期性自動提升樣品,以實現樣品濕重的連續測量,配合前期單獨測得的樣品干重,得到實時的樣品含水率;并在距地面1.2 m處同步地自動采集大氣溫度、大氣相對濕度、風速和光照強度等氣象數據。實現含水率數據和氣象數據的自動實時監測工作。采樣周期可設定為0.05~24 h任意時長。可燃物濕重的最大稱重量程為1 000 g。儀器測得結果中可燃物濕重的精度為± 0.1 g、溫度的精度為± 0.5 ℃、相對濕度的精度為± 0.1%、風速的精度為± 0.1 m/s、光照強度的精度為± 0.1 Lx、降水量精度為± 0.1 mm。

圖1 儀器結構示意圖Fig.1 Schematic diagram of the instrument structure

2.2 樣地設置

本研究于2020年11月1日至2020年11月16日及2021年5月1日至2021年5月21日進行樣地調查和地表可燃物含水率數據監測。在所選取的3種林分(蒙古櫟林、蒙古櫟-樟子松混交林、樟子松林)中分別設置1塊20 m×20 m的標準樣地。在每塊樣地中隨機設置1個樣點,布設森林可燃物含水率自動測量儀器,同時在每臺儀器周圍隨機選取3個樣點人工實測含水率。具體樣地信息見表1。

表1 樣地基本信息

2.3 數據來源

由于每個林分內地勢平緩,坡度變化不大,所以選擇的樣點均為平地。采集樣地當年地表凋落物作為實驗樣品,帶回實驗室,在105 ℃的條件下烘干24 h至恒質量,并記錄為干質量。同時測量尼龍網兜(100目)和尼龍網袋(20目)在正常室內條件下的濕質量后,烘至恒質量,隨后稱質量,記錄為網兜和網袋的干質量,發現變化在0.02 g以內,所以后續測量時不考慮網兜和網袋吸水。將樣品分別裝于尼龍網兜和網袋內,然后將裝有樣品的尼龍網兜與設備的稱重裝置相連接,水平放置于地面之上,下部與土壤相鄰接觸,上部無覆蓋物,使之能夠與周圍環境進行水汽交換。以30 min為步長通過拉力傳感器提升尼龍網兜,獲取樣品濕重和網兜的總重量并記錄,由于不考慮網兜的吸水,所以樣品濕重即為測得的總重量減去之前所測得的網兜干重。隨后將尼龍網袋隨機放置在儀器周圍,要求與上述相同,每個樣地設置3組重復,每天8:00—16:00每隔3 h在樣地人工實測樣品濕重+網袋的重量,并進行記錄。以3個樣點的算數平均值代表可燃物含水率數據。

氣象數據是通過儀器與重量同步(30 min為步長)測量獲得的,包括空氣溫度(T)、空氣相對濕度(H)和風速(W)等氣象參數。由于網兜在稱重過程中,會因過大的風速左右擺動,導致所得數據的誤差偏大,所以將儀器測量時的風速閾值設置為2 m/s,當風速大于閾值時不進行稱重。

2.4 地表可燃物含水率計算

地表可燃物含水率根據公式(1)進行計算

(1)

式中:M為地表可燃物的含水率,%;WH為地表可燃物的濕質量,g;WD為地表可燃物的干質量,g。

2.5 數據處理

將2020年所得的含水率數據與氣象數據進行基本的統計分析,并分析儀器含水率與人工實測含水率的相關性,繪制對比折線圖??梢缘玫蕉唠m然差異顯著,但都在0.01級別極顯著相關。然后,對儀器測量和人工實測含水率的差值同儀器測量含水率Wi以及氣象數據進行Pearson相關性分析和Spearman相關性分析。由于可燃物含水率對氣象因子的響應具有一定的滯后性[26],所以選擇了當期氣象因子以及前1、2、3、4、5 h的氣象因子,分別以下角標0、-1、-2、-3、-4、-5表示,如前2 h的溫度表示為T-2;因為風速閾值的設定導致研究期間大部分風速小于1m/s,且變化幅度較小,對含水率的影響偏低,為了使其對含水率的影響更加顯著,將風作為1個二進制變量進行研究[4]。隨后采用向前逐步回歸的方式對具有顯著相關關系的變量進行篩選,選出構建模型的變量,建立線性回歸模型;并且,考慮到差值與所選變量之間并不一定呈線性關系,還選擇非線性回歸模型中的二次多項式、三次多項式、指數函數和對數函數來建立校正模型[27-28],并計算模型的平均絕對誤差(MAE)和平均相對誤差(MRE),比較模型精度,選擇每種林型中最優非線性模型進行下一步研究。分別將通過2種方法所得校正后的儀器含水率作為縱坐標,人工實測含水率為橫坐標繪制1∶1誤差散點圖。隨后對校正后的儀器含水率與人工實測含水率進行配對樣本t檢驗,判斷二者是否存在顯著差異[24],并計算二者的平均絕對誤差和平均相對誤差,評價模型的校正效果。最后帶入2021年實驗數據,檢驗線性和非線性校正模型的外推性,確定最終的校正方法。

本研究所有的數據處理和繪圖均采用R 4.1.1和SPSS 25.0軟件完成。

3 結果與分析

3.1 可燃物含水率和氣象因子動態變化

研究期間氣象要素的基本統計分析情況見表2??梢钥闯?,在研究期間溫度的最大值為12.400 ℃,最小值為-6.900 ℃,平均值為4.344 ℃,通過75%分位數可以得出研究期間的溫度大部分都低于7.100 ℃。相對濕度的平均值為60.829%,通過其25%分位數可以得出研究期間相對濕度較大,大部分相對濕度高于48.300%。風速的平均值為0.218 m/s,中值為0.100 m/s,其75%分位數為0.400 m/s,標準差為0.319,其在研究期間的變化較小。

表2 氣象要素統計表

圖2給出了研究期間人工實測含水率數據和儀器測量含水率數據的對比圖以及相關系數。其中蒙古櫟林、蒙古櫟與樟子松混交林和樟子松林的人工實測含水率的平均值分別為17.225%、13.065%、12.881%;儀器測量含水率的平均值分別為28.622%、17.742%、24.298%,人工測量值與儀器測量值之間存在明顯的差異,但是從圖2中可以看出,實測數據和儀器數據的變化趨勢基本一致,Pearson相關性分析和Spearman相關性分析結果都在0.01級別極顯著相關。但是3種林型的儀器數據都明顯大于人工實測含水率數據,二者存在顯著差異。這種誤差都屬于同一方向并且不具有相互抵償性,增大樣本量并不能消除誤差。因此認為儀器測量所得數據與人工實測含水率數據之間的誤差屬于系統誤差,需要對一種方法數據進行校正,才可使準確率提高,消除顯著差異。

圖2 人工實測含水率數據與儀器測量含水率數據對比圖Fig.2 The comparison between the data of the moisture content measured manually and by the instrument

3.2 氣象因子的相關性分析

3種林型中,儀器測量含水率和人工實測含水率的差值與儀器含水率數據和氣象數據分別進行Pearson和Spearman相關性分析,結果如圖3所示。由圖3可以看出,2種方法所得結果中儀器測量含水率與3種林型的差值都具有極顯著的正相關關系。此外,Pearson相關性分析的結果中,蒙古櫟林還與1 h前的濕度具有顯著正相關關系;樟子松林與當期溫度具有極顯著負相關關系,與當期風速和2 h前的風速具有顯著正相關關系,且相關性隨距離采樣時間的減少而增加,與3 h前的溫度呈顯著負相關。Spearman相關性分析結果中,除儀器測量含水率外,僅樟子松林中的當期風速還呈顯著正相關關系。

3.3 儀器校正模型

3.3.1 線性回歸模型

表3給出了儀器測量含水率Y與人工實測含水率的差值通過逐步回歸所得的線性擬合模型。由表3可知,3種林型的模型都包括儀器含水率數據,另外蒙古櫟林模型還包括前1 h空氣相對濕度,樟子松林模型還將前3 h的溫度選入模型。樟子松林的模型擬合效果最好,其R2為0.892,大于其余2種林型的模型;MAE和MRE分別為0.012和0.098,均小于其余2種林型的模型。擬合出的線性回歸模型的R2變化范圍為0.509~0.892;MAE的變化范圍為0.012~0.020;MRE的變化范圍為0.098~0.224。

a.蒙古櫟林;b.蒙古櫟-樟子松混交林;c.樟子松林;**,在0.01級別相關性顯著; *,在0.05級別相關性顯著。

表3 儀器測量含水率與人工實測含水率差值的線性擬合模型Tab.3 The linear fitting model of the difference between the instrumental and manual measured moisture content

3.3.2 非線性回歸模型

非線性回歸模型(二次多項式、三次多項式、指數函數和對數函數)的建模結果見表4。由表4可知,在蒙古櫟林中,三次多項式模型的R2均大于其余模型,MAE和MRE 均小于其余模型。這表明對于蒙古櫟林,三次多項式模型為最優非線性模型。在蒙古櫟-樟子松混交林中,三次多項式模型的R2與二次多項式模型相等都大于其余模型,MAE小于其余模型。這表明,對體蒙古櫟-樟子松混交林,三次多項式模型為最優非線性模型。在樟子松林中,指數函數模型的R2均大于其余模型,MAE和MRE雖小于二次多項式模型和三次多項式模型,但差別較小并不顯著。所以,對于樟子松林,指數函數模型為最優模型。綜上,選擇各林型中的最優模型進行下一步的研究,下文中提及的非線性模型均為各林型中的最優模型。

3種林型中,樟子松林的模型擬合效果最好,其R2為0.942,大于其余2種林型的最優模型。MAE和MRE分別為0.007和0.063,均小于其余2種林型的最優模型。3種林型的最優非線性回歸模型的R2變化范圍為0.565~0.942;MAE的變化范圍為0.009~0.015;MRE的變化范圍為0.077~0.486。

表4 儀器測量含水率與人工實測含水率差值的非線性擬合模型Tab.4 The non-linear fitting model of the difference between the instrumental and manual measured moisture content

3.3.3 模型誤差分析

通過線性和非線性模型所擬合出的差值對儀器測量含水率進行校正,然后利用校正后的儀器數據與人工實測含水率數據繪制1∶1散點圖,同時計算二者的平均絕對誤差,如圖4所示。由圖4可以看出,蒙古櫟林中分別通過線性模型和非線性模型進行校正后的儀器值與實測值的MAE分別為0.016和0.014,混交林中分別為0.017和0.015,樟子松林中分別為0.012和0.009,變化范圍為0.009~0.017。無論是線性回歸模型還是非線性回歸模型其模型擬合的準確度由大到小都是樟子松林、蒙古櫟林、蒙古櫟-樟子松混交林。非線性回歸模型的擬合效果在每個林型都要優于線性回歸模型,其擬合線更接近于1∶1線。而且根據3.3.1和3.3.2小節可以得出,2種方法所構建的回歸模型中非線性模型的R2、MAE和MRE都要稍稍優于線性模型,但二者的差異并不明顯。

圖4 校正后的儀器值與實測值1∶1散點圖Fig.4 1∶1 scatter plot of the corrected instrumental value and the measured value

3.4 校正后儀器測量含水率數據的檢驗

表5給出了各林型校正后儀器數據與實測數據進行配對樣本t檢驗的結果。由表5可以看出,所有校正后的儀器數據都與實測數據具有顯著的正相關性。3種林型在通過線性回歸模型得到的校正后儀器數據與實測數據的相關性都要小于非線性回歸模型的。2種方法所有樣本檢驗的結果都為雙尾顯著性水平的P大于0.05。這就說明這些校正后的儀器數據與實測數據之間沒有顯著差異,通過校正,儀器數據的準確性提高。

3.5 模型外推性研究

將2021年所得數據代入到上述模型中,然后進行配對樣本t檢驗,判斷模型的外推性,其結果見表6。由表6可以看出,3種林型通過線性模型進行校正后,儀器數據與實測數據都具有顯著的正相關關系,t檢驗的結果都為P>0.05,表明二者之間沒有顯著差異,外推性良好;3種林型通過非線性模型進行校正后,蒙古櫟林和樟子松林的校正后儀器數據與實測數據具有顯著正相關關系,混交林為顯著負相關關系。t檢驗的結果都為P<0.05,具有顯著差異。這說明相比于通過線性模型進行校正后的儀器數據,非線性模型校正后與實測數據的差距較大,校正效果和外推性都較線性模型差。并且儀器數據與實測數據在理論上應具有相同的變化趨勢,二者呈顯著正相關。而通過非線性模型校正后混交林的校正后儀器數據與實測數據呈顯著負相關??梢缘贸?,非線性模型的外推性并不理想,校正后仍然具有顯著差異。所以,選擇線性回歸模型作為儀器測量含水率數據的校正方法。

表5 校正后儀器數據與實測數據樣本檢驗Tab.5 Sample inspection of instrumental data and measured data after calibration

表6 模型外推性檢驗Tab.6 Model extrapolation examine

4 結論與討論

本研究表明,儀器測量的可燃物含水率數據顯著高于人工實測的數據,二者具有顯著差異。但是,他們的變化趨勢相似且都在0.01級別顯著相關。通過對誤差的性質和造成原因進行分析發現儀器測量的含水率數據與人工實測含水率之間的誤差屬于系統誤差,需要進行校正以提高儀器數據的準確性,滿足林火預測預報的需要。

造成系統誤差的原因是網兜上所附著的水分和土壤顆粒等雜質。在野外條件下,林內環境復雜、濕度較大且降雨頻繁。裝可燃物用的尼龍網兜雖然自身的吸水性較差,正常室內環境下的質量與干質量的差值都在0.02 g以內,但網兜的目數較大、網格密集。網格上以及網兜邊緣的縫隙間會附著一些水分;網兜底部與土壤表面接觸,底部所附著的水分會吸附一些土壤顆粒等雜質,使其黏著在網兜底部,導致質量增大。儀器的工作原理是通過拉力傳感器推桿以固定的步長將裝托可燃物的網兜提起并進行稱質量,稱質量時無法清除網兜上水分等雜質。而計算過程中僅將儀器獲取的總質量減去網兜的干質量作為樣品的濕質量,所以樣品濕質量中仍包含網兜上所附著的水分和土壤顆粒等。但在人工實測過程中,為保證其準確性,在稱質量前會將網袋上所附著的水分和雜質盡量抖落。因此儀器測得的含水率數據中誤差主要來自網兜上所附著的水分和土壤顆粒等雜質。此外,尼龍網兜的目數不同其對可燃物與周圍環境之間的水汽交換造成的影響可能也不同,目數越大影響越大。在試驗過程中我們還發現可燃物的腐爛分解可能也會對測得的含水率數據產生影響??扇嘉镌谧匀画h境中會不斷地腐爛分解[29],廖利平等[30]的研究表明杉木(Cunninghamialanceolata)與火力楠(Micheliamacclurei)凋落物的1∶1混合物在一年時間分解了30%左右。在今后的研究中,應考慮這些方面的影響,設計試驗對其進行定性定量的研究。

氣象因子對儀器測量含水率數據與人工實測數據之間的顯著差異具有一定影響,但相對而言關系較小。本研究中,除儀器含水率數據,還分別將前1 h的相對濕度和前3 h的溫度選入校正模型。這與之前的研究結果相似[31-32],溫度和相對濕度是影響可燃物含水率的主要氣象因子,二者在影響可燃物含水率時可以起到協同作用[33],且含水率對氣象因子的響應具有一定的滯后性[26]。風速沒有被選入模型,這是由于本研究中儀器風速閾值的設置,導致風速變化較小,其對含水率的影響較小[34]。對風速數據進行二進制處理使其對含水率的影響更為顯著[4],但結果沒有發生變化。氣象因子雖被選入模型,但通過相關性分析可以看出,氣象因子與儀器數據和實測數據之間差值的關系相對較小,對校正模型的影響也相對較小。這就表明儀器數據與人工實測數據之間顯著差異的主要是由于儀器測量存在系統誤差導致,與氣象因子的關系較小。

本研究中沒有考慮降雨和土壤濕度的影響,因為當可燃物含水率大于35%時就基本不會燃燒[35],且雨后可燃物含水率存在較大的觀測誤差[36]。土壤濕度是可燃物含水率變化的重要因素[37-38]。但Rakhmatulina等[4]研究發現,土壤濕度對1 h可燃物含水率的影響雖然很顯著但是很小,土壤濕度只能解釋6%的可燃物含水率變化。此外,不同高度的溫度和相對濕度會對可燃物含水率的預測結果造成影響,在可燃物表面附近測得溫度和相對濕度的準確度相對較高[10]。而本研究所使用的溫濕度等氣象因子都是在距地面1.2 m處測得,在后續應進一步補充研究。

線性回歸模型具有良好的校正效果和外推性,可以對儀器測量含水率數據進行校正,校正后的數據滿足精度要求。本研究分別構建線性和非線性回歸模型對儀器測量含水率進行校正。所有校正模型都將儀器測量含水率數據選入其中,另外蒙古櫟林模型還包括前1 h相對濕度,樟子松林模型還將前3 h的溫度選入模型。2種校正模型中非線性回歸模型的精度要略優于線性回歸模型,但二者的差異并不顯著。校正后的儀器數據與人工實測數據都能通過配對樣本t檢驗,沒有顯著差異;二者的MAE均在1%附近,校正效果較好,滿足林火預測預報的精度要求[10-11]。但是,在外推性方面,只有線性回歸模型具有較強的外推能力,可以滿足校正后儀器測量含水率數據與實測數據沒有顯著差異;而非線性回歸模型的外推性并不理想,校正后仍然具有顯著差異。這可能是因為在構建非線性模型的過程中,由于數據量較少或模型的復雜度過高等問題出現了過擬合現象[39],導致非線性模型的外推性只能適用于建模數據而外推性表現很差,在后續的校正研究中要針對此問題進一步分析研究。

本研究表明,可燃物含水率自動測量儀器存在一定的誤差,但可以通過線性回歸模型進行校正,校正后的儀器含水率數據滿足林火預測預報的精度要求。但是,本研究沒有考慮到可燃物的腐爛分解造成的影響;并且只在哈爾濱地區的典型林分內進行試驗,沒有考慮地形地勢以及降雨的影響,具有一定的局限性。今后還應針對上述問題,對儀器的校正進行更進一步的研究。本研究結果可以為校正含水率自動測量儀器所得數據提供理論依據,提高了儀器測量含水率數據的準確性,具有良好的推廣性和實際應用價值。

猜你喜歡
測量模型
一半模型
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
把握四個“三” 測量變簡單
滑動摩擦力的測量和計算
滑動摩擦力的測量與計算
測量的樂趣
3D打印中的模型分割與打包
測量
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 国产真实自在自线免费精品| 国产精品成人不卡在线观看| 国产欧美精品一区二区| 98超碰在线观看| 国产在线专区| 在线国产欧美| 亚洲成人免费看| 久久精品视频一| 乱系列中文字幕在线视频| 97精品国产高清久久久久蜜芽| 婷婷综合缴情亚洲五月伊| 亚洲天堂精品在线| 99视频只有精品| 无码电影在线观看| 一级毛片免费播放视频| 色天天综合| 亚洲熟女中文字幕男人总站| 婷婷综合亚洲| 情侣午夜国产在线一区无码| 少妇人妻无码首页| 久久免费视频6| 国产网站免费观看| 成人国产小视频| 一级成人a做片免费| 欧美色图久久| 毛片久久网站小视频| 激情综合网址| 无码一区中文字幕| 国产性生交xxxxx免费| 国产办公室秘书无码精品| 午夜三级在线| 亚洲精品第一在线观看视频| 欧美日韩理论| 国产色网站| 国产乱码精品一区二区三区中文 | 波多野结衣第一页| 亚洲综合色吧| 一本色道久久88| 人人妻人人澡人人爽欧美一区| 国产欧美日韩资源在线观看| 日韩精品一区二区三区大桥未久| 国产精品尹人在线观看| 亚洲日韩在线满18点击进入| 四虎影视8848永久精品| 中文字幕在线欧美| 国产精品露脸视频| 丰满人妻被猛烈进入无码| 日本免费新一区视频| 国内黄色精品| 国产成人毛片| 国产精品天干天干在线观看| 亚洲第一色视频| 国产成人一区| 色婷婷电影网| 国产日韩久久久久无码精品| 国产在线一区视频| 最新国产你懂的在线网址| 激情综合网激情综合| 2024av在线无码中文最新| 一级不卡毛片| 久久精品这里只有国产中文精品| 狠狠ⅴ日韩v欧美v天堂| 免费高清自慰一区二区三区| www.国产福利| 久久综合色视频| 在线视频精品一区| 91精品视频在线播放| 99九九成人免费视频精品| 亚洲欧美一区在线| 精品欧美一区二区三区久久久| 中文字幕亚洲精品2页| 国产免费久久精品99re丫丫一| 欧美精品亚洲精品日韩专区| 成人av专区精品无码国产| 国产精品尤物铁牛tv| 久久这里只有精品2| 亚洲色无码专线精品观看| 亚欧成人无码AV在线播放| 久久国产精品电影| 国产成人1024精品下载| 亚洲欧洲日韩综合色天使| 亚洲不卡网|