• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    稻田土壤-作物系統(tǒng)模型參數(shù)敏感性分析與模型驗證

    2020-07-07 06:10:00史鑫蕊胡克林
    農(nóng)業(yè)機械學報 2020年5期
    關鍵詞:水力學硝化氮素

    史鑫蕊 梁 浩 周 豐 胡克林

    (1.中國農(nóng)業(yè)大學土地科學與技術學院,北京100193;2.農(nóng)業(yè)農(nóng)村部華北耕地保育重點實驗室,北京100193;3.河海大學農(nóng)業(yè)工程學院,南京210098;4.北京大學城市與環(huán)境學院,北京100871;5.教育部地表過程分析與模擬重點實驗室,北京100871)

    0 引言

    土壤-作物系統(tǒng)模型通過計算機程序和數(shù)學方程描述作物生長過程及碳氮元素在土壤-作物-大氣連續(xù)體中的運移和轉化規(guī)律[1]。與田間試驗相比,不僅可以節(jié)省大量的人力物力,而且可以通過情景分析模擬不同的田間管理、氣候變化、土壤屬性以及作物品種等因子對作物生長、水分消耗以及碳氮運移和轉化過程的影響,在水肥優(yōu)化管理、品種改良以及應對氣候變化等方面發(fā)揮了重要作用[2]。目前比較成熟的土壤-作物系統(tǒng)模型主要有EPIC、WOFOST、DNDC、RZWQM、APSIM 和DSSAT 等,在小麥、玉米、水稻等作物上已進行了廣泛的驗證和應用,但是這些模型應用于不同地區(qū)時均需要對模型參數(shù)進行校準和驗證[3]。由于輸入?yún)?shù)眾多,模型參數(shù)率定過程困難。另外,參數(shù)校驗過程中會產(chǎn)生“異參同效”現(xiàn)象,大大降低了模型校驗的可靠性[4-5]。敏感度分析能夠從眾多輸入?yún)?shù)中識別出敏感度高的參數(shù),通過簡化或固定低敏感參數(shù),可以減小模型校驗的工作量、降低不確定性[6-8]。

    模型敏感度分析主要有局部和全局敏感度分析兩種方法。局部敏感度法每次只改變一個參數(shù),評價單個參數(shù)變化對模擬結果的響應,忽略了參數(shù)間的交互作用對輸出結果的影響,僅適用于過程簡單的線性模型[9]。全局敏感度法綜合考慮各個參數(shù)及參數(shù)間的交互作用對輸出結果的影響,更適合于分析具有非線性的土壤-作物模型參數(shù)的敏感度[10]。常用的全局敏感度法主要有Morris、EFAST和Sobol'法等。DEJONGE 等[11]用Morris 和Sobol'兩種方法研究了不同灌溉條件下CERES-Maize 模型的參數(shù)敏感性,兩種方法得出的結果具有很強的相關性。宋明丹等[1]比較了Morris 和EFAST 兩種方法得出的CERES-Wheat 模型作物參數(shù)的敏感度,兩者結果基本一致,但Morris 方法計算工作量明顯較小。EFAST 和Sobol'均是基于方差分析的敏感度法,精度高且穩(wěn)定性強,但運算量大[12]。Morris法可以快速地從大量的模型輸入?yún)?shù)中識別出部分敏感度較高的參數(shù),運算量較小,已被廣泛應用于水文、生態(tài)和環(huán)境類模型的敏感度分析中[13-15]。俞雙恩等[16]采用Morris 法分析了DRAINMOD-S 模型6 個參數(shù)對土壤剖面含鹽量的的靈敏度。高穎會等[17]采用Morris 方法篩選出對SWMM 模型區(qū)域洪峰流量和徑流系數(shù)比較敏感的水文水力參數(shù)。JABLOUN 等[18]通過Morris 方法得到了Daisy 模型在模擬小麥-玉米輪作系統(tǒng)中作物產(chǎn)量和氮素淋洗的敏感參數(shù)。前人對各種模型的參數(shù)敏感度進行了大量研究,但大部分研究主要集中在旱地作物系統(tǒng)上,對于水田系統(tǒng)的研究還較少。

    LIANG 等[19]在參考和借鑒國外主流模型的基礎上開發(fā)了適用于水稻的WHCNS_Rice(Soil water heat carbon nitrogen simulator for rice)模型,該模型考慮了中國的氣候、土壤、環(huán)境和農(nóng)業(yè)管理特點,相比國外模型更適用于中國高強度集約化的管理模式,已成功用于模擬長江中游地區(qū)常規(guī)淹水和覆膜旱作條件下的水稻生長及稻田水氮遷移過程[20]。水田與旱地的環(huán)境條件及作物生長過程機理不同,水田考慮的過程更為復雜,輸入?yún)?shù)更多,而且具有很強的非線性,這就導致模型參數(shù)校準過程更為困難,大大限制了水田模型的推廣應用。此外,水田與旱地系統(tǒng)模型參數(shù)敏感性的異同點尚不明確。因此,本文以WHCNS_Rice 模型為例,應用長江中游地區(qū)2 年的水稻田間試驗數(shù)據(jù),基于Morris 和Sobol'兩種方法對該模型進行參數(shù)敏感度分析,旨在為提高水田系統(tǒng)模型的校準效率及其應用提供技術指導。

    1 材料與方法

    1.1 數(shù)據(jù)來源

    試驗地位于長江中游地區(qū)的湖北省荊州農(nóng)業(yè)氣象試驗站(30°21'N,112°9'E),種植制度為水稻-油菜輪作。研究區(qū)屬濕潤季風氣候,年平均氣溫16℃,年均降雨量1 095 mm。于2017 年和2018 年水稻生長季進行了田間試驗,對水稻生長、水分消耗及氮素去向等指標進行了動態(tài)觀測。水稻品種為“豐兩優(yōu)香一號”,2 年移栽日期分別為6 月6 日和6 月9 日,收獲日期分別為9 月14 日和9 月18 日,種植密度為9 萬株/hm2。小區(qū)面積150 m2(25 m×6 m),3次重復。根據(jù)當?shù)剞r(nóng)民管理習慣,施用純氮230.5 kg/hm2,純磷26.6 kg/hm2,純鉀40.4 kg/hm2。氮肥分3 次施用,分施比例為5∶3∶1,其中基肥和第2 次追肥施用復合肥,第1 次追肥施用尿素。水稻生長季內(nèi)水分管理采取淹水-曬田-再次淹水-乳熟期后自然落干至收獲[21]。

    試驗前開挖土壤剖面,取樣測定0 ~90 cm 土壤剖面各層的基本理化性質。試驗地土壤類型為水耕人為土,pH 值8.1,土壤有機碳質量比7.9 g/kg,全氮質量比1.2 g/kg。具體的土壤物理與水力學性質見表1。水稻生育期內(nèi)白天每2 h 測定一次田面水高度。通過流量計記錄每次發(fā)生灌水和徑流時的水量,取水樣采用流動分析儀(AA3 型,Bran +Luebbe,德國)測定銨態(tài)氮和硝態(tài)氮含量。2017 年采用動態(tài)箱法測定氨揮發(fā)量,水稻關鍵生長期取植株樣干燥測定秸稈和籽粒干物質質量,收獲時測產(chǎn),同時用碳氮元素分析儀(Flash2000 型,Thermo,美國)測定植株全氮含量。其他的測定指標及方法見文獻[21]。

    表1 土壤剖面基本物理與水力學性質Tab.1 Soil basic physical and hydraulic properties in soil profile

    1.2 WHCNS_Rice 模型

    WHCNS_Rice 模型以天為步長,由氣象數(shù)據(jù)和作物生物學參數(shù)驅動,可用于模擬水稻生長、稻田水分運動、氮素運移和轉化以及熱傳導等過程。模型采用Penman-Monteith 公式估算參考作物蒸散量。采用Green-Ampt 模型和Richards 方程分別模擬非飽和條件下土壤水分入滲和再分布過程,飽和條件下采用達西定律。每天的田面水高度通過水分平衡方程計算,公式為

    式中 Pn、Pn-1——第n 天、第n-1 天田面水高度,mm I

    n——第n 天灌溉量,mm

    P′n——第n 天降雨量,mm

    ETn——第n 天蒸散量,mm

    Rn——第n 天徑流量,mm

    Fn——第n 天水分入滲通量,mm

    稻田徑流的計算與旱地不同,稻田通常有一定高度的田埂,將此高度作為最大田面水允許高度。在強降雨條件下,需先將水稻田池蓄滿,超過田埂之后才會產(chǎn)生徑流。

    土壤有機質周轉以及碳氮循環(huán)直接參考Daisy模型。土壤氮素運移過程使用對流-彌散方程描述,源匯項考慮了有機質礦化、生物固持、尿素水解、氨揮發(fā)、硝化、反硝化和作物吸收等過程。采用改進的PS123 作物模型對作物生長發(fā)育進程、干物質生產(chǎn)和分配及作物產(chǎn)量形成過程進行模擬。詳細的模型原理見文獻[19 -20]。

    1.3 全局敏感性分析法

    1.3.1 Morris 方法

    Morris 是一種基于篩選分析的全局敏感度分析方法,利用OAT(One factor at a time)搜索途徑,通過計算輸入?yún)?shù)對輸出結果的基礎效應進行靈敏度分析。該方法有效平衡了計算的效率和準確性,對于輸入?yún)?shù)多且運算負荷大的模型具有很好的適用性。每個輸入?yún)?shù)對輸出結果的影響程度Ei[22]計算公式為

    其中

    式中 y——參數(shù)i 變化前的輸出值

    y*——參數(shù)i 變化后的輸出值

    Δi——參數(shù)i 的變化量Δ 介于1/(p -1)和1 -1/(p -1)之間,p 表示參數(shù)的水平。

    由于Morris 方法參數(shù)取值的隨機性,易造成取值誤差,所以需進行r 次重復,模型總運行次數(shù)為r(L+1)次,L 表示參數(shù)個數(shù)。根據(jù)Morris 方法建議,采樣因子r 設置為10,輸入?yún)?shù)個數(shù)L 為28,共運行模型290 次。最后計算每個參數(shù)對模擬值的影響程度Ei的平均值μ 和標準差σ。μ 越大,說明參數(shù)對輸出結果越敏感,而σ 表示的是參數(shù)間的交互作用,σ 越大表示該參數(shù)與其他參數(shù)交互作用越強。

    1.3.2 Sobol'方法

    Sobol'是一種基于方差的全局敏感性分析法,通過分解輸出變量的方差來定量化評價各輸入?yún)?shù)和參數(shù)間交互作用對輸出變量的影響程度[23]。假設y 為模型輸出變量;Xi為模型輸入?yún)?shù)(共m個)。y 的方差D(y)可分解為

    式中 Di——參數(shù)i 產(chǎn)生的方差

    Dij——參數(shù)i 和j 相互作用產(chǎn)生的方差

    Dijk——參數(shù)i、j 和k 相互作用產(chǎn)生的方差

    D1,2,…,m——m 個參數(shù)共同作用產(chǎn)生的方差

    對于參數(shù)i,可用一階敏感度指數(shù)Si反映單個參數(shù)的敏感性,全階敏感度指數(shù)STi表示單個參數(shù)與其他所有參數(shù)的共同影響,具體公式為

    式中 D~i——除了參數(shù)i 之外的參數(shù)的方差

    SALTELLI 等[23]建議,Sobol'方法取樣個數(shù)設置為n(m+2),其中m 為參數(shù)個數(shù)(28 個),n 取值應大于100。本研究計算了n 為1 024 時的Sobol'敏感度結果,需運行模型29 696 次。

    1.3.3 參數(shù)選擇、分布和范圍

    WHCNS_Rice 模型輸入?yún)?shù)主要分為土壤水力學參數(shù)、氮素轉化參數(shù)和作物參數(shù)3 類。首先根據(jù)實測的土壤理化性質數(shù)據(jù),基于土壤傳遞函數(shù)獲得水力學參數(shù)初始值[24],作物參數(shù)和氮素轉化參數(shù)初始值參考相關文獻數(shù)據(jù)或者直接采用默認值。結合文獻中其他模型敏感度分析結果[10],選擇28 個參數(shù)進行敏感性分析(表2),其中水力學參數(shù)每層5個共15 個、氮素轉化參數(shù)5 個、作物參數(shù)8 個。由于各參數(shù)的分布概率未知,假定所有的參數(shù)服從均勻分布。模型參數(shù)取值范圍采用初始值的±10%作為其區(qū)間的上限和下限,采用蒙特卡洛方法對每個參數(shù)隨機采樣。輸入?yún)?shù)組合的產(chǎn)生以及最后的結果分析均采用SimLab 2010 軟件完成[25]。

    表2 模型輸入?yún)?shù)Tab.2 Model input parameters

    1.4 模型參數(shù)率定

    根據(jù)模型敏感度分析結果,本文基于2017 年的大田實測數(shù)據(jù),采用試錯法對WHCNS_Rice 模型中的敏感參數(shù)進行率定。首先對土壤水力學參數(shù)進行調整,使模型模擬的田面水高度、蒸散量、徑流量等與實測值基本吻合。然后率定氮素轉化參數(shù),使氨揮發(fā)量、反硝化量、氮淋洗量、徑流損失量等模擬值與實測值吻合,若缺少實測值應使模擬值接近文獻參考值。最后,調整作物參數(shù),使產(chǎn)量、干物質質量、吸氮量等指標的模擬值與實測值盡量吻合。由于參數(shù)間存在交互作用,不同參數(shù)調整過程會相互影響,因此該過程需要周而復始,最終大部分模擬值與實測值誤差最小時模型率定過程結束。模型率定結束后,固定所有輸入?yún)?shù),根據(jù)2018 年實際田間管理和氣象數(shù)據(jù)運行模型,用實測的田面水高度、蒸散量、作物吸氮量和干物質質量等來驗證模型。若模型模擬值與實測值誤差在可接受范圍內(nèi),說明模型輸入?yún)?shù)合理,可用來模擬該地區(qū)水稻生長及稻田水分運動和氮素轉化過程。

    2 結果與分析

    2.1 Morris 敏感度分析

    2.1.1 水分輸出項對不同參數(shù)的敏感度

    根據(jù)模型輸入?yún)?shù)初始值,各參數(shù)設定±10%的變化范圍,基于Morris 敏感度分析方法,通過SimLab 軟件共產(chǎn)生了290 組模型組合參數(shù)。不同參數(shù)組合下分別運行WHCNS_Rice 模型,得到了290 組不同的輸出結果,最后采用SimLab 軟件進行參數(shù)敏感度分析。

    圖1 為模型模擬的蒸散量(ET)、水分滲漏量、徑流量和田面水高度對各輸入?yún)?shù)的敏感度結果。ET 對中期作物系數(shù)Kmid最敏感,其次為Kend和Kini,但是σ 較低,說明與其他參數(shù)交互作用較弱(圖1a)。Tsum是第4 敏感參數(shù),通過影響作物生育期天數(shù)間接影響ET,與其他參數(shù)有較強的交互作用。其余參數(shù)μ 均小于10 mm,但是F、Smax、θs等參數(shù)的σ 均較高,說明部分土壤水力學參數(shù)和作物參數(shù)有較強的交互作用。

    圖1 ET、滲漏量、田面水高度及徑流量對WHCNS_Rice 輸入?yún)?shù)的Morris 敏感度均值和標準差Fig.1 Average Morris mean effect and spread of WHCNS_Rice parameter with respect to ET,drainage,ponding water depth and runoff

    對于田面水高度、滲漏量和徑流量,敏感度最高的參數(shù)均為θs1,其μ 和σ 均較高(圖1)。其次為F1,具有較高的μ,但σ 相對θs1偏低。第3 敏感參數(shù)是Ks2,其μ 較小,但σ 較高。表層土壤飽和含水率和田間持水率很大程度上決定了土壤的蓄水能力,Ks2為犁底層土壤飽和導水率,決定了土壤水分的入滲速率。土壤蓄水率和導水率共同決定了土壤中水分的運動,從而對田面水高度、滲漏量和徑流量產(chǎn)生顯著影響?;谒制胶庠恚珽T、土壤儲水量、田面水高度、滲漏量和徑流量,這些水分平衡項之間存在相互影響、此消彼長的關系。土體儲水量的增加必然導致滲漏量降低,ET 會影響田面水高度的變化,進而對徑流量產(chǎn)生影響。

    2.1.2 氮素損失項對不同參數(shù)的敏感度

    WHCNS_Rice 模型稻田氮素損失項主要包括氮淋洗、反硝化、氨揮發(fā)和徑流損失等4 個氮素損失項。從模型敏感度分析結果來看,不同土層飽和含水率和田間持水率均對氮淋洗有較大的影響,氮淋洗對θs1和F1最敏感(圖2a)。其余敏感度較高的幾個參數(shù)按μ 值由大到小依次為θs2、θs3、F3、F2、W2,可以看出氮淋洗對各個土層的飽和含水率及田間持水率均比較敏感,這些參數(shù)值越大表示土壤蓄水能力越強,可減少水分滲漏,從而減小氮素隨水分下移造成的淋洗損失。其余參數(shù)敏感度均低于5 kg/hm2,對氮淋洗的影響較小。

    圖2 氮淋洗、反硝化、徑流和氨揮發(fā)對WHCNS_Rice 輸入?yún)?shù)的Morris 敏感度均值和標準差Fig.2 Average Morris mean effect and spread of WHCNS_Rice parameter with respect to nitrogen leaching,denitrification,runoff and ammonia volatilization

    對反硝化作用影響最大的參數(shù)是F1,但其σ 較低(圖2b)。其次是θs1和θs2,且具有較高的σ。反硝化作用是在厭氧條件下發(fā)生的,水力學參數(shù)通過影響土壤中的水分狀況進而對反硝化產(chǎn)生影響。反硝化經(jīng)驗系數(shù)Ad是第4 敏感參數(shù),與其他參數(shù)交互作用較弱,Ad越大,通過反硝化損失的氮素越多,對反硝化輸出項有最直接的影響。

    氮徑流損失對輸入?yún)?shù)的敏感性分析結果與氮淋洗相似,θs1和F1敏感度最大,其次是θs3和θs2(圖2c)。高蓄水能力下可以減少降雨和灌溉時產(chǎn)生的徑流損失。下層土壤的飽和含水率、田間持水率、飽和導水率等參數(shù)敏感度均值μ 相對較低,但非線性效應較強,主要通過影響水分的下滲速率對徑流產(chǎn)生間接作用。

    氨揮發(fā)對θs1和θs2最敏感,且具有較強的非線性效應,上層土壤的飽和含水率通過影響土壤中的銨態(tài)氮濃度間接影響氨揮發(fā)(圖2d)。對氨揮發(fā)影響較大的第3 個敏感參數(shù)是Kv,Kv越大,氨揮發(fā)速率越快。Smax是第4 敏感參數(shù),Smax決定了覆蓋度,覆蓋度影響了NH3從水面向大氣中的運動過程。不僅如此,覆蓋度是影響土面或水面溫度的重要原因,而氨揮發(fā)與溫度密切相關,高溫會加快氨逸出。

    2.1.3 作物生長指標對不同參數(shù)的敏感度

    從干物質質量、產(chǎn)量、吸氮量和葉面積指數(shù)(LAI)對模型輸入?yún)?shù)的敏感度結果來看,這些作物生長指標均對作物參數(shù)Tsum和Smax具有較高的敏感度(圖3)。Tsum控制了作物生長的物候期,Tsum變大作物生育期延長,反之生育期變短,Smax決定了作物冠層覆蓋度,因此這兩個作物參數(shù)對作物生長及產(chǎn)量形成均有較大影響(圖3a、3b)。土壤水力學參數(shù)θs1和θs2對干物質質量、產(chǎn)量和吸氮量影響較大,并與其他參數(shù)有較強的交互作用。由于水稻土一般具有犁底層,決定了水分入滲和養(yǎng)分淋失速率,因此Ks2對水稻生長也有一定的影響。土壤水力學參數(shù)通過影響稻田土壤中的水分狀況及氮素供應水平間接影響作物生長過程。

    圖3 干物質質量、產(chǎn)量、吸氮量和LAI 對WHCNS_Rice 輸入?yún)?shù)的Morris 敏感度均值和標準差Fig.3 Average Morris mean effect and spread of WHCNS_Rice parameter with respect to dry matter,yield,crop uptake and LAI

    對LAI 影響較大的3 個參數(shù)依次是Smax、Tsum、Smin,且σ 較高(圖3d)。Smax和Smin是與LAI 直接相關的參數(shù),直接影響LAI。Tsum通過控制生育期的長短對LAI 產(chǎn)生作用,三者共同決定了LAI。LAI 對Kmid的敏感度較低,但具有較高的σ,與其他參數(shù)交互作用較強。對干物質質量影響明顯的4 個參數(shù)依次為Tsum、Smax、θs1、θs2,其敏感度均值μ 均大于1 000 kg/hm2,θs2、Tsum和θs1的σ 較大。與前4 個參數(shù)相比,F(xiàn)1、Ks2、F3、Smin對干物質質量的影響相對較小,其μ 介于300 ~500 kg/hm2,干物質質量對其余參數(shù)不敏感。

    對作物產(chǎn)量影響較大的參數(shù)依次為θs1、θs2、Tsum、F1、Smax、Ks2、F3(圖3b)。產(chǎn)量對θs1和θs2敏感度最高,且交互作用強。從前面的結果可以看出,θs1和θs2對各項水分消耗及氮素損失均有較大的影響,進而會對作物生長期間的水分和氮素供應狀況產(chǎn)生作用,最終影響產(chǎn)量形成。作物吸氮量的敏感參數(shù)與產(chǎn)量基本相同,依次是θs1、Tsum、θs2、Smax、F1、Ks2、F3、W2。

    2.2 Sobol'敏感度分析

    Sobol'方法共產(chǎn)生了29 696 組模型參數(shù)組合,分別運行WHCNS_Rice 模型輸出相應的結果,然后用SimLab 2010 軟件進行全局敏感度分析。圖4 為28 個輸入?yún)?shù)對不同作物生長、水分及氮素平衡輸出項的一階及全階敏感度指數(shù)。一階敏感度指數(shù)劃定超過閾值0.1 的為敏感參數(shù)??梢钥闯觯州敵鲰椫姓羯⒘康拿舾袇?shù)分別為Kmid和θr1,滲漏量、徑流量及田面水高度的敏感參數(shù)均為θs1和F1。氮素損失項中,氮淋洗對F1敏感,氮徑流和氮反硝化對θs1和F1敏感,氨揮發(fā)對θr1和Ks1敏感。作物生長指標中,干物質質量、吸氮量和LAI 均對Tsum和Smax敏感,而產(chǎn)量只對Tsum敏感,Smax對其影響不大。整體上Sobol'一階敏感度指數(shù)得出的敏感參數(shù)與Morris 方法基本一致,但是氮素平衡項中氮淋洗和氨揮發(fā)無法獲得對其影響較大的參數(shù),因此改用手動調參,發(fā)現(xiàn)參數(shù)Kv和Ad分別對氨揮發(fā)和氮反硝化輸出結果有一定的影響,由此可見,Morris 得到的結果準確可靠。

    從全階敏感度指數(shù)看,大部分輸出項很難篩選出明顯的敏感參數(shù)。不同參數(shù)對蒸散量、滲漏量、徑流量和田面水高度的全階敏感度指數(shù)分別介于0.81 ~1.00、0.55 ~0.89、0.27 ~0.71、0.44 ~0.85之間,差異均不明顯。氮素損失項中,θs1和F1對氮徑流和氮反硝化的全階敏感度指數(shù)明顯高于其余參數(shù),而各參數(shù)對氨揮發(fā)和氮淋洗的全階敏感度指數(shù)差異不明顯。但作物生長輸出項的敏感參數(shù)比較明顯,除Tsum和Smax,θs1和θs2也有較高的敏感度,其次為F1和F2,這也與Morris 的結果一致(圖3)。

    某一參數(shù)全階敏感度指數(shù)與一階敏感度指數(shù)的差值即為該參數(shù)與其余參數(shù)交互作用對模型輸出項的影響程度。差值越大,表明該參數(shù)交互作用越強,這與Morris 方法中的σ 相對應,解釋了作物生長指標中部分土壤水力學參數(shù)敏感度較高的原因(圖3)。

    圖4 WHCNS_Rice 模型參數(shù)對不同輸出項的一階敏感度指數(shù)及全階敏感度指數(shù)Fig.4 First-order and total sensitivity of WHCNS_Rice model parameters for different objective functions

    與Morris 方法相比,Sobol'方法較好地考慮了參數(shù)間的交互作用,這是導致兩種方法結果略有不同的主要原因,但是Sobol'分析需要大批量運算才能得出有效的結果。另外,本文同時也計算了參數(shù)組合分別為928 和3 712 次時的敏感度結果,其敏感度指數(shù)出現(xiàn)大量的負值,無法獲得準確的敏感參數(shù)。此外,Sobol'產(chǎn)生的大量參數(shù)組合中個別參數(shù)在模型計算中無法輸出結果,導致結果為空。這種情況下,采用與之相近的上一組模型結果作為輸出值,這也造成了Sobol'分析結果的誤差。因此,綜合考慮計算量、可靠性及精度,Morris 不僅可以有效篩選出敏感參數(shù),而且計算量較小,是適合于稻田土壤-作物系統(tǒng)WHCNS_Rice 模型參數(shù)敏感度分析的有效方法。

    2.3 模型校準與驗證

    綜上模型敏感度分析結果,輸入?yún)?shù)的變化會對模擬結果造成很大的不確定性,合理的參數(shù)值對模擬結果的準確度十分重要。在模型調參過程中,應重點關注土壤水力學參數(shù)中的θs、F 和Ks,這幾個參數(shù)對水分運動、氮素運移和轉化以及作物生長均有顯著影響。其次是作物參數(shù),應重點關注Tsum、Smax和作物系數(shù)K。最后是氮素轉化參數(shù)中的Ad和Kv。以上參數(shù)是對模型敏感度較高的參數(shù),由于參數(shù)存在很強的非線性效應,且參數(shù)間存在交互作用,需要反復調整才能達到比較理想的效果。調參過程中對于敏感度較低的參數(shù),比如Rmax等,應根據(jù)實際情況進行設置,會使模型模擬效果更好。

    基于敏感度分析結果,使用2017 年實測數(shù)據(jù)校準模型。固定非敏感參數(shù),通過“試錯法”重點校準高敏感參數(shù),使模型模擬的干物質質量、吸氮量、田面水高度和蒸散量與實測值盡量吻合,最后采用2018 年實測數(shù)據(jù)驗證模型,校驗后的參數(shù)值見表2。

    從圖5 可以看出,模型校準年份(2017 年)的干物質質量、作物吸氮量、田面水高度和ET 的模擬值與實測值的相關系數(shù)r 分別為0.998、0.979、0.502和0.852,與2018 年模型驗證相對應的相關系數(shù)分別為0.998、0.985、0.879 和0.856。其中2017 年田面水高度模擬值與實測值的相關系數(shù)較低,主要是因為移栽后15 ~27 d 田面水高度模擬值和實測值相差較大所致。由于該時段并未發(fā)生大的降雨,只有第22、23 天有少量降雨,故有可能是人為觀測誤差造成的,但此誤差在可接受范圍內(nèi)。校準(驗證)的線性回歸方程的斜率分別為1.024 3、1.072 1、0.938 3 和0.975 0(1.000 1、1.120 5、1.099 9 和1.086 5),均接近1,表明模擬值與實測值顯著相關(P <0.01)(圖5)??傮w來說,模型模擬的各項指標與實測值之間具有很好的一致性,校驗后的WHCNS_Rice 模型可用于模擬該地區(qū)水稻生長過程、水分消耗以及稻田氮素的遷移和轉化過程。

    圖5 校準和驗證年份模擬和實測的干物質質量、作物吸氮量、田面水高度和ET 的相關性Fig.5 Measured/simulated values with linear regression line for dry matter,crop nitrogen uptake,ponding water depth and ET in 2017 and 2018

    3 討論

    土壤-作物系統(tǒng)模型是實現(xiàn)數(shù)字化農(nóng)田作物生長監(jiān)測、預測和決策支持的重要手段[12]。模型的參數(shù)校準是模型應用的重要前提和關鍵步驟[26],也是模型功能拓展的基礎,同時輸入?yún)?shù)的合理性和區(qū)域代表性決定了模擬結果的準確性。由于模型參數(shù)較多,在缺乏先驗知識的情況下,模型調參過程非常困難,限制了模型的推廣應用。敏感性分析可以從大量的參數(shù)中識別出敏感參數(shù),減少模型校驗的工作量和不確定性,進而提高模型校準的效率和精度。

    土壤-作物模型輸入?yún)?shù)大致可以分為3 大類:土壤水力學參數(shù)、作物參數(shù)和碳氮轉化參數(shù),不同類型的參數(shù)對模型輸出結果的影響有明顯差異。梁浩等[10]對WHCNS 模型小麥和玉米作物的參數(shù)敏感度進行了分析,結果表明土壤含水率、土壤硝態(tài)氮含量、作物產(chǎn)量和LAI 的綜合響應對作物參數(shù)最敏感,其次是土壤水力學參數(shù),而氮素轉化參數(shù)敏感度較低,這與本研究參數(shù)敏感度結果略有不同。本研究中除LAI 和干物質質量外,表現(xiàn)為土壤水力學參數(shù)敏感度高于作物參數(shù),尤其是水氮各輸出項。這可能是由水田和旱地土壤性質及環(huán)境條件的顯著差異造成的。水田與旱地土壤中的水氮運移有很大不同,土壤水力學參數(shù)是影響土壤中水分運動和氮素遷移的主要參數(shù)[27]。淹水稻田生育期內(nèi)很少出現(xiàn)水分脅迫情況,但是由于生育期內(nèi)頻繁降雨和大量灌溉,很容易發(fā)生大量的徑流和水分滲漏損失,其損失量就取決于土壤的水力學性質,尤其是犁底層飽和導水率,決定了稻田水分下滲速率。不僅如此,水分運動直接影響了土壤中氮素的遷移與轉化過程,從而會影響土壤的供氮水平。DEJONGE 等[11]的研究結果表明,在充分灌溉條件下,CERES-Maize 模型模擬的產(chǎn)量、蒸散量和LAI 對作物參數(shù)最敏感,而虧缺灌溉下田間持水率和萎蔫含水率是主要的敏感參數(shù),由此可見,土壤中的水分狀況是影響參數(shù)敏感性的重要因素。

    關于作物參數(shù)對不同模型輸出項的敏感度研究已經(jīng)有大量報道。ZHAO 等[28]研究了APSIM 模型輸出的產(chǎn)量、生物量、開花期和成熟期對作物參數(shù)的敏感度,發(fā)現(xiàn)產(chǎn)量和生物量均對積溫比較敏感。VANUYTRECHT 等[6]采用AquaCrop 模型模擬分析了4 個國家不同作物產(chǎn)量對參數(shù)的敏感度,結果表明水稻產(chǎn)量對積溫、最大覆蓋度敏感度較高。本研究中蒸散量、干物質質量、產(chǎn)量、吸氮量及LAI 對作物參數(shù)敏感度較高,特別是Tsum和Smax兩個參數(shù),這與上述研究結果一致。除Tsum和Smax外,本研究作物參數(shù)中作物系數(shù)是第3 敏感參數(shù),尤其是對ET有較大影響。由于WHCNS_Rice 模型中潛在作物蒸散量(ETp)是由參考作物蒸散量(ET0)與作物系數(shù)乘積得到的[29],因此不同時期作物系數(shù)的大小會對ET 產(chǎn)生直接影響。這與LIANG 等[30]的研究結果不太一致,其結果中對春玉米產(chǎn)量影響較大的第3 個參數(shù)是最大同化率Amax,而本研究中水稻生長對Amax不敏感,這主要是由于C3 和C4 作物光合作用機理不同造成的,在模型中采用了不同的光合作用公式。

    相對于土壤水力學參數(shù)和作物參數(shù),模型輸出項對氮素轉化參數(shù)的敏感度相對較低,只有氨揮發(fā)量和反硝化量分別對Kv和Ad敏感。而梁浩等[10]的研究結果表現(xiàn)為最大硝化速率Vn和硝化半飽和系數(shù)Kn敏感度較高,這與本研究的結果不一致。主要是由于水田和旱地的差異造成的,大量研究表明旱地硝化作用比較強烈,土壤硝態(tài)氮含量高,而銨態(tài)氮含量較低,因此與硝化作用相關的參數(shù)敏感度高。而水田土壤中硝態(tài)氮含量明顯低于銨態(tài)氮[31],硝化能力較弱,氨揮發(fā)和氮反硝化作用比較強烈,導致兩者的敏感參數(shù)正好相反。水田與旱地由于土壤水分條件不同,造成了氮素遷移和轉化過程差異較大,相應的敏感參數(shù)也有較大差異,尤其是氮素轉化參數(shù),因此模型重點校準參數(shù)也應根據(jù)土壤水分狀況調整。

    不同的作物品種、環(huán)境條件及田間管理均會對模型參數(shù)的敏感度產(chǎn)生較大影響。LIANG 等[30]研究表明不同水肥管理條件下模型參數(shù)敏感度有較大的變化。ASSENG 等[32]發(fā)現(xiàn)不同的氣候條件對產(chǎn)量模擬值造成了很大的不確定性。ZHAO 等[28]研究指出不同的施肥量比氣候條件和土壤類型的改變對模型參數(shù)敏感度的影響更大。本文僅研究了某一特定條件下的參數(shù)敏感度,沒有分析不同的土壤類型、氣候條件及水肥管理等差異造成的不確定性。為了擴大WHCNS_Rice 模型應用范圍,應進一步全面系統(tǒng)地考慮環(huán)境條件、作物品種、管理措施等不同情況下的參數(shù)敏感度及不確定性。

    4 結論

    (1)采用Morris 和Sobol'兩種全局敏感度分析方法,分別針對作物生長、水分運動及氮素損失3 部分輸出項,具體分析了稻田土壤-作物系統(tǒng)WHCNS_Rice 模型輸出項對土壤水力學參數(shù)、氮素轉化參數(shù)和作物參數(shù)的敏感性。結果表明,整體上Sobol'一階敏感度指數(shù)得出的敏感參數(shù)與Morris 方法基本一致,但是綜合考慮計算量及精度,Morris 方法可以快速并有效地篩選出模型敏感性參數(shù)。模型各個輸出項均對土壤水力學參數(shù)有較高的敏感度,尤其是對飽和含水率、田間持水率和飽和導水率的敏感度較高。其次是作物參數(shù),其中作物生育期總有效積溫、最大比葉面積和作物系數(shù)對產(chǎn)量、干物質質量、葉面積指數(shù)和蒸散量有較大影響,對水分運動和氮素遷移過程影響較小。相反,作物生長過程各輸出項對土壤水力學參數(shù)和作物參數(shù)均有較高的敏感性。各輸出項對氮素轉化參數(shù)敏感度相對最低,只有氨揮發(fā)一階動力學系數(shù)和反硝化經(jīng)驗系數(shù)分別對氨揮發(fā)和反硝化過程有一定的影響。

    (2)水田與旱地水分條件的不同,導致水田系統(tǒng)土壤水力學參數(shù)敏感度高于作物參數(shù),而旱地正好相反。另外,旱地系統(tǒng)由于硝化作用強烈,因此與硝化反應相關的參數(shù)(最大硝化速率Vn和硝化半飽和系數(shù)Kn)對氮素輸出項有較大影響,而在水田中與氨揮發(fā)和氮反硝化過程相關的參數(shù)對輸出結果影響較大。

    (3)基于模型參數(shù)敏感度篩選結果,根據(jù)稻田實測數(shù)據(jù)對WHCNS_Rice 模型參數(shù)進行了校驗。結果表明,模型模擬值與實測值具有較好的一致性,該模型可用于模擬長江中游地區(qū)水稻生長及稻田水氮遷移過程。該方法大大提高了模型校準效率,為土壤-作物系統(tǒng)模型的推廣應用提供了技術支持。

    猜你喜歡
    水力學硝化氮素
    飽和紫色土初始態(tài)和穩(wěn)定態(tài)細溝水力學特征研究*
    土壤學報(2022年1期)2022-03-08 08:52:10
    MBBR中進水有機負荷對短程硝化反硝化的影響
    二維水力學模型在紅光大橋洪水影響評價中的應用
    基于管網(wǎng)理論的人口遷移動力學模型構建
    科技視界(2016年27期)2017-03-14 23:09:34
    自排式沉沙池上游渠道水力學特性研究
    厭氧氨氧化與反硝化耦合脫氮除碳研究Ⅰ:
    海水反硝化和厭氧氨氧化速率同步測定的15N示蹤法及其應用
    楸樹無性系苗期氮素分配和氮素效率差異
    基于光譜分析的玉米氮素營養(yǎng)診斷
    氮素運籌對玉米干物質積累、氮素吸收分配及產(chǎn)量的影響
    免费观看av网站的网址| 嫩草影院精品99| 夫妻性生交免费视频一级片| 国产在线一区二区三区精| 午夜福利在线观看吧| 亚洲av男天堂| 精品一区二区三区人妻视频| 又黄又爽又刺激的免费视频.| 欧美xxxx黑人xx丫x性爽| 97超视频在线观看视频| 日本猛色少妇xxxxx猛交久久| 亚洲av男天堂| 欧美激情在线99| 免费黄色在线免费观看| 夜夜爽夜夜爽视频| 大陆偷拍与自拍| 亚洲精品日韩av片在线观看| 深爱激情五月婷婷| 97精品久久久久久久久久精品| 久久久久久久亚洲中文字幕| 午夜福利在线观看免费完整高清在| 久久国产乱子免费精品| 日日干狠狠操夜夜爽| 久久精品国产亚洲网站| 97超碰精品成人国产| 精品一区二区三区人妻视频| 日本猛色少妇xxxxx猛交久久| 禁无遮挡网站| 夜夜爽夜夜爽视频| 在线观看av片永久免费下载| 九色成人免费人妻av| 一区二区三区乱码不卡18| 亚洲精品国产av蜜桃| 亚洲av中文字字幕乱码综合| 久久久久免费精品人妻一区二区| 中文字幕免费在线视频6| 一二三四中文在线观看免费高清| 狂野欧美激情性xxxx在线观看| 国产精品女同一区二区软件| 一区二区三区四区激情视频| 亚洲综合精品二区| 国产片特级美女逼逼视频| 免费黄色在线免费观看| 一级av片app| av黄色大香蕉| 午夜免费观看性视频| 亚洲国产日韩欧美精品在线观看| 久久久久国产网址| 日本免费在线观看一区| 日韩av不卡免费在线播放| 日本午夜av视频| 国产 亚洲一区二区三区 | 亚洲成人中文字幕在线播放| 日日撸夜夜添| 精品一区在线观看国产| 国产精品久久久久久av不卡| 水蜜桃什么品种好| 如何舔出高潮| 亚洲人成网站高清观看| 免费看不卡的av| 男女啪啪激烈高潮av片| 丰满乱子伦码专区| 一区二区三区免费毛片| 99热这里只有精品一区| 欧美最新免费一区二区三区| 亚洲av免费在线观看| 午夜免费激情av| 免费黄网站久久成人精品| 麻豆成人av视频| av网站免费在线观看视频 | 丰满人妻一区二区三区视频av| 久久久久久久久大av| 国产一区亚洲一区在线观看| 欧美性猛交╳xxx乱大交人| 观看美女的网站| 丝袜美腿在线中文| 中文乱码字字幕精品一区二区三区 | 亚洲欧美中文字幕日韩二区| 蜜桃亚洲精品一区二区三区| 久久这里有精品视频免费| 国产伦一二天堂av在线观看| 亚洲av日韩在线播放| 亚洲真实伦在线观看| 免费播放大片免费观看视频在线观看| 淫秽高清视频在线观看| 日韩一区二区三区影片| 天天一区二区日本电影三级| 99久久中文字幕三级久久日本| 午夜老司机福利剧场| 最近最新中文字幕大全电影3| 你懂的网址亚洲精品在线观看| av免费在线看不卡| 2021少妇久久久久久久久久久| 最近最新中文字幕大全电影3| 国产成人精品久久久久久| 中文精品一卡2卡3卡4更新| 插阴视频在线观看视频| 可以在线观看毛片的网站| 男女啪啪激烈高潮av片| 深爱激情五月婷婷| 日本三级黄在线观看| 精品不卡国产一区二区三区| 国产国拍精品亚洲av在线观看| 亚洲国产色片| 欧美+日韩+精品| 日本午夜av视频| 国产成人a∨麻豆精品| 高清毛片免费看| 99视频精品全部免费 在线| 免费电影在线观看免费观看| 久久久久精品性色| 欧美人与善性xxx| 婷婷色综合大香蕉| 97热精品久久久久久| 国产亚洲av片在线观看秒播厂 | 欧美性感艳星| 九九爱精品视频在线观看| 三级男女做爰猛烈吃奶摸视频| 日韩大片免费观看网站| 人妻系列 视频| 免费看光身美女| av国产免费在线观看| 精品久久久久久久末码| 国产探花在线观看一区二区| 亚洲精品久久久久久婷婷小说| 日本av手机在线免费观看| 国产精品久久久久久久久免| 国产有黄有色有爽视频| 夫妻午夜视频| 亚洲成人一二三区av| 人妻制服诱惑在线中文字幕| 777米奇影视久久| 国产成人aa在线观看| av在线亚洲专区| 亚洲成人中文字幕在线播放| 男插女下体视频免费在线播放| 日韩电影二区| 五月伊人婷婷丁香| 亚洲最大成人手机在线| 亚洲欧美成人精品一区二区| 男人舔女人下体高潮全视频| 成人av在线播放网站| 国产激情偷乱视频一区二区| 中文欧美无线码| 国产高清有码在线观看视频| 国产免费又黄又爽又色| 黄片无遮挡物在线观看| 午夜福利视频精品| 欧美变态另类bdsm刘玥| 日本色播在线视频| 久久久色成人| 中国美白少妇内射xxxbb| 黄片wwwwww| 精品久久久久久久末码| 国产午夜精品一二区理论片| 久久精品夜夜夜夜夜久久蜜豆| 日本一二三区视频观看| a级一级毛片免费在线观看| 国产爱豆传媒在线观看| 久久午夜福利片| 国产日韩欧美在线精品| 亚洲自拍偷在线| 麻豆国产97在线/欧美| 亚洲av免费高清在线观看| 看黄色毛片网站| 亚洲熟女精品中文字幕| 国产精品嫩草影院av在线观看| 国产午夜精品一二区理论片| 免费av毛片视频| 能在线免费观看的黄片| 99re6热这里在线精品视频| 亚洲国产精品专区欧美| 亚洲成人久久爱视频| 小蜜桃在线观看免费完整版高清| 国产爱豆传媒在线观看| 亚洲国产精品sss在线观看| 国产日韩欧美在线精品| 国产成人精品福利久久| ponron亚洲| 亚洲av男天堂| 91精品一卡2卡3卡4卡| av天堂中文字幕网| 午夜福利高清视频| 我的女老师完整版在线观看| 日日摸夜夜添夜夜爱| 亚洲一区高清亚洲精品| 午夜免费观看性视频| 永久免费av网站大全| 国产欧美另类精品又又久久亚洲欧美| videossex国产| 国产乱人偷精品视频| 免费无遮挡裸体视频| 哪个播放器可以免费观看大片| 欧美一区二区亚洲| 国产亚洲一区二区精品| 久久久精品免费免费高清| 免费观看在线日韩| 97超视频在线观看视频| 亚洲欧洲日产国产| 国产中年淑女户外野战色| 国产又色又爽无遮挡免| 看非洲黑人一级黄片| 大香蕉97超碰在线| 亚洲不卡免费看| 亚洲经典国产精华液单| 国产一区二区亚洲精品在线观看| 亚洲成人久久爱视频| 亚洲精华国产精华液的使用体验| 亚洲一级一片aⅴ在线观看| 午夜激情欧美在线| 午夜福利在线观看免费完整高清在| 深爱激情五月婷婷| 免费少妇av软件| 91av网一区二区| 久久国内精品自在自线图片| 人妻一区二区av| 七月丁香在线播放| 亚洲国产成人一精品久久久| 日韩av在线免费看完整版不卡| 91精品伊人久久大香线蕉| 免费观看精品视频网站| 日本wwww免费看| 老司机影院成人| 国产探花极品一区二区| av又黄又爽大尺度在线免费看| 啦啦啦啦在线视频资源| av.在线天堂| 午夜亚洲福利在线播放| 久久久精品免费免费高清| 亚洲av成人精品一二三区| 亚洲欧美一区二区三区国产| 你懂的网址亚洲精品在线观看| 偷拍熟女少妇极品色| 一级爰片在线观看| 26uuu在线亚洲综合色| 简卡轻食公司| 你懂的网址亚洲精品在线观看| 偷拍熟女少妇极品色| 久久久久久久亚洲中文字幕| 91av网一区二区| 久久久久网色| 国产有黄有色有爽视频| 一个人观看的视频www高清免费观看| 午夜免费观看性视频| 亚洲怡红院男人天堂| 国产精品美女特级片免费视频播放器| 一级毛片黄色毛片免费观看视频| 91久久精品国产一区二区成人| av专区在线播放| 日韩制服骚丝袜av| 久久草成人影院| 永久免费av网站大全| 国产免费视频播放在线视频 | 99久久精品热视频| 久久99热6这里只有精品| 国产午夜福利久久久久久| 国产黄片视频在线免费观看| 免费av观看视频| 日韩三级伦理在线观看| 欧美日韩亚洲高清精品| 男人舔奶头视频| 久久久久九九精品影院| 成年版毛片免费区| 搡老妇女老女人老熟妇| 国产亚洲精品av在线| xxx大片免费视频| 99热这里只有精品一区| 只有这里有精品99| 亚洲av免费高清在线观看| 欧美不卡视频在线免费观看| 亚洲精品成人av观看孕妇| 亚洲欧美成人精品一区二区| 高清av免费在线| 伊人久久国产一区二区| 亚洲自拍偷在线| 能在线免费观看的黄片| 成人毛片60女人毛片免费| 亚洲婷婷狠狠爱综合网| 99re6热这里在线精品视频| 有码 亚洲区| 日日啪夜夜爽| 日韩精品有码人妻一区| 如何舔出高潮| 91av网一区二区| 免费黄色在线免费观看| 久久热精品热| 亚洲精品日本国产第一区| av福利片在线观看| 卡戴珊不雅视频在线播放| 国产成人午夜福利电影在线观看| 亚洲av免费在线观看| 熟妇人妻不卡中文字幕| 国产不卡一卡二| 国产乱来视频区| 日本-黄色视频高清免费观看| 女人十人毛片免费观看3o分钟| kizo精华| 男女边吃奶边做爰视频| 美女被艹到高潮喷水动态| 国产精品一及| 99久久人妻综合| 日韩欧美精品v在线| 最近的中文字幕免费完整| 免费观看在线日韩| 不卡av一区二区三区| 亚洲成色77777| 天美传媒精品一区二区| 国产成人午夜福利电影在线观看| 国产片内射在线| 男人添女人高潮全过程视频| 中文精品一卡2卡3卡4更新| 极品少妇高潮喷水抽搐| 99热网站在线观看| 亚洲国产精品国产精品| 精品一区二区三区四区五区乱码 | 天堂8中文在线网| 久热这里只有精品99| 精品午夜福利在线看| av线在线观看网站| 一区二区av电影网| 亚洲av在线观看美女高潮| 一区在线观看完整版| 国产精品国产三级国产专区5o| 免费女性裸体啪啪无遮挡网站| 欧美精品国产亚洲| 精品国产超薄肉色丝袜足j| 国产精品二区激情视频| 亚洲伊人久久精品综合| 夜夜骑夜夜射夜夜干| 午夜免费男女啪啪视频观看| 少妇被粗大猛烈的视频| 国产成人aa在线观看| 最近最新中文字幕大全免费视频 | 最近中文字幕高清免费大全6| 在线观看www视频免费| 满18在线观看网站| 女性被躁到高潮视频| 熟女电影av网| 咕卡用的链子| 肉色欧美久久久久久久蜜桃| 日韩欧美精品免费久久| 在线天堂最新版资源| 精品国产国语对白av| 性高湖久久久久久久久免费观看| 亚洲,一卡二卡三卡| 国产精品二区激情视频| 日本爱情动作片www.在线观看| 免费在线观看完整版高清| 欧美精品av麻豆av| 18禁国产床啪视频网站| 国产精品女同一区二区软件| 99国产精品免费福利视频| 成人漫画全彩无遮挡| 成年动漫av网址| 欧美精品国产亚洲| 午夜免费观看性视频| 高清av免费在线| 亚洲三级黄色毛片| 97在线视频观看| 午夜福利乱码中文字幕| 丰满乱子伦码专区| 久久亚洲国产成人精品v| 久久ye,这里只有精品| 人妻 亚洲 视频| 欧美精品国产亚洲| 美女中出高潮动态图| 国产成人精品婷婷| 美女中出高潮动态图| 国产成人精品婷婷| 秋霞伦理黄片| 国产视频首页在线观看| 国产毛片在线视频| 欧美日韩视频高清一区二区三区二| 亚洲男人天堂网一区| 久久综合国产亚洲精品| 99re6热这里在线精品视频| 久久久久久久久免费视频了| 日韩,欧美,国产一区二区三区| 久久人人爽av亚洲精品天堂| 亚洲国产日韩一区二区| av卡一久久| 国产白丝娇喘喷水9色精品| 啦啦啦视频在线资源免费观看| 欧美激情极品国产一区二区三区| 亚洲国产欧美网| 国产高清不卡午夜福利| 亚洲欧美清纯卡通| 高清黄色对白视频在线免费看| 亚洲精品日本国产第一区| 国产福利在线免费观看视频| 国产精品99久久99久久久不卡 | 男女国产视频网站| a级片在线免费高清观看视频| 亚洲国产精品一区三区| 欧美日本中文国产一区发布| 热re99久久国产66热| 国产白丝娇喘喷水9色精品| 丝袜美足系列| 免费观看av网站的网址| 亚洲一区中文字幕在线| 电影成人av| 秋霞在线观看毛片| 国产精品一国产av| 999久久久国产精品视频| 久久久精品国产亚洲av高清涩受| 精品亚洲乱码少妇综合久久| 亚洲国产欧美在线一区| 国产午夜精品一二区理论片| 观看美女的网站| 国产欧美亚洲国产| 日韩欧美精品免费久久| 国产成人精品久久二区二区91 | 999精品在线视频| 最黄视频免费看| 欧美成人午夜免费资源| 99久久中文字幕三级久久日本| 免费黄频网站在线观看国产| 欧美中文综合在线视频| 18禁国产床啪视频网站| 亚洲,一卡二卡三卡| 亚洲成色77777| 一本久久精品| 天天躁夜夜躁狠狠久久av| 少妇 在线观看| 91国产中文字幕| 18禁动态无遮挡网站| 国产麻豆69| 涩涩av久久男人的天堂| 2018国产大陆天天弄谢| 蜜桃在线观看..| 天美传媒精品一区二区| 久久国产精品大桥未久av| 免费在线观看视频国产中文字幕亚洲 | 伦理电影免费视频| 国产一区二区三区综合在线观看| 色播在线永久视频| 大话2 男鬼变身卡| 韩国av在线不卡| 男人操女人黄网站| 极品少妇高潮喷水抽搐| 色哟哟·www| 国产伦理片在线播放av一区| 啦啦啦视频在线资源免费观看| 色网站视频免费| 超碰成人久久| 亚洲精品久久久久久婷婷小说| 亚洲精品中文字幕在线视频| 国产成人91sexporn| 国产免费又黄又爽又色| 久久久久人妻精品一区果冻| 在线天堂中文资源库| 亚洲欧美一区二区三区久久| 纵有疾风起免费观看全集完整版| 国产亚洲一区二区精品| 日本wwww免费看| 免费少妇av软件| 欧美精品av麻豆av| 伦理电影大哥的女人| 一区二区三区乱码不卡18| 交换朋友夫妻互换小说| 婷婷色av中文字幕| 精品一区二区三区四区五区乱码 | 免费观看无遮挡的男女| 精品酒店卫生间| 国产在视频线精品| 男女边摸边吃奶| 久热久热在线精品观看| 18禁裸乳无遮挡动漫免费视频| 久久精品aⅴ一区二区三区四区 | 国产精品国产av在线观看| xxxhd国产人妻xxx| 国产熟女欧美一区二区| 少妇人妻精品综合一区二区| 街头女战士在线观看网站| 中文字幕亚洲精品专区| 少妇的丰满在线观看| 午夜福利一区二区在线看| 午夜免费鲁丝| 国产有黄有色有爽视频| 精品少妇一区二区三区视频日本电影 | 国产女主播在线喷水免费视频网站| 国产日韩欧美视频二区| 欧美日韩精品成人综合77777| videos熟女内射| 色哟哟·www| 十八禁网站网址无遮挡| 午夜老司机福利剧场| 天天躁夜夜躁狠狠久久av| 国产一级毛片在线| 一本—道久久a久久精品蜜桃钙片| h视频一区二区三区| 日韩成人av中文字幕在线观看| 久久国产精品大桥未久av| 9热在线视频观看99| 两性夫妻黄色片| av在线app专区| 丰满乱子伦码专区| 看十八女毛片水多多多| 人妻少妇偷人精品九色| 日韩中字成人| 亚洲国产看品久久| 日本av手机在线免费观看| 国产有黄有色有爽视频| 麻豆精品久久久久久蜜桃| 亚洲欧美成人综合另类久久久| 国产日韩欧美在线精品| 国产综合精华液| 国产av精品麻豆| 亚洲五月色婷婷综合| 日韩不卡一区二区三区视频在线| 一区福利在线观看| 黄色毛片三级朝国网站| 欧美精品一区二区大全| 在线亚洲精品国产二区图片欧美| av网站在线播放免费| 国产探花极品一区二区| 欧美国产精品va在线观看不卡| 国产精品成人在线| 女人高潮潮喷娇喘18禁视频| 国产高清国产精品国产三级| 亚洲成人av在线免费| 国产 精品1| 最黄视频免费看| 亚洲三级黄色毛片| 国产av国产精品国产| 永久网站在线| 婷婷色av中文字幕| 亚洲成国产人片在线观看| 亚洲国产精品999| 性色avwww在线观看| 国产极品天堂在线| 十八禁高潮呻吟视频| 国产人伦9x9x在线观看 | 欧美精品亚洲一区二区| 视频区图区小说| 久久精品aⅴ一区二区三区四区 | 一本一本久久a久久精品综合妖精 国产伦在线观看视频一区 | 国产又色又爽无遮挡免| 国产免费一区二区三区四区乱码| 这个男人来自地球电影免费观看 | 免费高清在线观看日韩| 久久人人爽人人片av| 日本wwww免费看| 亚洲人成网站在线观看播放| 王馨瑶露胸无遮挡在线观看| 18在线观看网站| 国产日韩欧美亚洲二区| 亚洲av电影在线观看一区二区三区| 最近中文字幕2019免费版| 日韩一区二区三区影片| 中国三级夫妇交换| 国产精品二区激情视频| 色婷婷久久久亚洲欧美| 精品国产乱码久久久久久小说| 久久婷婷青草| 色播在线永久视频| 欧美日韩成人在线一区二区| 丝袜喷水一区| 校园人妻丝袜中文字幕| 久久99一区二区三区| 三上悠亚av全集在线观看| 色哟哟·www| 国产成人精品在线电影| 韩国av在线不卡| 成人国语在线视频| 日本黄色日本黄色录像| 欧美日韩精品网址| 午夜av观看不卡| 国产成人精品一,二区| 亚洲 欧美一区二区三区| 男的添女的下面高潮视频| 伊人久久大香线蕉亚洲五| 伦理电影免费视频| 精品国产乱码久久久久久男人| 韩国高清视频一区二区三区| 精品午夜福利在线看| 岛国毛片在线播放| 国产精品99久久99久久久不卡 | 一本久久精品| 亚洲熟女精品中文字幕| 亚洲精品久久成人aⅴ小说| 亚洲国产欧美网| 久久精品亚洲av国产电影网| 免费少妇av软件| 久久国产亚洲av麻豆专区| 婷婷色麻豆天堂久久| 欧美av亚洲av综合av国产av | 久久久久久人妻| 国产精品一区二区在线不卡| 你懂的网址亚洲精品在线观看| 国产成人欧美| 黑人猛操日本美女一级片| 国产精品av久久久久免费| 老汉色∧v一级毛片| 免费久久久久久久精品成人欧美视频| 波野结衣二区三区在线| 一区二区三区激情视频| 久久久久视频综合| 女性生殖器流出的白浆| 三级国产精品片| 自拍欧美九色日韩亚洲蝌蚪91| 亚洲国产日韩一区二区| 免费黄色在线免费观看| 可以免费在线观看a视频的电影网站 | 99精国产麻豆久久婷婷| 欧美亚洲日本最大视频资源| 国产亚洲午夜精品一区二区久久| 最近中文字幕2019免费版| 丝袜在线中文字幕| 男女国产视频网站| 欧美xxⅹ黑人| 国产深夜福利视频在线观看| 国产综合精华液| 久久久亚洲精品成人影院| 欧美日韩一区二区视频在线观看视频在线|