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

基于OLI數據的信陽市境內淮河流域水質遙感反演

2021-01-16 02:49:09張宏建周健皇甫款
人民長江 2021年12期
關鍵詞:水質模型

張宏建 周健 皇甫款

摘要:運用遙感影像進行連續、大范圍的水質監測是目前研究水環境問題的熱點之一。為探索更多可以利用遙感監測的水質指標,以快速獲取水質指標時空分布情況,為水環境治理提供數據支持,以河南省信陽市境內淮河流域為研究對象,根據實測水文監測數據建立了OLI數據與水質指標(WQI)之間的統計模型,包括溶解氧(DO)、氟化物、氨氮(NH3-N)、總磷(TP)以及pH等,以探尋更多能夠利用OLI數據進行監測的水質參數。結果表明:① OLI數據可以用于監測氨氮、總磷、pH值、氟化物、溶解氧的空間分布;② 研究區內整體水質偏堿性,氟化物含量較低,不存在過度氟污染;③ 干支流交匯處以及城區河段氨氮、總磷污染較重,溶解氧含量低,需加強水質管理。

關 鍵 詞:OLI數據; 水質指標; Landsat8影像; 水質反演; 遙感反演; 信陽境內淮河流域

中圖法分類號: X832;X87

文獻標志碼: A

DOI:10.16232/j.cnki.1001-4179.2021.12.008

0 引 言

“綠水青山就是金山銀山”是習近平總書記關于生態文明建設最著名的科學論斷之一,揭示了保護環境就是保護生產力,改善環境就是發展生產力的客觀規律。據不完全統計,中國淮河、海河、松花江、黃河、遼河水系的水質都受到了不同程度的污染[1]。受到污染的河流由于匯流作用間接影響了湖泊、海洋等水體的水質情況[2-4]。促進水生態文明建設,必須加強內陸水體水質監測,改善內陸水環境問題。傳統的水質監測方法即為布設監測斷面對水體環境中存在的污染因素或者污染物進行監測,這種監測方式需投入大量人力、物力、財力[5-7]。除此之外中國水域面積較大,相對來講目前全國所設置的固定監測斷面或水文站的數量十分有限,這也使得傳統的監測方法的局限性更加明顯。

衛星遙感影像覆蓋范圍廣、成像速度快、信息量豐富且性價比高,能有效彌補傳統的水質監測方法費時、費力、難以同步反映整個區域的缺點,非常適合用于水體污染物濃度的監測。隨著遙感空間技術的發展,利用遙感監測水質信息的工作也逐漸展開。利用遙感對水質信息進行監測的技術逐漸從定性監測發展到定量監測[8-9]。國內外應用遙感監測水質指標中大多研究區在水面較為寬廣的湖泊或海岸[10-12]。水質監測中常用的遙感數據主要有HRV數據、Landsat系列數據、EO-1衛星數據、AVHRR數據、環境1號數據等。考慮到時間分辨率、空間分辨率、數據的易獲得性及光譜分辨率等因素,目前內陸水體定量遙感中比較常用的是免費、易獲得的中高分辨率影像數據[13]。Landsat8影像比之前的Landsat系列信噪比高,對于水質監測部門來說,是目前非常好用的數據源。

水污染受多種因素影響[12],常規監測的水質指標(WQI,Water Quality Index)有葉綠素a、濁度、總懸浮物濃度、溶解氧、氨氮、五日生化需氧量、高錳酸鹽指數和化學需氧量等[14-16]。本文研究了Landsat8遙感衛星的OLI(Operational Land Imager)數據與水質指標(WQI)之間的統計模型,包括溶解氧(DO)、氟化物、氨氮(NH3-N)、總磷(TP)以及pH值,重點在于分析Landsat8影像各波段或波段組合與各類指標的相關性,以及各水質指標之間的相關性,旨在探索更多適合利用遙感進行監測的水質指標,豐富水質反演方面的研究。在此基礎之上,尋找利用Landsat8影像監測WQI的最優波段組合,為Landsat8數據在水質遙感中的應用(尤其是不易直接進入或者采樣和化學分析成本較高的偏遠地區的應用)提供數據支撐。

1 材料與方法

1.1 Landsat8遙感數據介紹

本研究從地理空間數據云上面下載了2016年3月1日的 Landsat8衛星影像,云量1.21%,影像質量較高。Landsat8衛星攜帶有陸地成像儀(OLI)和熱紅外傳感器[14](TIRS),將波段5,6,4進行假彩色合成,水體邊界清晰,能有效區分陸地和水體,深淺的橙色和綠色是陸地,深淺藍色是水體。表1列出了Landsat8衛星的波段信息。本文主要使用了0.450~2.300 μm的波段數據。

1.2 Landsat8數據預處理

首先利用envi5.3對獲取的遙感影像進行預處理,包括輻射定標和大氣校正。輻射定標可消除傳感器本身的影響,它將記錄的原始DN值(像元亮度值)轉換為大氣外層表面反射率L(或稱為輻射亮度值)。

L=gain×DN+offset(1)

gain=Lmaxλ-LminλfullDNrange=Lmaxλ-LminλDNmaxλ-DNminλ(2)

offset=Lminλ(3)

式中:每個含有λ角標的參數表示不同波段影像取值不同,具體參數可以從衛星影像的頭文件中得到;L是某個波段光譜輻射亮度;gain為增量校正系數;offset為校正偏差量;DN是圖像灰度值;DNmax和DNmin為遙感器最大和最小灰度值;Lmax和Lmin分別為最大和最小灰度值所相應的輻射亮度。

大氣校正是使用多光譜遙感數據進行地表參數定量分析的前提,主要用于減少或消除氣溶膠和大氣分子的吸收和散射對地物反射率的影響。本研究主要是用FLAASH大氣校正方法對原始數據進行處理,處理結果如圖1所示,預處理過程均在ENVI5.3環境下進行。為了減少岸邊混合像元造成的偶然誤差,在選取各波段圖像上取樣點的DN值的時候采取3×3的矩陣,選取掩模內灰度值最小的像元(與其他地物相比,水體的反射率相對較低)進行取樣。

1.3 相關分析方法

相關分析常用的5種方法為圖表相關分析、協方差及協方差矩陣、相關系數、一元或多元回歸、信息熵和互信息。由于本文需要提煉模型進行預測,因此選用回歸的方式進行相關分析。回歸分析過程中主要采用單波段模型和波段組合模型。單波段模型即在單波段影像的水體反射率與各水質指標之間建立相關關系。波段組合模型中組合方式主要有差值模型、比值模型及一些水體指數模型[17-19],如公式(4)~(11)所示。

EWI=B3-B5-B6B3+B5+B6(4)

RNDWVI=3×B2+1.5×B3-B5-

2×B6+B7(5)

RNDWI=B6-B4/B6+B4(6)

AWEInsh=4×B3-B6-(0.25×B5+

2.75×B7)(7)

AWEIsh=B2+2.5×B3-1.5×

B5+B6-0.25×B7(8)

IWS=4×B6-B2B6-2×B6B2(9)

NDWI=(B3-B5)/(B3+B5)(10)

NWI=NDWI-RNDWVI+AWEInsh+AWEInshNDWI+RNDWVI+AWEInsh+AWEInsh(11)

式中:B2,B3,B5,B6,B7分別為Landsat8影像的Blue,Green,NIR,SWIR1,SWIR2的表觀反射率;

EWI為增強型水體指數,RNDWVI為修訂型歸一化植被指數,RNDWI為修訂型歸一化水體指數,NDWI為歸一化水體指數,NWI為新型水體指數,IWS、AWEInsh和AWEIsh為自動化水體提取指數。

水質反演模型有經驗統計模型和機器學習模型。經驗統計模型簡單有效、準確性差些,但對樣本數量無較高要求。人工神經網絡等機器學習方法的準確性高,但需要大量樣本數據進行預訓練,且無法提供具體的定量表達式。因此,筆者所構建的反演模型主要考慮了5種經驗統計模型,如表2所列。

1.4 遙感水質監測機理

一般來說清潔水體有固定的電磁強度,而一些污染指標如藻類、富營養化物質、有機污染物等的電磁強度與水體不同,一旦水體被污染,水體的表面特性就會發生變化,進而改變水表面對電磁波的反射強度,而這改變了的強度具體表現為遙感圖像上不同的地物水體反射率值。理論上來講,水體表面反射率與水體污染指標之間可能存在某種函數關系。

Rλ=F(x)(12)

式中:Rλ是λ波段的水體反射率,x是某污染指標的濃度。

遙感水質監測的方法主要有直接法和間接法。直接法通常是根據光學活性水質參數與遙感數據之間的內在聯系直接建立反演模型進行水質參數的反演。間接法通常適用于非光學活性水質參數,基于遙感反演得到的光學活性水質參數,利用非光學活性水質參數與光學活性水質參數之間強相關性,建立模型,間接進行非光學活性水質參數的反演。由于水體類型差異、季節性差異和地域差異,某些非光學活性水質參數與光學活性水質參數的相關性不固定,這時一般根據經驗關系利用直接法構建統計模型對非光學活性水質參數進行反演。

2 實驗與分析

2.1 實驗區概況

信陽市境內淮河流域屬于亞熱帶季風氣候,降雨多集中在5~10月份,年均降水量1 300 mm左右,自北向南遞增,氣候溫暖濕潤,相對濕度年均77%。淮河自西向東橫貫信陽市,境內河長363.5 km,地面落差178 m,全境98.2%面積屬于淮河流域。信陽市境內淮河流域的主要支流多位于淮河南岸,較大的支流有竹竿河、潢河、浉河、白露河等。市內建有5座大型水庫,13座中型水庫,866座小型水庫。本實驗選取了3座大型水庫(南灣水庫、石山口水庫、潑河水庫)、淮河干流及其3個主要支流(竹竿河、潢河、浉河)作為實驗區域,研究區及監測斷面位置示意如圖2所示。

2.2 水質數據

從信陽市水文局獲取到2016年淮河流域信陽市水功能區各監測斷面的實地監測資料。數據顯示整個研究區遠未達到水功能區水質控制目標(Ⅱ類水或Ⅲ類水),尤其是枯水期期間,因此實時高效地掌握該區域的水體污染情況很有必要。這里以枯水期期間2016年3月1日的監測數據為例,選取主要污染指標DO、氟化物、NH3-N、TP以及pH為研究對象。本文選取20個監測斷面的水質數據進行實驗,隨機選取16個監測斷面作為訓練集,其余斷面數據作為測試集。

2.3 光譜分析

將各水質指標按不同濃度分為若干組進行光譜特征分析,結果如圖3所示。隨著TP、NH3-N、氟化物濃度的升高,遙感反射率逐漸升高。隨著DO濃度的增加,遙感反射率逐漸減小,且這一特征在B4波段表現尤為明顯,B5、B6和B7波段表現不太明顯。以上現象均表明水體污染程度(水質指標的含量)可通過光譜反射率進行區分,初步認為利用遙感進行這幾類水質指標的監測是可行的。

2.4 水質敏感波段選擇

本文采用直接法對淮河流域信陽段水域進行監測。主要采用表2中的5種數學模型進行回歸分析,通過回歸分析的方法來確定水質參數的敏感波段,以影像數據作為自變量,水質數據作為因變量進行分析。最后將相關系數最高的數學模型作為該波段或波段組合與水質參數之間的相關系數,繪制相關矩陣(見圖4)。自變量的顯著性是根據曲線估計結果中各個自變量系數后面的Sig值判斷:當Sig值大于0.05時,自變量的顯著性較低,通常認為該自變量無法反演因變量的值并將其舍去;當Sig值小于0.05時,在95%的顯著水平下顯著;當Sig值小于0.01時,在99%的顯著水平下顯著。回歸模型的決策系數(R2)越接近1,兩者的擬合度就越高,決策系數越接近零,兩者的擬合度就越低。這里規定相關系數在0~0.2之間為極弱相關或不相關;相關系數在0.2~0.4之間為弱相關;相關系數在0.4~0.6之間為中等強度相關;相關系數在0.6~0.8之間為強相關;相關系數在0.8以上為極強相關。

為了提高水質指標與影像數據之間的擬合度,將影像數據進行波段組合,重新進行相關分析。表3展示了決策系數在0.6以上且Sig值小于0.05的波段組合。實驗發現NH3-N與B4/B2、(B2-B4)/(B2+B4)之間的相關系數均可達0.85;TP與(B3-B4)/(B3+B4)之間的決策系數最高,兩者擬合度在0.93以上;DO與B4/B6之間的決策系數最高,兩者之間的擬合度在0.79以上;氟化物與(B3-B5)/(B3+B5)之間的決策系數最高,兩者之間的擬合度在0.88以上;與pH擬合效果比較好的波段組合是RNDWI,相關系數達0.70。實驗表明利用波段組合的方法可明顯提高部分水質指標與影像數據之間的擬合度。

對各水質指標進行相關分析,相關系數矩陣如圖5所示。DO與TP、NH3-N之間均具有較強的相關性。理論上來講,聚磷菌是好氧細菌,TP含量高的時候,水底的聚磷菌迅速繁殖,繁殖過程消耗水中氧氣,從而降低DO含量。NH3-N是水體中的主要耗氧污染物之一,一定濃度下NH3-N含量過高會促進藻類等植物的生長,也會消耗水中的DO。因此,NH3-N和TP均與DO之間具有一定的負相關性,這與本實驗分析結果相同。pH、NH3-N、TP兩兩之間均具有較強的相關性,其中NH3-N和TP之間存在較強的正相關性,而pH與NH3-N(或TP)之間具有明顯的負相關關系,在進行數據審核時,可利用其相關關系判斷數據的合理性。

3 水質指標的空間分布

3.1 反演模型的構建

在反演模型的構建過程中,一定情況下可能會出現過擬合現象,導致訓練數據的擬合效果很好,而測試數據的擬合效果較差,為了避免這種情況,對表3中所列出的波段或波段組合分別進行反演模型的構建,并利用測試數據進行測試,從中選出測試結果精度較高的模型作為最優的反演模型對研究區進行反演。通過多次實驗,本文確定了適用于該研究法的反演模型,如公式(13)~(17)所示。從最終確定的反演模型特征中可以看出,該研究區最適宜的數學模型為三次擬合模型和二次擬合模型。從各水質參數實測值和估計值對比結果(見圖6)可以看出:本文所構建的模型能夠對水質參數實現一定精度的估計。最后通過對測試集的驗證,得出各模型平均絕對誤差分別為0.42 mg/L,0.14(無量綱),0.87 mg/L,0.04 mg/L,0.02 mg/L,平均相對誤差分別為36.37%,1.77%,10.65%,29.74%,8.90%。

INH3-N=-0.69+20.778×AWEIsh-109.27×

AWEIsh2+186.927×AWEIsh3(13)

IpH=3.134+61.862×B3-180.381×B32(14)

IDO=39.199-49.553×B4B6+

24.849×B4B62-3.936×B4B63(15)

I氟化物=0.259+1.19×NDWI-16.713×

NDWI2+47.438×NDWI3(16)

ITP=0.355-5.894×(B3-B4B3+B4)+19.437×

(B3-B4B3+B42+58.929×(B3-B4B3+B43(17)

3.2 水質指標空間格局特征

以淮河干流及其支流潢河為例,基于以上模型對水質進行反演,各水質參數的時空分布情況如圖7所示。整體上看,研究區內水質良好,大體趨勢為河岸水質略差于河心水質,入河口和流經城區居民點的河段水質略差。具體地,從空間分布上來看,pH在6.018~8.438之間波動,平均值為8.196,水體整體偏堿性。氟化物(以F計)最大值為0.864 mg/L,平均濃度為0.249 mg/L,均在地表水環境控制標準Ⅰ、Ⅱ、Ⅲ類水限值以內。NH3-N濃度在0.016~4.087 mg/L之間,均值為0.807 mg/L。

部分區域水體中的TP含量在0.400 mg/L以上,已經超過了Ⅴ類水的標準限值。NH3-N和TP在潢河入淮河口和流經居民區的河段水質較差,相對應地該區域內溶解氧含量較低,屬于水質污染比較嚴重的區域。

4 結 語

Landsat8因其時空分辨率相對較高且免費易獲取,在水質監測方面具有巨大的應用前景。通過以上實驗,探索了利用OLI數據監測各水質指標的可行性并分析了反演各水質指標的最優波段組合。豐富了狹長的內陸河流水質反演方面的研究內容,實驗結果對于監測站點少、采樣和化學分析總體成本較高的地區具有重要的應用價值。

目前在水質反演過程中仍然存在許多問題。由于一些水體污染的參數(如水體重金屬污染)無法通過遠程遙感確定,目前遙感監測技術無法完全取代傳統的監測方法。此外,對于水質遙感監測來說,季節性差異和區域性差異較大,同一水質指標的反演模型參數也很難統一,仍需根據不同季節、不同區域的實際情況有針對性地構建適合當地的相關模型,為相關機構的監測和水污染防治提供數據支撐。

參考文獻:

[1] 張寶鋒,陳峰,田曉慶,等.2005~2017年中國七大水系水質變化趨勢分析[J].人民長江,2020,51(7):33-39.

[2] 丁夢嬌,丘仲鋒,張海龍,等.基于NPP-VIIRS衛星數據的渤黃海濁度反演算法研究[J].光學學報,2019,39(6):9-17.

[3] 王麗艷,史小紅,孫標,等.基于MODIS數據遙感反演呼倫湖水體COD濃度的研究[J].環境工程,2014,32(12):103-108.

[4] 馬榮華,戴錦芳.結合Landsat ETM與實測光譜估測太湖葉綠素及懸浮物含量[J].湖泊科學,2005,17(2):97-103.

[5] 王佳楠.基于BP神經網絡的高錳酸鹽指數預測研究[D].西安:長安大學,2017.

[6] 岳佳佳,龐博,張艷君,等.基于神經網絡的寬淺型湖泊水質反演研究[J].南水北調與水利科技,2016,14(2):26-31.

[7] 黃妙芬,宋慶君,毛志華,等.應用CDOM光學特性估算水體COD:以遼寧省盤錦市雙臺子河和遼東灣為例[J].海洋學報(中文版),2011,33(3):47-54.

[8] 王皓,趙冬至,王林.水質遙感研究進展[J].海洋環境科學,2012,31(2):285-288.

[9] KOWN Y,BAEK S,LIM Y,et al.Monitoring coastal chlorophyll-a concentrations in coastal areas using machine learning models[J].Water,2018,10(8):1020.

[10] DRNHFER K,OPPELT N.Remote sensing for lake research and monitoring:recent advances[J].Ecological Indicators,2016,64:105-122.

[11] 于祥.渤海表面非光學活性水質參數MODIS遙感定量反演技術研究[D].煙臺:中國科學院煙臺海岸帶研究所,2017.

[12] 彭保發,陳哲夫,李建輝,等.基于GF-1影像的洞庭湖區水體水質遙感監測[J].地理研究,2018,37(9):1683-1691.

[13] 張兵,申茜,李俊生,等.太湖水體3種典型水質參數的高光譜遙感反演[J].湖泊科學,2009,21(2):182-192.

[14] 袁金國,牛錚,王錫平.基于FLAASH的Hyperion高光譜影像大氣校正[J].光譜學與光譜分析,2009,29(5):1181-1185.

[15] 閆大江.基于Landsat系列遙感影像的不同渾濁度內陸地表水體提取研究[D].西安:西北大學,2018.

[16] WANG Y,XIA H,FU J,et al.Water quality change in reservoirs of Shenzhen,China:detection using LANDSAT/TM data[J].Science of the Total Environment,2004,328(1-3):195-206.

[17] FEYISA G L,MEILBY H,FENSHOLT R,et al.,Automated water extraction index:a new technique for surface water mapping using landsat imagery[J].Remote Sensing of Environment,2014,140:23-35.

[18] 王小平,張飛,ABDUWASIT G,等.艾比湖流域地表水水質指標與水體指數關系研究[J].環境科學學報,2017,37(3):900-909.

[19] 萬建鵬,官云蘭,葉素倩,等.基于綜合權重水體指數的水體提取研究:以鄱陽湖為例[J].東華理工大學學報(自然科學版),2015,38(2):206-210.

(編輯:劉 媛)

Analysis on water quality monitoring indicators by remote sensing based on OLI data:

case of Huaihe River Basin in Xinyang City

ZHANG Hongjian1,ZHOU Jian1,HUANGFU Kuan2

(1.Xinyang Water Conservancy Survey and Design Institute,Xinyang 464000,China; 2.School of Water Conservancy Science and Engineering,Zhengzhou University,Zhengzhou 450000,China)

Abstract:

Using satellite images for continuous and large-scale water quality monitoring is currently one of the hotspots in water environmental issue research.The purpose of this study is to explore the water quality indicators that can be monitored well by remote sensing images,so as to quickly obtain the temporal-spatial distribution of water quality indicators and provide data support for water environmental governance.Based on the measured hydrological monitoring data of Huaihe River in Xinyang territory in Henan Province,a statistical model of OLI data and water quality indicators (WQI) was established,including dissolved oxygen (DO),fluoride,ammonia nitrogen (NH3-N),total phosphorus (TP) and pH and so on,to explore all the water quality parameters that can be monitored with OLI data.The results show that the OLI data can be used to monitor the spatial distribution of ammonia nitrogen,total phosphorus,pH,fluoride and dissolved oxygen.The overall water quality in the study area was alkaline and the fluoride content is low,and there was no excessive fluorine pollution.The urban river sections and confluences of the main stream and the tributary rivers were heavily polluted by ammonia nitrogen and total phosphorus,and the dissolved oxygen contents of these two sections were low.The water quality management of the Xinyang City needs to be strengthened.

Key words:

OLI data;water quality indicators;Landsat8 image;water quality inversion;remote sensing;Huaihe River Basin in Xinyang City

猜你喜歡
水質模型
一半模型
水質抽檢豈容造假
環境(2023年5期)2023-06-30 01:20:01
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
一月冬棚養蝦常見水質渾濁,要如何解決?這9大原因及處理方法你要知曉
當代水產(2019年1期)2019-05-16 02:42:04
這條魚供不應求!蝦蟹養殖戶、垂釣者的最愛,不用投喂,還能凈化水質
當代水產(2019年3期)2019-05-14 05:42:48
圖像識別在水質檢測中的應用
電子制作(2018年14期)2018-08-21 01:38:16
3D打印中的模型分割與打包
濟下水庫徑流水質和垂向水質分析及評價
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 狠狠亚洲五月天| 性69交片免费看| 国内熟女少妇一线天| 麻豆精品视频在线原创| 色悠久久久| 久久国产拍爱| 中文字幕一区二区人妻电影| 欧美福利在线观看| 国产成人亚洲欧美激情| 99热这里只有精品5| 国产又色又爽又黄| 性视频一区| 欧美日韩va| 欧美色香蕉| 在线观看欧美国产| 久久久久夜色精品波多野结衣| 中文字幕免费视频| 亚洲区第一页| 久久久久久尹人网香蕉| 婷婷色狠狠干| 亚洲综合狠狠| 国产黄在线观看| 99r在线精品视频在线播放| 亚洲有无码中文网| 亚洲综合第一页| 亚洲中文在线视频| 毛片基地美国正在播放亚洲| 久热99这里只有精品视频6| 午夜国产在线观看| AV老司机AV天堂| 2021天堂在线亚洲精品专区| 日韩毛片免费观看| 日韩福利在线视频| 71pao成人国产永久免费视频| 黄色一级视频欧美| 色呦呦手机在线精品| 97狠狠操| 2021国产精品自拍| 中文字幕亚洲电影| 欧美影院久久| 91福利在线观看视频| 一级黄色网站在线免费看| 亚洲福利网址| 老司国产精品视频| 在线看片中文字幕| 97国产一区二区精品久久呦| 都市激情亚洲综合久久| 亚洲男人天堂网址| 伊伊人成亚洲综合人网7777| 亚洲第一色网站| 1024国产在线| 在线观看国产一区二区三区99| 亚洲无码A视频在线| 国产成人高清在线精品| 亚洲第一色网站| 四虎国产精品永久一区| 国产99视频免费精品是看6| 成人免费网站久久久| 国产精品精品视频| 一区二区理伦视频| 国产黄色免费看| 国产精品性| 亚洲国产一成久久精品国产成人综合| 国产区成人精品视频| 国产青榴视频| 色综合a怡红院怡红院首页| 精品无码人妻一区二区| 亚洲婷婷在线视频| 看你懂的巨臀中文字幕一区二区 | av无码一区二区三区在线| 久久国产亚洲偷自| 一级成人欧美一区在线观看| 亚洲天堂免费在线视频| jijzzizz老师出水喷水喷出| 国产精品无码一区二区桃花视频| 91视频免费观看网站| 亚洲国产黄色| 中文国产成人精品久久一| 色天堂无毒不卡| 伊人大杳蕉中文无码| 国产高清不卡| 一区二区在线视频免费观看|