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

基于LSTM降水量預測的咸陽“旱腰帶”地區干旱趨勢分析

2023-07-31 13:49:16馬篩艷王薇裴莉莉吳玉龍劉艷
安徽農業科學 2023年13期

馬篩艷 王薇 裴莉莉 吳玉龍 劉艷

摘要 以咸陽“旱腰帶”地區13個氣象站2010—2019年氣象監測資料為基礎,通過分析該地區降水、氣溫及干旱的特征,結合1960年以來縣氣象站資料,模擬“旱腰帶”地區近60年降水量,構建“旱腰帶”地區降水量數據集;基于長短期記憶網絡(LSTM),以歷史降水量數據為輸入對“旱腰帶”地區的降水量進行預測;根據地區干旱分級標準,由預測得到的降水量結果計算得到不同地區的干旱級別,從而實現對干旱發展趨勢的分析。結果表明,基于降水量預測數據得到的干旱分級與真實情況相比,精度達到85%以上,能夠為“旱腰帶”地區環境乃至旱區改善及生態修復提供理論基礎和決策依據。

關鍵詞 干旱發展趨勢;干旱等級;降水量預測;致旱原因;LSTM

中圖分類號 S162? 文獻標識碼 A? 文章編號 0517-6611(2023)13-0192-06

doi:10.3969/j.issn.0517-6611.2023.13.044

Analysis of Drought Trend in Xianyang “Drought Belt”Area Based on LSTM Precipitation Forecast

MA Shai-yan1,WANG Wei1,PEI Li-li2 et al

(1.Xianyang Institute of Agricultural Meteorology,Xianyang,Shaanxi 712034;2.School of Information Engineering, Changan University, Xian,Shaanxi 710064)

Abstract Based on the meteorological monitoring data of 13 meteorological stations in the “drought belt” area of Xianyang from 2010 to 2019, by analyzing the characteristics of precipitation, temperature and drought in this area,combined with data from county meteorological stations since 1960, the precipitation in the “drought belt” area in the past 60 years was simulated, and the precipitation data set in the “drought belt” area was constructed. Based on the long short-term memory (LSTM) network, the precipitation in the “drought belt” region was predicted with historical precipitation data as input. According to the regional drought classification standard, the drought level of different regions could be calculated from the predicted precipitation results, so as to realize the analysis of the drought development trend.The results showed that the drought classification based on the precipitation forecast data had an accuracy of 85% compared with the real situation, which could provide a theoretical basis and decision-making basis for environmental improvement and ecological restoration in the “drought belt” area.

Key words Drought development trend;Drought level;Precipitation forecast;Causes of drought;LSTM

基金項目 咸陽市重大科技專項“東莊水庫區域生態環境變化與種植布局調整氣象分析研究”(2020K01-33)。

作者簡介 馬篩艷(1977—),女,陜西興平人,副研級高級工程師,碩士,從事農業氣象預報技術方法研究。

收稿日期 2022-07-23

氣象干旱主要是降水異常偏少導致的水分短缺現象,而干旱災害則是干旱自然現象和人類活動共同作用的結果[1]。干旱是危害農牧業生產的第一大自然災害,旱區氣候特點及干旱發展趨勢是旱區生態修復的理論基礎和決策依據[2],對干旱和旱區氣候條件的研究歷來都是學者們關注的熱點。

張強等[3]結合干旱研究的熱點問題和發展趨勢,分析了我國干旱研究的不足,提出了我國未來干旱研究需加強典型干旱頻發區綜合性干旱科學試驗研究。李耀輝等[4]研究表明氣候變化使我國驟發性干旱顯著增加,且在未來幾十年有可能持續下去。咸陽“旱腰帶”地區是構建北山生態屏障防止北方沙塵的最后一道防線,研究咸陽“旱腰帶”地區氣候特點,充分利用氣候資源,發展適應氣候特點的林果經濟,不僅能提高農民收入,也有利于該地區生態的逐漸修復,對實施大西安綠色城鎮化戰略意義重大[5-11]。一些學者們也已經開展了相關方面的研究,許瑋等[12]分析了“旱腰帶”地區近幾年氣溫及降水的變化特征并得出氣候呈干暖化趨勢;馬延慶等[13-16]分別研究了“旱腰帶”地區葡萄、平歐大果榛子、樹上干杏、石榴等經濟林果氣候適應性及栽培技術和管理方法等,以期充分利用氣候資源,發展適宜性高的雜果,實現鄉村振興和產業脫貧。

近年來,機器學習模型在降水量預測方面應用趨向成熟,預測呈現較高的準確度。李正方等[17]利用杭州某地10年間的部分降雨數據驗證了基于粒子群優化的CART算法對降水量預測的準確性及穩定性。舒濤等[18]通過分析小波神經網絡在趨勢預測中的優越性,結合NARX動態神經網絡對拉薩市近10年7月份的降水量進行預測及預測檢驗。朱曉霞等[19]通過分析蘭州市近年每月降雨數據與滑坡的時空分布特征及兩者關系,建立有效降水量預測模型,對黃土滑坡進行預警。LSTM是 Hochreiter等[20]1997年提出的基于循環神經網絡理論的具有記憶長短期信息能力的神經網絡,主要是為了解決長序列訓練過程中的梯度消失和梯度爆炸問題。經過若干代的發展,已經在氣候環境、溫濕度、降水量預測方面得到了廣泛的應用。

目前對“旱腰帶”地區的研究基本以短期的氣象條件為主,缺乏對長期監測數據挖掘與未來干旱趨勢預測的相關研究,對“旱腰帶”各地差異及與周邊區域分析也鮮有涉及,使得“旱腰帶”地區環境改善及生態修復缺乏理論基礎和決策依據。該研究基于2010—2019年“旱腰帶”地區氣象觀測資料,結合1960年以來縣氣象站資料,通過分析“旱腰帶”地區及其所屬縣區降水量數據特點,構建“旱腰帶”地區降水量數據集,并建立LSTM模型,預測“旱腰帶”地區未來幾年的降水量,以期為政府在“旱腰帶”地區實施生態修復、鄉村振興等方面決策提供依據。

1 “旱腰帶”地區氣象要素特征分析

1.1 數據來源

該研究所用13個“旱腰帶”地區監測站點資料來源于中國氣象局cimiss數據庫,縣氣象站資料來源于地面測報文件,數據均通過mdos進行了質控。13站位于“旱腰帶”地區13個鄉鎮,其中涇陽縣4站、禮泉縣2站、乾縣4站、三原縣3站,咸陽“旱腰帶”地區地形和站點分布如圖1所示。

因自動站11月至次年3月不觀測降水量,故區域站與自動站年降水量統計時段均為4—10月(此時段為年降水量的主要分布時段,具有一定代表性,滿足分析的需要)。因“旱腰帶”地區區域站多為兩要素站,只觀測降水量和氣溫,所以該研究干旱等級采用降水距平百分率指標計算[21],計算公式如下:

Pa=p-×100%(1)

式(1)中,Pa為降水量距平百分率(%);p為某時段的降水量(mm);為常年同期氣候平均降水量(mm)。

氣溫模擬采用線性回歸,最小二乘法計算回歸系數,確定回歸方程,并對方程進行顯著性檢驗。

1.2 “旱腰帶”地區氣象要素特征分析

1.2.1 氣象要素年變化特征。

從圖2可以看出,2010—2019年咸陽“旱腰帶”地區年平均氣溫14.1 ℃,年降水量422.6 mm,極端最高氣溫42.8 ℃,極端最低氣溫-15.9 ℃,日最大降水量99.4 mm。10年間“旱腰帶”地區氣溫呈緩慢上升趨勢,降水量呈減少趨勢。

1.2.2 氣象要素月變化特征。

從圖3可以看出,近10年咸陽“旱腰帶”地區月平均氣溫曲線呈單峰分布,1月平均氣溫最低,為-0.1 ℃,7月平均氣溫最高,為27.0 ℃。降水主要集中在7—9月,其中9月降水量最多,為105.8 mm,占年降水量的25%,4月降水量最少,為33.9 mm。

1.3 “旱腰帶”地區干旱情況特征分析

為了定量描述“旱腰帶”地區干旱程度,計算了各站近10年的干旱指數,對照氣象干旱等級表(表1),得到各站干旱等級(表2)。從歷年各站點干旱情況(表2和圖4)來看,2010—2019年“旱腰帶”地區各站出現了不同干旱類型的頻次。

從時間分布(表2和圖4)上看,2010、2011、2017、2019年 “旱腰帶”各站以無旱或輕旱為主,未出現過中旱以上等級。2016年出現3站次特旱,2站次重旱,是所有年份中特旱及重旱出現最多的一年,其次是2013年,出現了2站次特旱;2010—2014年,以無旱或者輕旱為主,特旱重旱偶有發生,5年間共出現3站次特旱,1站次重旱;2015—2019年出現4站次特旱,5站次重旱,特旱、重旱區域及發生頻次明顯增加,表現出特旱與無旱或輕旱交替出現的極端變化特征,值得注意的是,2016和2018年特旱區域較大,2017和2019年均沒有出現中旱及以上等級的干旱。

從空間分布(表2)上看,三原大程10年間發生3次重旱、1次特旱,干旱最為嚴重;乾縣梁山特旱發生頻次最高,10

年出現2次;涇陽口鎮,禮泉煙霞、昭陵,乾縣陽洪、臨平10年間以無旱或輕旱為主,未出現過中旱以上等級的干旱。

綜上所述,時間分布上,與前5年(2010—2014年)相比,近5年(2015—2019年)特旱及重旱發生區域及頻次均有明顯擴大趨勢,同時表現出特旱與無旱或輕旱交替出現的極端變化特征。空間分布上,三原陵前、馬額、大程3站特旱及重旱發生頻次較高,干旱尤為嚴重;禮泉煙霞、昭陵2站未發生過中旱以上等級干旱,干旱最為緩和;涇陽及乾縣內部差異較大,部分地區干旱等級較高,部分地區則較低。

2 “旱腰帶”地區降水量數據集構建

2.1 縣站降水量與“旱腰帶”站點降水量對應分析

對2010—2019年“旱腰帶”各站與對應縣站降水量進行分析,發現各站與所屬縣站降水量變化趨勢較為一致,因此利用1960年以來縣站降水量變化趨勢來表征“旱腰帶”各地降水量變化趨勢是可行的。從1960—2019年降水量變化(圖5)可以看出,自1960年以來,涇陽、乾縣、禮泉3站降水量均呈減少趨勢,其中涇陽減少最為明顯,約為0.2 mm/10 a,其次是乾縣,禮泉無明顯減少趨勢。三原站1970年建站,1988—2005年機構改革業務調整,取消了觀測,所以此時段資料缺測,但從其他時段來看,變化趨勢與其他3站基本相同。

2.2 數據集構建

原始數據集包含區域站日降水量與4縣站月降水量,時空尺度對齊以及數據修復是構建數據集的關鍵步驟。

2.2.1 時空尺度對齊。

4縣站降水量數據包含涇陽、禮泉、乾縣、三原4縣的數據,區域站降水量數據包含安吳、云陽、口鎮、王橋、煙霞、昭陵、注泔、梁山、陽洪、臨平、陵前、馬額、大程13個站點的數據,通過地理位置信息關聯,將區域站降水量數據與4縣站降水量數據相匹配,構建新的數據集,例如,云陽監測站位于涇陽縣,因此云陽監測站2010年之前的數據可以用涇陽縣的降水量數據代替;梁山監測站位于乾縣,因此梁山監測站2010年之前的數據可以用乾縣的降水量數據代替。

4縣站降水量數據集起始時間從1955年至1959年不等,同時這些年數據較早,數據質量并不高,存在大量缺失的情況,因此選取1960年為降水量數據起始時間。區域站降水量數據集均為2010—2020年,其中2020年數據不完整,選取2010—2019年數據。將1960—2009年的4縣站降水量數據與2010—2019年的區域站降水量數據對應整合,構建時間為1960—2019年的“旱腰帶”區域監測站降水量數據集。

此外,4縣站降水量數據為每年4—10月的月降水量數據,區域站降水量數據為完整的全年日降水量數據,因此這里對區域站數據進行降采樣,即提取每年4—10月的降水量數據。最后將4縣站和區域站的降水量數據按年尺度進行求和,得到年降水量。

2.2.2 數據篩選與修復。

同時該數據集中包含缺失值,需要首先對數據進行篩選與修復,減少對預測精度的影響。對于有部分數據缺失的區域站,對整體的年降水量及干旱等級影響較小,將缺失值補充為0。對于有大量數據缺失或數據連續缺失過多的區域站,直接將其刪除。最終選取涇陽云陽、乾縣梁山、乾縣注泔、乾縣陽洪、乾縣臨平5個區域站近60年的年降水量數據構建數據集。

3 基于LSTM模型降水量預測的“旱腰帶”地區干旱趨勢分析

3.1 LSTM模型

LSTM網絡屬于遞歸神經網絡(recurrent neural network,RNN),傳統的前饋神經網絡只能從輸入節點接收數據,數據只能從輸入層到隱藏層,然后再到輸出層。前饋神經網絡和RNN之間的關鍵區別在于,它們除了在內部狀態空間上執行操作外,還能夠在輸入空間上進行操作。需要強調的是,RNN模型有梯度消失和梯度爆炸2個缺點。Hochreiter等[20]通過引入輸入和輸出門的概念,開發了LSTM架構來克服這些問題。LSTM模型的工作原理類似于RNN。然而,它們的區別在于,除了被稱為細胞[22]的內部處理單元外,循環神經元的使用門被稱為“遺忘門”“更新門”和“輸出門”。每個門都負責特定的任務,“遺忘門”忘記部分過去的信息,“輸入門”記住部分現在的信息,然后將過去的記憶與現在的記憶合并后通過“輸出門”決定最終輸出的部分。

(1)遺忘階段。

這個階段主要是對上一個節點傳進來的輸入進行選擇性忘記,計算公式如下:

ft=σ(Wf×[Ht-1,Xt]+bf)(2)

式(2)中,ft表示遺忘門輸出,Ht-1表示上一個節點的輸出,Xt表示輸入,Wf表示遺忘門權重,bf表示遺忘門偏置,σ表示 Sigmoid 層。

(2)輸入階段。

這個階段將該節點的輸入有選擇性地進行“記憶”,計算公式如下:

it=σ(Wi×[Ht-1,Xt]+bi)(3)

式(3)中,it表示輸入門輸出,Ht-1表示上一個節點的輸出,Xt表示輸入,Wi表示輸入門權重,bi表示輸入門偏置,σ表示Sigmoid層。

當前候選節點狀態計算公式如下:

t=tanh(Wc×[Ht-1,Xt]+bc)(4)

式(4)中,t表示當前候選節點狀態,Ht-1表示上一個節點的輸出,Xt表示輸入,Wc表示權重,bc表示偏置,tanh表示tanh層。

對前2個階段得到的結果進行計算,計算公式如下:

Ct=ft×Ct-1+it×t(5)

式(5)中,Ct表示傳輸給下一個節點的狀態,Ct-1表示上一個節點的狀態。

(3)輸出階段。

這個階段將決定哪些信息將會被當成當前狀態的輸出,輸出門計算邏輯:

ot=σ(Wo×[Ht-1,Xt]+bo)(6)

式(6)中,ot表示輸出門輸出,Ht-1表示上一個節點的輸出,Xt表示輸入,Wo表示輸出門權重,bo表示輸出門偏置,σ表示 Sigmoid 層。

當前節點最終輸出對上一階段得到的Ct采用tanh 函數進行了放縮,計算公式如下:

Ht=ot×tanh(Ct)(7)

式(7)中,Ht表示當前節點的最終輸出。

3.2 降水量預測結果與分析

以1960—2019年“旱腰帶”區域監測站降水量數據集為模型輸入,通過構建LSTM模型,以5年的數據為一組,對下一年數據進行預測。使用訓練好的LSTM模型預測2016—2019年的降水量,并與真實數據進行對比。以R2、均方根誤差方差(RMSE)和平均絕對誤差(MAE)分別作為精度和誤差的評價指標,部分預測結果如表3所示。

其中,乾縣梁山監測站與乾縣注泔監測站的預測數據與真實數據對比如圖6所示。從圖6可以看出,該模型對于降水量的預測與真實數據的擬合程度較高,說明模型能夠較為準確地對降水量未來發展趨勢進行預報,從而實現對干旱等級發展趨勢的分析。在對乾縣注泔區域站和乾縣梁山區域站降水量的預測中,4年降水量預測結果的平均絕對誤差(MAE)分別為28.18和53.63 mm(表3)。

3.3 基于降水量預測結果的干旱趨勢分析

根據干旱等級計算公式,選取1981—2010年的平均降水量作為同期氣候平均降水量,將2016—2019年預測降水量對應的干旱等級與真實干旱等級相比較,對比結果如表4所示。

該研究選取5個站點2016—2019年的數據用于測試模型精度,20條數據中17條數據預測正確,準確度達到85%。從區域站角度分析預測結果(表4),乾縣臨平、乾縣注泔4年的干旱等級預測準確率達到100%,涇陽云陽、乾縣梁山、乾縣陽洪4年的干旱等級預測準確率均為75%。從干旱等級角度分析預測結果,輕旱、無旱等級預測準確率達到100%;模型只是較少出現的重旱、特旱等級難以完全預測準確。

4 結論

基于2010—2019年咸陽“旱腰帶”地區及周邊區域的氣象觀測數據,分析“旱腰帶”地區及周邊區域的干旱特征,結果發現,在時間分布方面,近5年特旱及重旱發生區域及頻次均有明顯的增加趨勢,同時表現出特旱與無旱或輕旱交替出現的極端變化特征。在空間分布方面,三原縣特旱及重旱發生頻次較高,干旱尤為嚴重,禮泉干旱最為緩和,涇陽及乾縣內部差異較大,部分地區干旱等級較高,部分地區則較低。各地主要致旱原因略有差異,但降水量下降是共同特點。

以咸陽“旱腰帶”地區周邊區縣的氣象站資料模擬該地區近60年降水量,建立基于LSTM降水量預測模型,對干旱趨勢進行預測分析。雖然降水量預測值與真實值存在一定差異,但是在進行干旱等級計算后,基于降水量預測數據得到的干旱分級與真實情況相比,精度達到85%,能夠為“旱腰帶”地區環境改善及生態修復提供理論基礎和決策依據。

該研究主要分析了“旱腰帶”地區氣象干旱等級及基于LSTM降水量預測模型的干旱發展趨勢,對咸陽“旱腰帶”地區干旱發展有一定指示意義,但由于區域自動站多為兩要素站,僅觀測降水、氣溫兩項,故該研究中干旱級別采用降水距平百分率來劃分,而降水距平百分率在表征干旱時只考慮降水因素,未考慮蒸發等因素,是該研究的不足及有待改進之處。

參考文獻

[1] 張劍明,段麗潔.2013年夏季湖南省持續高溫干旱變化特征及其成因分析[J].氣象與環境學報,2018,34(4):45-51.

[2] 馬永忠,蔡福,趙先麗,等.花后持續干旱對玉米生理參數及產量的影響[J].氣象與環境學報,2019,35(2):102-106.

[3] 張強,姚玉璧,李耀輝,等.中國干旱事件成因和變化規律的研究進展與展望[J].氣象學報,2020,78(3):500-521.

[4] 李耀輝,周廣勝,袁星,等.干旱氣象科學研究:“我國北方干旱致災過程及機理”項目概述與主要進展[J].干旱氣象,2017,35(2):165-174.

[5] 蔡福,明惠青,謝艷兵,等.東北地區春玉米關鍵生育期干旱對根系生長的影響[J].氣象與環境學報,2018,34(2):75-81.

[6] 杜和平.咸陽市渭北旱腰帶水土流失危害及防治對策[J].農業科技與信息,2016(16):44-45,47.

[7] 李健偉,朱引團.渭北“旱腰帶”的呼喚:三秦環保世紀行渭北“旱腰帶”石灰石開采影響生態環境問題追蹤報道[J].法治與社會,2012(10):16-18.

[8] 朱鋆.生態文明引領下的環境脆弱地區多維度發展路徑:以咸陽市旱腰帶核心地區為例[J].規劃師,2017,33(S2):164-168.

[9] 中共咸陽市委 咸陽市人民政府關于印發《咸陽市“旱腰帶”特困片區扶貧攻堅規劃(2011年~2020年)》的通知[A].2014.

[10] 杜馳曄,王芳,孫麗楠,等.“旱腰帶”地區產業扶貧的調查與思考[R].2020.

[11] 馮雪.以五大發展理念為引領推進旱腰帶地區農業強 農村美 農民富[J].新西部,2018(32):26-27.

[12] 許瑋,趙婧,方永俠,等.咸陽渭北旱腰帶地區氣候特征分析[J].咸陽師范學院學報,2020,35(4):73-76.

[13] 馬延慶,范建璋,劉長民,等.陜西咸陽渭北“旱腰帶”優質葡萄基地生態條件分析與對策[J].中外葡萄與葡萄酒,2007,32(5):33-36.

[14] 劉建海,阮班錄,仝玉琴,等.平歐大果榛子在咸陽旱腰帶地區栽培管理技術[J].陜西農業科學,2020,66(12):95-98.

[15] 何曉光,陳明杰,孫丙寅,等.渭北“旱腰帶”樹上干杏栽培技術研究[J].陜西農業科學,2020,66(11):64-66.

[16] 吳雙,姜麗霞,李宇光,等.基于自然災害風險理論的黑龍江省玉米干旱風險評價[J].氣象與環境學報,2019,35(6):139-144.

[17] 李正方,杜景林,周蕓.基于改進CART算法的降雨量預測模型[J].現代電子技術,2020,43(2):133-137,141.

[18] 舒濤,葉唐進,李俊杰,等.降雨量及疊加預測方法研究[J].高原氣象,2021,40(1):169-177.

[19] 朱曉霞,張力,楊樹文.降雨引發的蘭州黃土滑坡時空規律分析和臨界降雨量預測[J].中國地質災害與防治學報,2019,30(4):24-31.

[20] HOCHREITER S,SCHMIDHUBER J.Long short-term memory[J].Neural computation,1997,9(8):1735-1780.

[21] 楊霏云,鄭秋紅,羅蔣梅,等.實用農業氣象指標[M].北京:氣象出版社,2015.

[22] SUTSKEVER I,VINYALS O,LE Q V.Sequence to sequence learning with neural networks[C]//The 27th international conference on neural information processing systems(NIPS14).Cambridge,MA:MIT Press,2014.

主站蜘蛛池模板: 亚洲欧美日韩中文字幕在线| 国产一区二区人大臿蕉香蕉| 国产欧美日韩视频一区二区三区| 亚洲天堂首页| 5555国产在线观看| 久久精品视频一| 国产在线视频二区| 国产精品人成在线播放| 人妻精品全国免费视频| 亚洲日韩Av中文字幕无码| 99re视频在线| 国产杨幂丝袜av在线播放| 国产真实乱人视频| 毛片网站在线播放| 香蕉综合在线视频91| 国产波多野结衣中文在线播放| 国产精品美女在线| 日本不卡在线| 欧美综合区自拍亚洲综合天堂 | 亚洲美女操| 成人免费一区二区三区| 六月婷婷激情综合| 四虎影院国产| 国产呦精品一区二区三区下载| 欧美成a人片在线观看| 婷婷五月在线| 免费网站成人亚洲| 欧美精品亚洲二区| 亚洲a免费| 人妻一本久道久久综合久久鬼色| 亚洲无码37.| 一本色道久久88综合日韩精品| 99久久无色码中文字幕| 一级毛片免费观看不卡视频| 亚洲av无码成人专区| 91久久国产热精品免费| 国产精品永久久久久| 欧美福利在线观看| 欧美日韩中文国产va另类| 国产精品久久久久鬼色| 欧美一级视频免费| a级毛片毛片免费观看久潮| 亚洲日本中文综合在线| 在线永久免费观看的毛片| 国产精品久久久精品三级| 999国产精品| 国产精品久久久精品三级| 毛片免费观看视频| 超碰色了色| 久久情精品国产品免费| 午夜a级毛片| 久久精品国产在热久久2019| 青青青草国产| 澳门av无码| 天天色综合4| 亚洲欧美日韩色图| 国内精品久久人妻无码大片高| 成AV人片一区二区三区久久| 欧美日本一区二区三区免费| 99久久精品国产综合婷婷| 91精品国产无线乱码在线| 久久99国产乱子伦精品免| 高清无码手机在线观看| 黄色一及毛片| 国产精品尤物在线| 国产成人综合亚洲网址| 玖玖免费视频在线观看| 天天综合网色中文字幕| 午夜福利网址| 国产欧美高清| 欧美成人亚洲综合精品欧美激情| 日本91视频| 666精品国产精品亚洲| 无码精品国产VA在线观看DVD| 四虎精品黑人视频| 日本免费福利视频| 91蜜芽尤物福利在线观看| 成人午夜视频在线| 天堂va亚洲va欧美va国产 | 亚洲最大福利网站| 国产91久久久久久| 国产成人91精品|