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

基于多模式集成輸出天氣變量的參考作物騰發量預報

2023-05-15 03:37:18常曉敏李攀登魏科宇左廣宇
農業工程學報 2023年5期
關鍵詞:模型

常曉敏,李攀登,魏科宇,左廣宇

基于多模式集成輸出天氣變量的參考作物騰發量預報

常曉敏1,李攀登1,魏科宇1,左廣宇2

(1. 太原理工大學水利科學與工程學院,太原 030024; 2. 太原理工大學電氣與動力工程學院,太原 030024)

參考作物騰發量(reference evapotranspiration, ET0)是農業生產中一項重要的參數,對評估未來的干旱程度和實現農業精細化管理具有重要意義。為進一步提高ET0的預報精度,該研究將多模式集成方法應用于ET0的預報,運用遺傳算法-回歸型支持向量機對歐洲中期天氣預報中心、美國國家環境預報中心、日本氣象廳和韓國氣象廳4個中心全球集合預報模式輸出的天氣變量進行多模式集成處理,基于最優的模式和方案使用Penman-Monteith公式對山西運城站未來1~7 d的ET0進行預報,并對其在站點附近農業試驗田的適用性進行驗證。結果表明,多模式集成能夠調和單一模式在氣象預報中的優劣,從而提高ET0預報的精度和長預見期下的穩定性;在ET0預報中,多模式方案的性能明顯優于原始單一模式,由最優模式和方案組成的重組方案預報性能最好,具有最小的均方根誤差、平均絕對百分比誤差,分別為0.65~0.81mm/d和19.43%~23.78%,以及最高的決定系數(0.83~0.89)。在對試驗田未來1~7 d的ET0預報中,重組方案仍表現出良好的預報性能,均方根誤差、平均絕對百分比誤差不超過0.83 mm/d和34.57%。該研究能有效提升數值天氣預報在運城站下屬鄉鎮地區的適應性,為當地農業實際生產提供準確的ET0預報信息,對于農業需水預測以及水資源優化管理具有重要意義。

騰發量;支持向量機;遺傳算法;Penman-Monteith公式;氣象參數;ET0預見期

0 引 言

作物騰發量(evapotranspiration, ET)是指在作物生長過程中土壤蒸發和作物蒸騰所消耗的水量,它對評估未來的干旱程度和實現農業精細化管理具有重要意義。雖然可通過蒸滲儀等儀器對ET進行測量,但其使用條件十分嚴格,不適用于農業實際生產。更常規的做法是通過作物系數和參考作物騰發量(reference evapotranspiration, ET0)間接估計ET。因此,ET0是否能準確計算決定了ET的估算精度。

聯合國糧農組織規定FAO-56 Penman-Monteith(P-M)模型為計算ET0的標準方法[1]。該方法以氣象資料為輸入,在不同區域和氣候條件下都有著較高的準確率,常作為校準其他模型的標準方法[2-3]。目前已有大量基于P-M模型估算ET0的研究。徐俊增等[4]基于實測ET0值,評估了11種日ET0計算方法的計算結果,指出P-M方法與實測值最為接近。袁小環等[5]分析了P-M模型計算值與實測值在不同天氣以及不同尺度下的差異性,得出P-M模型在北京地區有較好的適用性。但是,由于區域間氣象條件的差異,P-M模型在不同地區表現出不同的適用性[6]。強小嫚等[7]在評估P-M模型、FAO-17 Modified-Penman(FAO-MP)模型、Modified Penman(PBX-MP)模型以及ASCE Penman-Monteith(ASCE-PM)模型在陜西關中地區的適用性時,指出ASCE-PM模型更適用于該地區ET0的計算,而P-M模型預報精度相對較差。此外,使用P-M模型估算ET0時,需要輸入較多的氣象觀測數據,計算較為復雜,對于缺乏氣象觀測設備的地區而言,P-M模型并不適用[8]。在此基礎上,相關學者提出利用少數天氣變量預估ET0的簡化模型,如基于氣溫預報的Hargreaves-Samani模型[9],基于氣溫和輻射的Priestley-Taylor模型[10]、Hargreaves模型[11]等。但是這些模型由于數據輸入少,對ET0的評價不夠全面,使得ET0估算精度較低[12]。因此,多數學者仍使用P-M模型作為ET0的計算方式。而如何獲取準確的天氣預報數據已成為研究的關鍵所在。

近年來,隨著數值天氣預報的主要過程和現象在分辨率、參數化和物理表示方面的改進,使天氣預報的性能不斷提升[13],基于天氣預報數據的間接方法被廣泛應用于ET0日預報。數值天氣預報較高的空間分辨率提供了捕捉ET0空間變化的高分辨率能力[14],在不同地區和預見期的ET0預報中都具有較高的精度[15-17]。但是,不同的數值模式在分辨率、初始場、資料同化技術、動力框架以及物理參數化方案等方面具有明顯差異,各模式的預報能力存在較大差異[18-19]。相關的研究表明[20-21],將具有不同物理、數值和初始條件的多個單一數值天氣預報模式進行多模式集成的方法,能夠傳遞物理過程和模式的不同表達,并通過組合多個模式的結果來減少單一模式初始場的不確定性及系統偏差,通常具有比任何單一模式更高的性能。

現階段,基于多模式集成的ET0預報研究還較少:RUIZ-áLVAREZ等[21]使用Hargreaves-Samani(H-S)模型和溫度數據對ET0進行預報,未能發揮出數值天氣預報多氣象變量的優勢。MEDINA等[22]評估了個別天氣變量對ET0預報性能的影響,但未分析個別天氣變量在不同預見期下的預報性能,認為直接對原始ET0進行偏差修正在計算上更有效。然而,YANG等[23]的研究證明直接校準用原始預報輸入變量構建的原始ET0預報方法不能處理從輸入變量到ET0預報的誤差傳播,與之相比,基于修正輸入變量的ET0校準預報具有更低的偏差、更高的相關系數和更高的性能。多模式集成在ET0預報中的應用進一步提高了ET0的預報精度,但不同的多模式集成和修正方法的改進效果存在差異,特別是該方法在中國的應用研究還未見報道,其在不同地區的適用性還有待研究。

為了進一步提高ET0的預報精度,本文基于歐洲中期天氣預報中心(European Centre for Medium-Range Weather Forecasts, ECMWF)、美國國家環境預報中心(National Centers for Environmental Prediction, NCEP)、日本氣象廳(Japan Meteorological Agency, JMA)和韓國氣象廳(Korea Meteorological Administration, KMA)4種天氣預報模式,使用遺傳算法-回歸型支持向量機(genetic algorithm-support vector regression, GA-SVR)對上述模式在山西運城站輸出的天氣變量進行多模式組合,并分別評估單一模式和多模式集成方案輸出的天氣變量的預報性能,選取各天氣變量最佳的預報方案。在此基礎之上,基于各方案輸出的天氣變量,使用FAO-56 P-M模型對山西運城站未來1~7 d的ET0進行預報,并對其預報性能和年內變化趨勢進行分析。此外,使用研究站點附近農業生產地區的氣象數據計算得到當地的實測ET0,用于評價ET0最佳預報方案在運城站下屬鄉鎮的準確性和實用性,以滿足當地農業實際生產的需要,為農業水資源的優化管理提供依據。

1 數據和方法

1.1 研究區概況

本文以山西運城站(35°7′N,111°4′E)為研究對象,研究區的地理位置如圖1所示。該研究區地處黃土高原東部,華北平原西部,地勢較為復雜。從氣候上看,該地區屬于半干旱大陸性季風氣候,四季分明,降雨主要集中于夏季,年平均降水量為500~600 mm,年均氣溫為11~13 ℃,年輻射總量為5 016~5 852 MJ/(m2·a),在全國屬于中等水平。試驗田位于運城市芮城縣,地處山西省的最南端。年平均氣溫為14 ℃左右,年平均降水量和無霜期分別在500 mm和200 d以上。芮城縣因其適宜的自然和地理條件,成為了中國小麥的優質產區。2020年芮城縣小麥種植面積32 046 hm2,夏糧總產19萬t以上,平均單產5 939 kg/hm2,小麥產量和質量穩居山西省第一梯隊,多年蟬聯“全國產糧大縣”稱號。因此,提高當地參考作物騰發量的預報精度,對該地區的糧食安全生產具有重要意義。

圖1 研究區地理位置

1.2 數據來源及處理

氣象數據來源于運城站2019年3月1日—2021年3月7日的歷史實測數據和數值天氣預報數據,以及試驗田2021年3月的實測數據:1)歷史實測數據來自國家氣象科學數據中心(http://www.nmic.cn),主要為運城站歷史逐日的日最高氣溫、日最低氣溫、日照時數、平均相對濕度和10 m高處的平均風速。2)數值天氣預報數據來自TIGGE數據集(https://apps.ecmwf.int/datasets/data/tigge)4個中心全球集合預報模式(ECMWF、NCEP、JMA、KMA)的控制預報產品,預見期為1~7 d,預報數據包括2 m露點溫度、過去6 h內2 m最高氣溫、過去6 h內2 m最低氣溫、凈短波太陽輻射、10 m緯向風分量和10 m經向風分量。數值天氣預報統一選取世界時間12:00作為時間起點,對應當日北京時間20:00,與歷史實測數據的觀測時間相一致。3)試驗田的氣象環境數據來自于農田小型氣象站,各傳感器參數如表1所示。

表1 農田小型氣象站傳感器參數

注:為實際風速。

Note:is the actual wind speed.

從TIGGE數據集獲得的原始數值天氣預報為二進制格點數據(GRIB2)格式,通過wgrib2.exe解碼得到相應的數據文本。該文本包含了33°~37°N、109°~113°E范圍內4個模式的數值天氣預報數據,其分辨率為0.5°×0.5°。為獲得運城站的天氣預報數據,使用反距離權重插值法(inverse distance weight,IDW)將文本數據插值到運城站。

1.3 P-M公式

用FAO-56 P-M公式[1]計算ET0如下所示:

數值天氣預報和歷史實測數據缺少實際水汽壓、凈太陽輻射以及2 m風速等數據,需要用其他氣象變量進行代替。根據FAO-56文件推薦的計算方式,可將平均相對濕度以及2 m露點溫度轉換為實際水汽壓;通過日照時數及凈短波太陽輻射得到太陽輻射(s),進而得出n值;將10 m緯向風分量、10 m經向風分量轉換為10 m平均風速,再由式(2)計算[1]得到2 m平均風速。對數值天氣預報中每日內4個時段的6 h內2 m最高氣溫和最低氣溫進行排序,得到每日的2 m最高氣溫和最低氣溫預報值。

式中10為10 m平均風速,m/s;2為2 m平均風速,m/s。

1.4 遺傳算法-回歸型支持向量機

本文選取回歸型支持向量機(support vector regression,SVR)作為多模式集成方法。SVR作為一種具有廣泛適應性的機器學習方法,預測性能明顯優于簡單集合平均法[21],同時具有非線性映射、自組織性,泛化能力強等優勢[24],能滿足對多種氣象參數預報精度的要求。SVR是將支持向量機用于回歸分析,即采用非線性映射()將低維樣本集{(x,y|=1,?,)}(其中,x為輸入項,y為輸出項)中的特征向量映射到更高維的空間中,利用不敏感損失函數,在這個空間中尋找一個最優的回歸平面(),使所有特征向量到該平面的距離最短。當采用線性回歸時,回歸函數()可表示為

式中αα為拉格朗日乘數;為偏置向量;為所支持的向量上限;為輸入項;x為第個輸入樣本的列向量。

在非線性回歸中,需在式(3)中引入一個核函數(,i)。核函數的選取是SVR建模的關鍵,而高斯徑向基核函數(radial basis function,RBF)在SVR中應用最為廣泛,同時在大樣本和小樣本的應用中都有著較好的性能。因此本文選取該核函數作為SVR模型的核函數。RBF函數計算式如下:

式中為核參數,>0。

將式(4)代入式(3),可得到SVR非線性回歸函數:

式中(,x)為引入的核函數。

在SVR模型中,懲罰因子對SVR模型的泛化能力影響較大,核參數對樣本在高維特征空間映射的復雜程度影響較大,這2個參數的取值共同決定著SVR模型性能的優劣。但是在實際應用中,SVR模型尋找懲罰因子和核參數時會出現“過學習”的現象,會影響模型的預報能力。為解決這一問題,本文采用遺傳算法(genetic algorithm, GA)合理選擇和,以提高SVR模型的性能。遺傳算法是一種基于達爾文進化論的智能優化算法,具有強大的全局并行尋優能力。

采用遺傳算法與回歸型支持向量機相結合的方法(GA-SVR)對FAO-56 P-M公式的輸入變量進行多模式集成。將多模式方案中各單一模式在1~7 d預見期內的各氣象變量(2 m最高氣溫、2 m最低氣溫、實際水汽壓、凈短波太陽輻射、10 m風速)作為GA-SVR的輸入項,將對應的歷史實測氣象變量作為輸出項。選取2019年3月1日-2020年2月29日共51 240組氣象數據構建模型的訓練集,用于模型在不同預見期及不同天氣變量中的參數調優工作;選取2020年3月1日-2021年2月28日共51 100組氣象數據作為測試集,用于評估各模式方案對天氣變量在1~7 d預見期下的預報性能。通過對輸入輸出數據進行訓練,構建結構風險最小的GA-SVR模型,然后給定單一模式下的氣象變量即可得到相應的多模式集成輸出值,如圖2所示。

注:C為懲罰因子,g為核參數。

1.5 多模式集成方案

將4個單一模式組成不同的多模式集成方案作為GA-SVR模型輸入項預報數據的來源,根據測試集結果對多模式集成方案的預報性能進行評價,以獲得各氣象參數最優的多模式集成方案。單一模式和多模式集成方案的劃分如表2所示。由于KMA模式下凈短波太陽輻射數據存在異常,故不對KMA模式下的太陽輻射數據進行多模式集成。

表2 模式集成方案

注(Note): ECMWF, European Centre for Medium-Range Weather Forecasts; NCEP, National Centers for Environmental Prediction; JMA, Japan Meteorological Agency; KMA, Korea Meteorological Administration.

1.6 統計分析指標

為驗證每個多模式集成方案是否能提升各天氣變量和ET0的預報性能,需要對各多模式集成方案輸出的天氣變量和通過這些變量所得出的ET0的預報精度進行評價。在對氣溫、實際水汽壓、風速、太陽輻射和ET0的預報中,重點關注預報的準確性與可靠性。因此,本文選取平均絕對百分比誤差(APE)、均方根誤差(MSE)、決定系數(2)3個統計指標對天氣變量和ET0的預報性能進行分析和評估。準確性主要由APE和MSE來衡量。APE是相對誤差度量值,它使用絕對值來避免正負誤差相互抵消,可以衡量預測誤差的平均大小。MSE能更好地反映樣本的離散程度,對異常值更加敏感。上述2個指標越小,預報的準確性越高。可靠性主要由2來衡量。2能夠反映出預報值與實測值之間的擬合程度。2越接近于1,預報值與實測值的擬合優度越好。各統計指標計算式見文獻[25-26]。

2 結果與分析

2.1 氣象變量預報性能分析

2.1.1 氣溫預報性能分析

日最高氣溫、日最低氣溫、實際水汽壓、太陽輻射和平均風速的預報性能統計指標如圖3所示。由圖3a~圖3f可知,日最高氣溫(max)和日最低氣溫(min)的2都隨預見期的增長呈下降趨勢,但在1~7 d的預見期內2都不小于0.90,表明氣溫預報值和觀測值之間存在較強的線性關系。在對max的預報中,單一模式的APE和MSE分別在21.17%~30.60%和3.31~4.97 ℃之間,多模式集成方案的APE和MSE分別在9.61%~20.39%和1.55~2.87 ℃之間。各單一模式的APE和MSE差別較大,且隨預見期的增長變化幅度較小。各多模式集成方案的APE和MSE差別較小,且隨預見期的增長都呈穩定增長趨勢。多模式集成方案的APE和MSE在對應預見期下始終小于4個單一模式的APE和MSE,而ENJK方案的APE和MSE均值最小,為max預報性能最優的方案。

在對min的預報中,單一模式的APE和MSE分別在44.37%~72.23%和2.42~2.95 ℃之間,多模式集成方案的APE和MSE分別在44.37%~72.95%和1.84~2.52 ℃之間。隨著預見期的增長,單一模式的MSE呈現先減小后增大的趨勢,而多模式集成方案的MSE呈現穩定增長的趨勢。不同預見期下,絕大多數多模式集成方案的APE和MSE均小于各單一模式下的APE和MSE,其中ENJK方案的APE和MSE均值最小,為min預報性能最優的方案。

從min和max預報性能的分析可以看出,單一模式下min的預報性能要優于max,這與其他研究者得出的結論相一致[27-28][63]。但通過多模式集成后,max預報的APE和MSE明顯降低。max預報性能提升較大的原因是max主要與每日云量和地表輻射強度有關[29],數值天氣預報邊界層不確定的物理參數使得max誤差較大,而多模式集成可結合各模式的初始場條件從而減小模式的系統偏差[30]。采用GA-SVR模型作為多模式集成方法能較好地提升氣溫的預報性能。

2.1.2 實際水汽壓預報性能分析

實際水汽壓(a)預報性能統計指標如圖3g~圖3i所示。單一模式的2在0.85~0.96之間,多模式集成方案的2在0.91~0.96之間,表明a預報值和觀測值之間存在較強的線性關系。在預報的準確性上,單一模式的APE和MSE分別在14.33%~31.42%和0.25~0.52 kPa之間,多模式集成方案的APE和MSE分別在11.89%~29.39%和0.17~0.27 kPa之間。單一模式中,各模式的APE和MSE差別較大,且隨預見期的增長變化幅度較小。在1~7 d預見期下,K模式的預報性能最好,而N模式預報性能最差。受此影響,多模式集成方案中未包含N的EJ方案的預報性能優于其余包含N的多模式方案。由于多模式集成能較好地降低系統性誤差,各多模式方案的預報性能仍均優于各單一模式。從總體上看,1~7 d預見期下,單一模式的APE和MSE均值范圍分別為15.72%~30.86%和0.27~0.50 kPa,多模式集成方案的APE和MSE均值范圍分別為15.38%~17.76%和0.21~0.25 kPa。由此可知,多模式集成能有效提升a的預報性能。

2.1.3 太陽輻射預報性能分析

太陽輻射(s)預報性能統計指標如圖3j~圖3l所示。單一模式的2在0.64~0.78之間,多模式集成方案的2在0.70~0.79之間。隨預見期的增長,各模式和方案的2基本呈下降趨勢。單一模式中,N的2略高于E和J,且在1~4 d內尤為突出。受此影響,在1~4 d內,NJ的2最高;超過4 d后,ENJ的2最高;而EJ的2在絕大多數預見期下均為多模式集成方案的最小,但相較原始單一模式,仍有小幅提升。以上現象表明多模式集成能較好地改進s預報精度。在預報的準確性方面,單一模式的APE和MSE分別在24.26%~32.73%和3.94~4.82 MJ/m2之間,多模式集成方案的APE和MSE分別在22.00%~25.59%和4.09~4.62 MJ/m2之間。隨著預見期的增長,N和J的APE和MSE均呈上升趨勢;而E和各多模式集成方案的APE和MSE均處于不斷波動之中,無明顯的變化趨勢。在1~4 d內,N模式的APE和MSE均處于相對較低的水平;但4 d后,N模式的APE和MSE增長迅速,此時E模式的APE和MSE為單一模式中最低。多模式集成方案中,EN方案能較好地改善3 d以后s的預報性能,但改善效果依然有限。而其余多模式集成方案的APE和MSE始終與單一模式差別較小。這說明多模式集成改善s預報性能的能力有限,且受原始單一模式預報性能的影響較大。事實上,估算太陽輻射所使用的Angstr?m-Prescott(AP)模型的系數取決于研究地點的氣候和地理特征[31]。在無實測數據來校正AP系數的情況下,使用FAO-56給出的推薦系數會在一定程度上降低太陽輻射的估算精度。而估算過程中引入的較多誤差,使得s具有更大的可變性,從而超出了多模式集成的糾偏能力。為此,本研究選取單一模式中預報性能最好的N方案來預報太陽輻射值。

注:RMSE為均方根誤差,MAPE為平均絕對百分比誤差,R2為決定系數。

2.1.4 風速預報性能分析

平均風速(10)預報性能統計指標如圖3m~圖3o所示。10預報的相關性為所有考慮的氣象變量中最低,其中單一模式的2在0.08~0.32之間,多模式集成方案的2在0.11~0.32之間。隨預見期的增長,各模式和方案的2都呈迅速下降趨勢,表明風速預報值與實測值的相關性較弱,同時多模式集成對其的改善效果也十分有限,符合風速非線性、隨機性的變化規律[32]。在準確性方面,單一模式的APE和MSE分別在34.72%~47.60%和0.89~1.17 m/s之間,多模式集成方案的APE和MSE分別在29.94%~33.63%和0.71~0.81 m/s之間。隨預見期的增長,各模式和方案的APE和MSE均呈上升趨勢。在1~7 d內,各多模式集成方案的APE和MSE差別較小,且都顯著低于各單一模式。單一模式對10在1~7 d預見期下的APE和MSE均值范圍分別為37.23%~43.24%和0.96~1.10 m/s,而多模式集成方案的APE和MSE均值范圍分別為31.08%~32.42%和0.75~0.77 m/s。由此可知,多模式集成能有效改善10的預報精度。其中,ENJK方案在不同預見期下的APE和MSE均值最小,為預報性能最優的方案。

2.2 ET0預報性能分析

2.2.1 各方案預報性能比較

為進一步確定多模式集成修正的氣象變量對ET0預報性能的改善效果,本文對測試集中單一模式和多模式集成方案輸出的氣象變量進行ET0預報。除此之外,根據各氣象變量預報性能分析結果,選取各氣象變量預報性能最優的模式或多模式集成方案組成新的重組方案(restructuring,RC)進行ET0預報。RC方案的氣象變量組成如表3所示。ET0預報性能統計指標如圖4所示。

表3 重組方案

圖4 參考作物騰發量(ET0)預報性能統計指標

由圖4可知,單一模式和多模式集成方案的2分別在0.76~0.88和0.81~0.89之間,且都隨預見期的增長呈下降趨勢,這與各氣象變量相關性的變化趨勢相一致。在單一模式中,N模式的相關性最好,而J模式最差,這主要與J模式max和s等氣象變量較低的相關性有關。受此影響,EJ方案為多模式集成方案中相關性最差的方案,而其余多模式集成方案的相關性差別較小。在所有方案中,RC方案為相關性最好的方案。

在預報精度上,單一模式的APE和MSE分別在在20.13%~26.50%和0.70~1.07 mm/d之間,多模式集成方案的APE和MSE分別在19.08%~23.80%和0.65~0.84 mm/d之間。在單一模式中,N和J的預報精度分別為最高和最低,但其隨預見期的增長都下降較快;E的預報精度始終處于中等水平,隨預見期的增長變化幅度較小。而在多模式集成方案中,各方案的預報精度明顯優于各單一模式,且隨預見期的增長變化幅度較小。其中,EJ方案的預報精度最低,這與E和J預報精度的變化特點相類似。而其余包含N的多模式集成方案在1~4 d都保持著較高的預報精度,這可能與N方案在s預報中的表現有關;同樣的,因EN方案能較好地改善中長預見期下s的預報性能,所以EN方案在5~7 d的預見期下依然表現優異。以上現象也印證了太陽輻射對ET0預報性能影響最大的觀點[22]。RC方案在1~7 d預報的APE和MSE分別為19.43%~23.78%和0.65~0.81 mm/d,2為0.83~0.89。與其他方案相比,該方案在大多數預見期下的APE和MSE都最小,預報精度最高,對指導農業生產的意義重大。

綜上所述,單一模式在ET0預報中存在著預報精度較低且隨預見期的增長預報精度下降較快的問題。盡管多模式集成方案在一定程度上受到單一模式的影響,但多模式集成能夠調和單一模式在氣象預報中的優劣,從而提高ET0預報的精度和長預見期下的穩定性,實現預報精度與預報時長的統一。

2.2.2 ET0預報季節性分析

為確定不同模型和方案在ET0預報中季節性變化的特點,本文基于測試集的評估結果,對重組方案(RC)、最優單一模式(N)和最優多模式集成方案(EN)在不同預見期(短、中、長即1、4、7 d)下ET0預報的年內變化趨勢進行分析,并按季節對預報精度進行評價。ET0年內變化趨勢和預報精度評價結果分別如圖5和表4所示。

圖5 不同預見期下ET0實測值和預報值在年內的變化趨勢

由圖5可知,不同預見期下ET0預報的年內變化趨勢基本相同,ET0隨季節變化的趨勢很明顯。從春季的3月開始,隨著氣溫的升高和日照時長的增加,ET0處于不斷升高的趨勢。但由于春季升溫過程中氣候變化復雜,大范圍寒潮的發生使得ET0不斷波動[33],最終在5月末到達最高值,該值也為全年的最高值。到6月進入夏季后,較春末的ET0略有下降。同時降水的頻率和時長增多,特別是極端天氣的出現,使得太陽輻射值較低,使得ET0大幅波動。在9月時,ET0開始出現下降的趨勢;在之后的9月末到11月初,一直在2 mm/d附近波動。直到11月中旬,ET0出現了下降,并最終維持在低于1 mm/d的水平。進入冬季后,由于氣溫低,同時天氣變化情況變化較小,使得ET0仍舊低于1 mm/d。在1月中旬,隨著氣溫的回升,ET0也出現了小幅的增長,在2月中下旬達到最高值4 mm/d。之后降雨和寒潮天氣的出現使得ET0再次降低。

從預報精度上來看(表4),在春季中,RC方案在不同預見期下各預報性能統計指標均為最優。在4、7 d的預見期下,EN方案與RC方案在APE和MSE上差別很小。而N模式在所有預見期下的預報性能都大幅劣于其他2個方案,特別是在4、7 d預見期下的2相對較低。在夏季中,ET0對太陽輻射的敏感系數達到頂峰[34],而夏季大幅度的天氣變化使得更難以對ET0進行準確預報。從總體上來看,夏季的APE和MSE均大幅高于春季,2也大幅降低。隨著預見期的增長,3種方案的預報性能均呈現降低的趨勢。進入秋季后,3種方案的MSE大幅下降,2也有所提升,而APE與夏季相比變幅不大。各方案之間的預報性能差別不大。隨著預見期的增長,各模式的預報性能均有小幅下降,EN方案的預報性能要優于N和RC方案。冬季中,不同預見期下,N與EN方案的預報性能差別不大,而RC方案的APE和MSE均為最小,2也最高,為冬季3個方案中最優的方案。

整體上,各方案對ET0預報性能最好的季節是秋季,其次是冬季、春季和夏季。RC方案除夏季的長預見期外,全年的預報性能都十分優異,在春季和冬季的預報性能為所有方案中最優,秋季的預報精度也都保持在較高的水平。EN方案在秋季中長預見期下有著很好的預報性能。N方案在短預見期下的預報性能與其他方案差別較小,但在中長預見期的預報性能下降較快,為預報性能最差的方案。

表4 不同季節各模型不同預見期下預報精度評價指標

2.2.3 ET0最優預報方案在試驗田的應用

采用多模式集成方法進行ET0預報,需要大量的歷史實測氣象數據進行訓練。但通常只能獲取國家級地面觀測站的氣象數據,而缺乏實際農業生產地區的實測氣象數據,所以無法直接對農業生產地區的ET0進行預報。但ET0的變化主要取決于氣候和地理因素[35]。為此,本文選取位于運城市芮城縣的冬小麥地作為試驗田,通過農田小型氣象站獲取試驗田氣象數據,利用風速、空氣溫濕度和太陽總輻射等氣象環境參數對該地的ET0進行估算,并將該值作為試驗田的ET0實際值對ET0預報的精度進行評估,以驗證運城站的ET0預報結果是否適合實際農業生產地區。ET0預報所需的氣象參數來自于預報性能最優的多模式集成方案——重組方案。選取試驗田2021年3月的氣象數據對ET0的預報性能進行分析。

ET0預報性能統計分析結果如圖6所示。試驗田在2021年3月1~7 d預見期下ET0預報的APE和MSE分別在21.68%~34.57%和0.56~0.83 mm/d之間,2則在0.35~0.62之間。隨預見期的增長,APE和MSE呈上升趨勢,2呈下降趨勢。與重組方案在測試集中同一月份的表現相比,試驗田的2明顯下降,相關程度降低,這主要是試驗田與運城站地理位置的差異造成的。但APE和MSE未見明顯增長。已有的研究表明[36],利用氣象參數預報ET0的MSE在1.0 mm/d左右,說明重組方案在試驗田的ET0預報中仍具有較高的精度,依然能夠滿足當地農業實際生產的需要。為減小地理位置差異對ET0預測精度的影響,在今后的研究中,可通過農田小型氣象站獲取試驗田足夠的氣象數據來對模型進行訓練,以此提高當地ET0的預測精度。

圖6 2021年3月ET0預報性能統計指標

3 結 論

本文使用遺傳算法-支持向量機(genetic algorithm- support vector regression, GA-SVR)對歐洲中期天氣預報中心、美國國家環境預報中心、日本氣象廳和韓國氣象廳4種天氣預報模式1~7 d預見期下的天氣變量進行多模式集成,基于優化后的天氣變量對運城站的參考作物騰發量(reference evapotranspiration, ET0)進行預報,分析其預報性能,并驗證了其在實際農業生產地區的適用性。主要結論如下:

1)多模式集成能夠較好地改進部分氣象變量的預報性能,對日最高氣溫的改進幅度最大,其次為實際水汽壓、10 m風速、日最低氣溫,對入射太陽短波輻射預報性能的改進不明顯。

2)多模式集成能夠調和單一模式在氣象預報中的優劣,從而提高ET0預報的精度和長預見期下的穩定性,實現預報精度與預報時長的統一。

3)根據各氣象變量優化結果組成的重組方案在ET0預報中表現最優,該方案在1~7 d預報的均方根誤差為0.65~0.81 mm/d,平均絕對百分比誤差為19.43%~23.78%,決定系數為0.83~0.89。

4)重組方案在對試驗田的ET0預報中2下降幅度較大,但均方根誤差、平均絕對百分比誤差未有明顯增長,其值不超過0.83 mm/d和34.57%,能夠滿足當地農業實際生產的需要。

如今,數值天氣預報已成為全球許多氣象中心的主要預測工具。本文提出的多模式集成方案一方面能顯著降低氣象參數和ET0預測的誤差,彌補單一模式集成造成的過度自信,保持預測值和實測值之間的相關性,另一方面也保證了數值天氣預報在一定范圍內的適應性(在預報的空間尺度上,數值天氣預報的最小預報單位為縣級,這使得區域天氣預報對不同的地形條件適應性差)。然而多模式集成方案對于ET0的預測尚有不完善之處。使用的全球集合預報模式數據和國家氣象科學數據中心的實測數據在部分氣象變量上存在差異,通過相關模型和公式進行推算容易引入較大的誤差,可能會超出多模式集成的糾偏能力,從而影響ET0的預報精度。在今后的研究中應通過建立的農田小型氣象站積累足夠的氣象環境數據來對模型進行訓練,從而提高ET0的預測精度。此外,由于試驗站點較少,本文提出的多模式集成方案目前僅對運城站周圍縣區、鄉鎮的ET0預測具有一定的適應性,對于其他農業生產地區和季節的ET0預報性能和適用性還有待驗證。

[1] ALLEN R G, PEREIRA L S, RAES D, et al. Crop evapotranspiration-Guidelines for computing crop water requirements[R]//FAO Irrigation and Drainage. Paper 56. Roma: FAO, 1998.

[2] 劉夢,羅玉峰,汪文超,等. 基于天氣預報的漳河灌區參考作物騰發量預報方法比較[J]. 農業工程學報,2017,33(19):156-162.

LIU Meng, LUO Yufeng, WANG Wenchao, et al. Comparison of three reference crop evapotranspiration forecasting methods based on short-term weather forecast in Zhanghe irrigation district[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(19): 156-162. (in Chinese with English abstract)

[3] 陳紹民,李曉麗,楊啟良,等. 基于機器學習的遮蔭設施內參考作物蒸散量估算[J]. 農業工程學報,2022,38(11):108-116.

CHEN Shaomin, LI Xiaoli, YANG Qiliang, et al. Estimation of reference evapotranspiration in shading facility using machine learning[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2022, 38(11): 108-116. (in Chinese with English abstract)

[4] 徐俊增,彭世彰,丁加麗,等. 基于蒸滲儀實測數據的日參考作物蒸發騰發量計算方法評價[J]. 水利學報,2010,41(12):1497-1505.

XU Junzeng, PENG Shizhang, DING Jiali, et al. Evaluation of methods for estimating daily reference crop evapotranspiration based on lysimeter grass experiments[J]. Journal of Hydraulic Engineering, 2010, 41(12): 1497-1505. (in Chinese with English abstract)

[5] 袁小環,滕文軍,張輝,等. 實測草坪蒸散量評價P-M模型在北京地區適用性[J]. 農業工程學報,2018,34(7):147-154.

YUAN Xiaohuan, TENG Wenjun, ZHANG Hui, et al. Suitability assessment of P-M model by measuring ET0of turfs in Beijing, China[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2018, 34(7): 147-154. (in Chinese with English abstract)

[6] 袁小環,楊學軍,陳超,等. 基于蒸滲儀實測的參考作物蒸散發模型北京地區適用性評價[J]. 農業工程學報,2014,30(13):104-110.

YUAN Xiaohuan, YANG Xuejun, CHEN Chao, et al. Applicability assessment of reference evapotranspiration models in Beijing based on lysimeter measurement[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2014, 30(13): 104-110. (in Chinese with English abstract)

[7] 強小嫚,蔡煥杰,孫景生,等. 陜西關中地區ET0計算公式的適用性評價[J]. 農業工程學報,2012,28(20):121-127.

QIANG Xiaoman, CAI Huanjie, SUN Jingsheng, et al. Adaptability evaluation for reference evapotranspiration (ET0) formulas in Guanzhong Region of Shaanxi[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2012, 28(20): 121-127. (in Chinese with English abstract)

[8] YASSIN M A, ALAZBA A A, MATTAR M A. Modelling daily evapotranspiration using artificial neural networks under hyper arid conditions[J]. Pakistan Journal of Agricultural Sciences, 2016, 53(3): 695-712.

[9] LUO Y, CHANG X, PENG S, et al. Short-term forecasting of daily reference evapotranspiration using the Hargreaves–Samani model and temperature forecasts[J]. Agricultural Water Management, 2014, 136: 42-51.

[10] 張寶忠,許迪,劉鈺,等. 多尺度蒸散發估測與時空尺度拓展方法研究進展[J]. 農業工程學報,2015,31(6):8-16.

ZHANG Baozhong, XU Di, LIU Yu, et al. Review of multi-scale evapotranspiration estimation and spatio-temporal scale expansion[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(6): 8-16. (in Chinese with English abstract)

[11] 賈悅,崔寧博,魏新平,等. 考慮輻射改進Hargreaves模型計算川中丘陵區參考作物蒸散量[J]. 農業工程學報,2016,32(21):152-160.

JIA Yue, CUI Ningbo, WEI Xinping, et al. Modifying Hargreaves model considering radiation to calculate reference crop evapotranspiration in hilly area of central Sichuan Basin[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2016, 32(21): 152-160. (in Chinese with English abstract)

[12] 馮克鵬,田軍倉,洪陽. 自尋優最近鄰算法估算有限氣象數據區潛在蒸散量[J]. 農業工程學報,2019,35(20):76-83.

FENG Kepeng, TIAN Juncang, HONG Yang. Method for estimating potential evapotranspiration by self-optimizing nearest neighbor algorithm[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2019, 35(20): 76-83. (in Chinese with English abstract)

[13] HAMILL T M, BATES G T, WHITAKER J S, et al. NOAA's second-generation global medium-range ensemble reforecast dataset[J]. Bulletin of the American Meteorological Society, 2013, 94(10): 1553-1565.

[14] LIU B, LIU M, CUI Y, et al. Assessing forecasting performance of daily reference evapotranspiration using public weather forecast and numerical weather prediction[J]. Journal of Hydrology, 2020, 590: 125547.

[15] ISHAK A M, BRAY M, REMESAN R, et al. Estimating reference evapotranspiration using numerical weather modelling[J]. Hydrological Processes, 2010, 24: 3490-3509.

[16] SILVA D, MEZA F J, VARAS E. Estimating reference evapotranspiration (ET0) using numerical weather forecast data in central Chile[J]. Journal of Hydrology, 2010, 382(1/2/3/4): 64-71.

[17] PERERA K C, WESTERN A W, NAWARATHNA B, et al. Forecasting daily reference evapotranspiration for Australia using numerical weather prediction outputs[J]. Agricultural and Forest Meteorology, 2014, 194: 50-63.

[18] 楊學勝. 業務集合預報系統的現狀及展望[J]. 氣象,2001,27(6):3-9.

YANG Xuesheng. The new development and the outlook of the operational ensemble prediction system[J].Meteorological Monthly, 2001, 27(6): 3-9. (in Chinese with English abstract)

[19] 智協飛,王田,季焱. 基于深度學習的中國地面氣溫的多模式集成預報研究[J]. 大氣科學學報,2020,43(3):435-446.

ZHI Xiefei, WANG Tian, JI Yan. Multimodel ensemble forecasts of surface air temperature over China based on deep learning approach[J]. Transactions of Atmospheric Sciences, 2020, 43(3): 435-446. (in Chinese with English abstract)

[20] HAGEDORN R, DOBLAS-REYES F J, PALMER T N. The rationale behind the success of multi-model ensembles in seasonal forecasting-I. Basic concept[J]. Tellus Series a-Dynamic Meteorology and Oceanography, 2005, 57(3): 219-233.

[21] RUIZ-ALVAREZ M, GOMARIZ-CASTILLO F, ALONSO- SARRIA F. Evapotranspiration response to climate change in semi-arid areas: using random forest as multi-model ensemble method[J]. Water, 2021, 13(2): 222.

[22] MEDINA H, TIAN D, SRIVASTAVA P, et al. Medium-range reference evapotranspiration forecasts for the contiguous United States based on multi-model numerical weather predictions[J]. Journal of Hydrology, 2018, 562: 502-517.

[23] YANG Q, WANG Q J, HAKALA K, et al. Bias-correcting input variables enhances forecasting of reference crop evapotranspiration[J]. Hydrology and Earth System Sciences, 2021, 25(9): 4773-4788.

[24] CHO D, YOO C, IM J, et al. Comparative assessment of various machine learning-based bias correction methods for numerical weather prediction model forecasts of extreme air temperatures in urban areas[J]. Earth and Space Science, 2020, 7(4):e2019EA000740.

[25] 曹雯,段春鋒,徐祥,等. 基于BCCP_CPSv2模式的淮河流域月參考作物蒸散概率訂正預報[J]. 農業工程學報,2022,38(23):61-69.

CAO Wen, DUAN Chunfeng, XU Xiang, et al. Probability correction in monthly reference crop evapotranspiration prediction in Huaihe River Basin using BCC_CPSv2 model[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2022, 38(23): 61-69. (in Chinese with English abstract)

[26] 江顯群,陳武奮,邵金龍. 基于公共天氣預報的參考作物騰發量預報[J]. 排灌機械工程學報,2019,37(12):1077-1081.

JIANG Xianqun, CHEN Wufen, SHAO Jinlong. Forecast of reference crop evapotranspiration based on public weather forecast[J]. Journal of Drainage and Irrigation Machinery Engineering (JDIME), 2019, 37(12): 1077-1081. (in Chinese with English abstract)

[27] LUO Y, TRAORE S, LYU X, et al. Medium range daily reference evapotranspiration forecasting by using ANN and public weather forecasts[J]. Water Resources Management, 2015, 29(10): 3863-3876.

[28] XIONG Y, LUO Y, WANG Y, et al. Forecasting daily reference evapotranspiration using the Blaney-Criddle model and temperature forecasts[J]. Archives of Agronomy and Soil Science, 2016, 62(6): 790-805.

[29] LOPEZ-FRANCA N, ZANINELLI P G, CARRIL A F, et al. Changes in temperature extremes for 21st century scenarios over South America derived from a multi-model ensemble of regional climate models[J]. Climate Research, 2016, 68(2/3): 151-167.

[30] 智協飛,季曉東,張璟,等. 基于TIGGE資料的地面氣溫和降水的多模式集成預報[J]. 大氣科學學報,2013(3):257-266.

ZHI Xiefei, JI Xiaodong, ZHANG Jing, et al. Multimodel ensemble forecasts of surface air temperature and precipitation using TIGGE datasets[J]. Transactions of Atmospheric Sciences, 2013(3): 257-266. (in Chinese with English abstract)

[31] MOUSAVI R, SABZIPARVAR A A, MAROFI S, et al. Calibration of the Angstrom-Prescott solar radiation model for accurate estimation of reference evapotranspiration in the absence of observed solar radiation[J]. Theoretical and Applied Climatology, 2015, 119(1/2): 43-54.

[32] YONG B, QIAO F, WANG C, et al. Ensemble Neural Network method for wind speed forecasting. Proceedings of the 33rd IEEE International Workshop on Signal Processing Systems. October 20-23, 2019[C]. Nanjing: IEEE, 2019: 31-36.

[33] 張祎瑋,李芬. 山西2020年春季氣候特征及其對農業生產的影響[J]. 農業災害研究,2020,10(4):55-57.

ZHANG Yiwei, LI Fen. The Spring Climate feature and its effect on Shanxi's agriculture in 2020[J]. Journal of Agricultural Catastrophology, 2020, 10(4): 55-57. (in Chinese with English abstract)

[34] 廖顯琴. 陜西省參考作物騰發量的時空變化及其預測[D]. 楊凌:西北農林科技大學,2010.

LIAO Xianqin. Temporal-Spatial Variations and Predictions of Reference Crop Evapotranspiration in Shaanxi Province[D]. Yangling: Northwest A&F University, 2010. (in Chinese with English abstract)

[35] SRIVASTAVA P K, HAN D, YADUVANSHI A, et al. Reference evapotranspiration retrievals from a mesoscale model based weather variables for soil moisture deficit estimation[J]. Sustainability, 2017, 9(11): 1971.

[36] 羅玉峰,李思,彭世彰,等. 基于氣溫預報和HS公式的參考作物騰發量預報[J]. 排灌機械工程學報,2013,31(11):987-992.

LUO Yufeng, LI Si, PENG Shizhang, et al. Forecasting reference crop evapotranspiration based on temperature forecast and Hargreaves-samani equation[J]. Journal of Drainage and Irrigation Machinery Engineering (JDIME), 2013, 31(11): 987-992. (in Chinese with English abstract)

Reference crop evapotranspiration forecast using multi-model integrated output weather variables

CHANG Xiaomin1, LI Pandeng1, WEI Keyu1, ZUO Guangyu2

(1.,,030024,;2.,,030024,)

Reference evapotranspiration (ET0) is one of the most important parameters to predict the drought degree in agricultural production and precision management. In this research, the multi-model ensemble method was applied to improve the forecasting precision for ET0. Genetic algorithm-support vector regression was utilized to forecast the weather variables output by European Centre for Medium-Range Weather Forecasts, National Centers for Environmental Prediction, Japan Meteorological Agency, and Korea Meteorological Administration models. The prediction accuracy and reliability of the model were evaluated by three indicators: root mean square error (MSE), mean absolute percentage error (APE), and determination coefficient (2). Penman-Monteith equation was used to forecast the ET0from 1 to 7 days in the future, according to the optimal models and schemes. Finally, the applicability of the improved model was verified in agricultural test field, including the Yuncheng Station in Shanxi Province of China. The results show that the multi-mode integration improved the prediction performance of air temperature, actual water vapor pressure, and wind speed under a single mode, with the largest improvement in the daily maximum temperature, followed by the actual water vapor pressure, wind speed, and daily minimum temperature. The prediction performance showed a downward trend with the increase in the prediction period. There was no trend of multi-mode integration on the prediction performance of solar radiation. The forecast accuracy of a single model in the ET0forecast decreased rapidly with the growth of the forecast period. By contrast, the multi-mode integration scheme improved the accuracy of ET0forecasting and the stability in a long forecast period, and then balanced the tradeoff of weather forecasting accuracy and duration. In ET0forecasts, the performance of the multi-model schemes was significantly better than that of the original single models. The best prediction performance was achieved in the forecasting precision of the recombination scheme with the optimal models. The smallestMSEandAPEwere 0.65-0.81 mm/d and 19.43%-23.78%, respectively, and the highest2was 0.83-0.89. There was excellent forecast performance of the reorganization scheme in the seasonal ET0forecast throughout the year, except for the long forecast period in summer. The forecast performance in spring and winter was the best of all the schemes, whereas, the forecast accuracy in autumn was also maintained at a high level. The reorganization scheme still showed better prediction performance in the ET0prediction of the test field in the next 1-7 days, with theMSEand the average absolute percentage error not exceeding 0.83 mm/d and 34.57%, respectively. The adaptability of numerical weather forecasts can be verified in the township areas under Yuncheng Station, providing accurate ET0forecast information for local agricultural production. It was of great significance for agricultural water demand prediction and optimal management of water resources. However, the ET0prediction performance and applicability of the multi-mode integration scheme need to be verified for the other agricultural production areas and seasons.

evapotranspiration; support vector machine; genetic algorithm; Penman-Monteith equation; meteorological parameters; ET0forecast period

10.11975/j.issn.1002-6819.202209108

S161

A

1002-6819(2023)-05-0079-11

常曉敏,李攀登,魏科宇,等. 基于多模式集成輸出天氣變量的參考作物騰發量預報[J]. 農業工程學報,2023,39(5):79-89.doi:10.11975/j.issn.1002-6819.202209108 http://www.tcsae.org

CHANG Xiaomin, LI Pandeng, WEI Keyu, et al. Reference crop evapotranspiration forecast using multi-model integrated output weather variables[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2023, 39(5): 79-89. (in Chinese with English abstract) doi:10.11975/j.issn.1002-6819.202209108 http://www.tcsae.org

2022-09-14

2022-12-10

山西省基礎研究計劃項目(202103021224054,20210302124318);山西省高等學校科技創新項目(2021L025)

常曉敏,博士,副教授,碩士生導師,研究方向為智慧水利。Email:305643669@qq.com

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产精品成人一区二区不卡| a毛片免费观看| 狼友视频国产精品首页| 亚洲精品国产日韩无码AV永久免费网| 欧美成a人片在线观看| 2019国产在线| 日韩欧美在线观看| 97影院午夜在线观看视频| 97色伦色在线综合视频| 国产精品久久久久婷婷五月| 999在线免费视频| 亚洲欧洲免费视频| 一区二区三区四区精品视频| 欧美成人综合视频| 国产精品久久国产精麻豆99网站| 久久鸭综合久久国产| 国产成人一区在线播放| 亚洲国产精品无码久久一线| 国产区精品高清在线观看| 欧美精品成人一区二区视频一| 亚洲伊人天堂| 久久久久国产精品熟女影院| 黄色不卡视频| 欧美激情第一欧美在线| 四虎永久免费地址| 亚洲人成网站色7777| 免费在线不卡视频| 国产成人精品第一区二区| 亚洲视频无码| 国产丝袜无码一区二区视频| 精品成人一区二区三区电影| 美美女高清毛片视频免费观看| 无码精品国产dvd在线观看9久| 青青国产视频| 丰满的少妇人妻无码区| 国产一区在线观看无码| 色香蕉网站| 国产一区二区精品高清在线观看| 国产精品九九视频| 欧美激情一区二区三区成人| 亚洲日韩精品无码专区97| 99性视频| 福利小视频在线播放| 日本午夜影院| 日韩久草视频| 一本大道东京热无码av| 黄色污网站在线观看| 精品国产91爱| 欧洲精品视频在线观看| 婷婷亚洲最大| 久久免费精品琪琪| 久久人与动人物A级毛片| 一级黄色欧美| 手机在线国产精品| 凹凸国产熟女精品视频| 亚洲视频一区在线| 国产精品亚洲欧美日韩久久| 91在线无码精品秘九色APP| 欧美97欧美综合色伦图| 亚洲av色吊丝无码| 日韩一区二区三免费高清| 亚洲伊人天堂| 久久不卡国产精品无码| 美女高潮全身流白浆福利区| 毛片久久久| 国产极品美女在线观看| 国产在线视频福利资源站| 欧美成人影院亚洲综合图| 国产精品毛片一区视频播| 亚洲成人动漫在线观看| 久久国产精品嫖妓| 99re热精品视频中文字幕不卡| 亚洲美女高潮久久久久久久| 在线免费a视频| 午夜视频免费试看| 狠狠做深爱婷婷久久一区| 欧美色伊人| 国产自在线播放| 91麻豆精品国产高清在线| 超碰精品无码一区二区| 91偷拍一区| 精品视频在线观看你懂的一区|