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

    采用AquaCrop作物生長(zhǎng)模型研究中國(guó)玉米干旱脆弱性

    2020-03-03 11:45:44朱秀芳侯陳瑤
    關(guān)鍵詞:種植區(qū)脆弱性校正

    徐 昆,朱秀芳,2,3,劉 瑩,侯陳瑤

    ·農(nóng)業(yè)信息與電氣技術(shù)·

    采用AquaCrop作物生長(zhǎng)模型研究中國(guó)玉米干旱脆弱性

    徐 昆1,朱秀芳1,2,3※,劉 瑩1,侯陳瑤1

    (1. 北京師范大學(xué)地理科學(xué)學(xué)部遙感科學(xué)與工程研究院,北京 100875;2. 北京師范大學(xué)地表過(guò)程與資源生態(tài)學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京 100875;3. 北京師范大學(xué)環(huán)境演變與自然災(zāi)害教育部重點(diǎn)實(shí)驗(yàn)室,北京 100875)

    干旱脆弱性評(píng)價(jià)作為干旱風(fēng)險(xiǎn)評(píng)估和災(zāi)損評(píng)估的重要環(huán)節(jié),在保障國(guó)家糧食安全和農(nóng)業(yè)可持續(xù)發(fā)展中具有重大意義。該文以中國(guó)5大玉米種植區(qū)為研究區(qū)域,以其中241個(gè)主要玉米種植城市為基本單元,采用擴(kuò)展傅里葉幅度檢驗(yàn)法選取出2個(gè)敏感參數(shù)(作物冠層形成后到衰老之前的作物系數(shù)和參考收獲指數(shù)),并在此基礎(chǔ)上對(duì)AquaCrop作物模型進(jìn)行逐市的參數(shù)標(biāo)定。利用參數(shù)標(biāo)定后的模型對(duì)不同灌溉條件下玉米受到的水分脅迫及相應(yīng)情景下的產(chǎn)量進(jìn)行模擬計(jì)算,分別建立了5個(gè)玉米種植區(qū)對(duì)應(yīng)的干旱脆弱性曲線。結(jié)果表明:5個(gè)區(qū)域的脆弱性曲線擬合結(jié)果均為S形曲線,當(dāng)干旱強(qiáng)度指標(biāo)達(dá)到0.2附近時(shí),產(chǎn)量損失率開始迅速增加;當(dāng)干旱強(qiáng)度指標(biāo)達(dá)到0.6左右時(shí),產(chǎn)量損失率接近最大值。擬合函數(shù)的決定系數(shù)2分別在0.47~0.98之間,曲線擬合結(jié)果較好,在中國(guó)區(qū)域性玉米干旱脆弱性研究與干旱風(fēng)險(xiǎn)評(píng)估領(lǐng)域具有一定的理論與應(yīng)用價(jià)值。

    干旱;水分脅迫;脆弱性曲線;AquaCrop模型;玉米

    0 引 言

    自20世紀(jì)50年代以來(lái),以變暖為主要特征的全球氣候變化不斷加劇,干旱、洪澇等極端氣象事件發(fā)生的頻率也隨之增加。農(nóng)業(yè)生產(chǎn)在很大程度上依賴于自然氣候條件,極易在氣象災(zāi)害中受到損失[1-2]。聯(lián)合國(guó)政府間氣候變化專門委員會(huì)(intergovernmental panel on climate change,IPCC)和聯(lián)合國(guó)糧農(nóng)組織(food and agriculture organization of the united nations,F(xiàn)AO)均在報(bào)告中指出,農(nóng)業(yè)是最容易受到氣候變化影響的產(chǎn)業(yè)之一[3-5]。而干旱作為因氣候變化引起的常見(jiàn)自然災(zāi)害之一,因其持續(xù)時(shí)間長(zhǎng)、影響復(fù)雜的特性,被認(rèn)為是最具破壞性的自然災(zāi)害之一。過(guò)去幾十年,中國(guó)干旱事件頻發(fā)[6-9]。2001-2016年,中國(guó)年均耕地干旱受災(zāi)面積約1.85×107hm2,占各類氣象災(zāi)害受災(zāi)總面積的近50%,其中絕收面積達(dá)到2.16×106hm2[10]。農(nóng)業(yè)干旱已成為威脅中國(guó)糧食安全和農(nóng)業(yè)可持續(xù)發(fā)展的重要因素[11-12]。因此,干旱風(fēng)險(xiǎn)評(píng)估在農(nóng)業(yè)管理中具有重要意義。

    作為對(duì)承災(zāi)體抵抗災(zāi)害能力的一種度量,脆弱性評(píng)估是風(fēng)險(xiǎn)評(píng)估的重要基礎(chǔ)。早期的脆弱性評(píng)估以定性研究為主[13-14],后來(lái)隨著模糊數(shù)學(xué)在風(fēng)險(xiǎn)評(píng)價(jià)中的應(yīng)用,脆弱性評(píng)估的定量方法得到不斷發(fā)展。目前,常用的定量脆弱性評(píng)價(jià)方法主要包含3種類型[15-16]:基于歷史災(zāi)害損失數(shù)據(jù)的脆弱性評(píng)價(jià)方法[17]、基于指標(biāo)的脆弱性評(píng)價(jià)方法[18]以及基于災(zāi)害損失曲線(脆弱性曲線)的脆弱性評(píng)價(jià)方法。

    脆弱性曲線已廣泛應(yīng)用于洪水[19-20]、地震[21-23]、臺(tái)風(fēng)[24-25]、滑坡[26-27]、雪崩[28-29]、冰雹[30]等多類災(zāi)種的脆弱性定量研究中。在農(nóng)業(yè)干旱的研究領(lǐng)域中,隨著多種作物生長(zhǎng)模型的迅速發(fā)展,模型構(gòu)建法成為脆弱性曲線研究的重要方法之一[31]。Yin等[32]基于EPIC作物生長(zhǎng)模型,利用水分脅迫指數(shù)構(gòu)建了全球35個(gè)國(guó)家和地區(qū)的玉米干旱脆弱性曲線,Guo等[16]在此基礎(chǔ)上加入了環(huán)境因素的考慮,建立了“產(chǎn)量損失-干旱指數(shù)-環(huán)境指標(biāo)”的三維脆弱性曲面。Wang等[15]、Yue等[33]分別基于EPIC模型構(gòu)建了中國(guó)小麥的區(qū)域性脆弱性曲線,Jia等[34-35]則分別對(duì)中國(guó)不同區(qū)域的玉米脆弱性曲線進(jìn)行擬合研究。

    上述研究對(duì)作物生長(zhǎng)模型的校正多采用國(guó)家、省份、流域、自然區(qū)劃等較大區(qū)域?yàn)槟P托U幕締卧?,難以保證單元內(nèi)部的自然地理?xiàng)l件的均一性。干旱指標(biāo)的構(gòu)建采用水分脅迫指數(shù)的逐日累加的方法,只考慮了干旱的累積效應(yīng)而忽略了干旱持續(xù)時(shí)間的影響,為解決以上問(wèn)題,該文以中國(guó)5大玉米種植區(qū)域內(nèi)的241個(gè)玉米種植城市為研究區(qū)域,以地級(jí)市為模型校正的基本單元,對(duì)AquaCrop作物模型參數(shù)進(jìn)行逐市的標(biāo)定,在此基礎(chǔ)上模擬了不同灌溉情景下玉米受到的水分脅迫及對(duì)應(yīng)的產(chǎn)量,使用水分脅迫的日均值代替?zhèn)鹘y(tǒng)的累加值來(lái)描述干旱強(qiáng)度,進(jìn)而構(gòu)建5個(gè)玉米種植區(qū)的干旱脆弱性曲線。

    1 研究區(qū)概況

    玉米是中國(guó)三大主要谷類作物之一。中國(guó)玉米種植范圍廣泛,主要集中在東北、華北和西南地區(qū),形成從東北到西南的一條長(zhǎng)而傾斜的玉米種植帶。根據(jù)氣候、土壤、地貌等地理?xiàng)l件以及耕作制度等各項(xiàng)因素,佟屏亞[36]將中國(guó)劃分為6個(gè)主要玉米種植區(qū),分別為北方春播玉米種植區(qū)、黃淮海夏播玉米種植區(qū)、西南山地玉米種植區(qū)、南方丘陵玉米種植區(qū)、西北灌溉玉米種植區(qū)以及青藏高原玉米種植區(qū)。

    該文在Earth Stat全球作物分布數(shù)據(jù)的基礎(chǔ)上,結(jié)合各省市、農(nóng)村的統(tǒng)計(jì)年鑒與調(diào)查年鑒,最終選取中國(guó)241個(gè)主要玉米種植城市作為研究對(duì)象。上述241個(gè)城市均分布在5個(gè)玉米種植區(qū)中:73個(gè)在北方春播玉米種植區(qū),66個(gè)在黃淮海夏播玉米種植區(qū),52個(gè)在西南山地玉米種植區(qū),34個(gè)在南方丘陵玉米種植區(qū),16個(gè)在西北灌溉玉米種植區(qū),如圖1所示。

    圖1 中國(guó)玉米產(chǎn)區(qū)分布

    2 數(shù)據(jù)與方法

    2.1 AquaCrop模型原理簡(jiǎn)介

    AquaCrop是2009年由FAO組織研發(fā)出的一款新型作物生長(zhǎng)模型,模型主要包含土壤水分平衡、作物生長(zhǎng)模擬和大氣組分3個(gè)基本模塊[37]。它是水分驅(qū)動(dòng)模型[38],通過(guò)控制土壤中可利用水的含量來(lái)影響作物產(chǎn)量。FAO灌溉與排水第33號(hào)文件給出了作物產(chǎn)量和水分響應(yīng)的轉(zhuǎn)換關(guān)系,如式(1)所示:

    式中Y和0分別為作物的潛在產(chǎn)量和實(shí)際產(chǎn)量,kg/m2;ET和ET0分別為作物潛在蒸散量和實(shí)際蒸散量,mm;k為產(chǎn)量對(duì)水分響應(yīng)的系數(shù)。

    AquaCrop模型對(duì)上述方程進(jìn)行了改進(jìn),將蒸散量進(jìn)一步分為土壤蒸發(fā)量和作物蒸騰量2部分,從而避免了非生產(chǎn)性用水(土壤蒸發(fā))與生產(chǎn)性用水(作物蒸騰)效應(yīng)的混淆。最終的生產(chǎn)量以生物量和收獲指數(shù)來(lái)表示,用以突出水分脅迫對(duì)二者各自的影響。改進(jìn)后的公式為[37]

    式中為最終作物產(chǎn)量,kg/m2;為生物量,kg/m2;HI為收獲指數(shù);WP為生物量水分生產(chǎn)效率,kg/(m2·mm);T為作物蒸騰量,mm。

    2.2 數(shù)據(jù)說(shuō)明

    AquaCrop作物模型所需的輸入數(shù)據(jù)主要包括4個(gè)部分:氣象數(shù)據(jù)、土壤數(shù)據(jù)、作物數(shù)據(jù)和管理數(shù)據(jù)(表1)。在本研究中,氣象數(shù)據(jù)來(lái)源于國(guó)家氣候中心和歐洲中期天氣預(yù)報(bào)中心,主要包含中國(guó)1979—2015年的氣溫、降水、風(fēng)速、太陽(yáng)輻射等;土壤數(shù)據(jù)來(lái)自ISRIC世界土壤信息,主要包含土壤質(zhì)地、土壤有機(jī)質(zhì)含量和土層剖面等信息;作物數(shù)據(jù)包含作物生長(zhǎng)發(fā)育參數(shù)、蒸發(fā)蒸騰參數(shù)、產(chǎn)量形成參數(shù)、脅迫參數(shù)等;管理數(shù)據(jù)則取自于模型默認(rèn)值。此外,各市的市級(jí)玉米單產(chǎn)數(shù)據(jù)來(lái)源于中國(guó)農(nóng)村統(tǒng)計(jì)年鑒以及全國(guó)各省統(tǒng)計(jì)年鑒,用于AquaCrop模型參數(shù)的校正和驗(yàn)證。

    2.3 研究方法

    2.3.1 模型校正與驗(yàn)證

    參數(shù)敏感分析可以有效減少模型校正過(guò)程中的數(shù)據(jù)處理量,大幅度提高工作效率。該文采用擴(kuò)展傅里葉幅度檢驗(yàn)法(extended Fourier amplitude sensitivity test, EFAST)來(lái)進(jìn)行敏感參數(shù)分析[39-40]。EFAST是一種基于方差分解的全局敏感性分析方法,目前已經(jīng)廣泛應(yīng)用于多種作物模型[41-44]。該方法通過(guò)分解模擬結(jié)果的方差獲得各參數(shù)的敏感指數(shù),參數(shù)的敏感性指數(shù)越大表示該參數(shù)對(duì)模型輸出結(jié)果的影響越大。研究中選取2個(gè)敏感性指數(shù)最大的參數(shù)作為待校正參數(shù)。

    在其他非敏感參數(shù)固定的情況下,將待校正的敏感參數(shù)以FAO作物參考手冊(cè)中提供的參考值為初始值,在0.7~1.3范圍內(nèi)以0.02為步長(zhǎng)進(jìn)行隨機(jī)變化,在所有參數(shù)組合下運(yùn)行模型并得到對(duì)應(yīng)的模擬產(chǎn)量,然后與統(tǒng)計(jì)產(chǎn)量進(jìn)行擬合,以二者的歸一化均方根誤差(normalized root mean square error,NRMSE)作為模型校正的評(píng)價(jià)指標(biāo)。通常認(rèn)為,當(dāng)NRMSE達(dá)到10%以下為優(yōu)秀,10%~20%之間為良好,20%~30%之間為尚可,大于30%則說(shuō)明擬合結(jié)果較差。

    完成模型校正工作后,將標(biāo)定的參數(shù)組合輸入模型模擬出各市玉米產(chǎn)量,并與統(tǒng)計(jì)數(shù)據(jù)進(jìn)行擬合,以決定系數(shù)(2)作為評(píng)價(jià)指標(biāo),對(duì)校正結(jié)果進(jìn)行驗(yàn)證。

    2.3.2 脆弱性曲線構(gòu)建

    利用AquaCrop作物模型中的灌溉管理功能進(jìn)行灌溉情景設(shè)置。情景1為完全灌溉情景,在作物模型中選擇灌溉方式為凈灌溉(net irrigation),設(shè)置灌溉水平為100%,即灌溉能夠完全滿足作物生長(zhǎng)需求。情景2為無(wú)灌溉情景,對(duì)應(yīng)灌溉方式中的雨養(yǎng)(rain-fed),即不灌溉,完全依靠自然降水提供作物生長(zhǎng)所需水分的情景。2種情景都處于歷史真實(shí)氣象條件、實(shí)際土壤條件和默認(rèn)田間管理?xiàng)l件下,是否存在水分脅迫是2種灌溉情景的唯一區(qū)別。故認(rèn)為情景1與情景2中模擬產(chǎn)量的差值即為剔除其他環(huán)境脅迫的影響、只由水分脅迫造成的損失,為歷史真實(shí)氣候條件下因干旱導(dǎo)致的產(chǎn)量損失值。

    表1 數(shù)據(jù)說(shuō)明

    脆弱性曲線是一種定量化表達(dá)承災(zāi)體脆弱性的方法,在該文中用來(lái)衡量作物在遭受干旱影響時(shí)可能導(dǎo)致的產(chǎn)量損失。故構(gòu)建脆弱性曲線的關(guān)鍵在于干旱強(qiáng)度指標(biāo)和產(chǎn)量損失指標(biāo)的設(shè)計(jì)和構(gòu)建。

    該文使用作物水分脅迫指數(shù)(crop water stress indicator,CWSI)來(lái)定義干旱強(qiáng)度指標(biāo),式(4)為水分脅迫指數(shù)的計(jì)算方式[45],用作物模型模擬中的實(shí)際蒸散量與潛在蒸散量的比值來(lái)表示。水分脅迫指數(shù)受到氣象條件、土壤性質(zhì)、作物遺傳參數(shù)以及耕作制度等多種因素的共同影響,取值范圍在0~1之間,指數(shù)越大代表作物受到的水分脅迫越強(qiáng)烈。

    式中CWSI為水分脅迫指數(shù),ET和ET0分別為實(shí)際蒸散量和潛在蒸散量,mm。

    考慮到干旱的累積效應(yīng)和干旱持續(xù)時(shí)間的影響,筆者將干旱強(qiáng)度指標(biāo)定義為水分脅迫指數(shù)的日均值:

    式中DHI表示干旱強(qiáng)度指數(shù),CWSI表示第天的水分脅迫指數(shù),表示生育期天數(shù),d。

    產(chǎn)量損失指標(biāo)用產(chǎn)量損失率表示:

    式中YLR表示玉米因水分脅迫導(dǎo)致的產(chǎn)量損失率,1和2分別表示情景1(完全灌溉)和情景2(無(wú)灌溉)條件下的玉米產(chǎn)量,t/hm2。

    從AquaCrop模型的模擬輸出結(jié)果中提取2個(gè)不同灌溉情景下各玉米種植城市歷年的產(chǎn)量,計(jì)算其因干旱導(dǎo)致的產(chǎn)量損失率及對(duì)應(yīng)的干旱強(qiáng)度指標(biāo)。在此基礎(chǔ)上擬合出5個(gè)玉米種植區(qū)的干旱脆弱性曲線。已有研究顯示[15,33-34,46],干旱強(qiáng)度指標(biāo)和產(chǎn)量損失率之間符合logistic函數(shù)關(guān)系。故本研究基于1979—2015年共37 a的長(zhǎng)時(shí)期歷史情境,使用logistic曲線實(shí)現(xiàn)脆弱性曲線擬合。

    3 結(jié)果與分析

    3.1 模型校正與驗(yàn)證

    由參數(shù)敏感性分析得到的2個(gè)敏感參數(shù)分別為作物冠層形成后到衰老之前的作物系數(shù)K和參考收獲指數(shù)HI0(表2)。

    表2 參數(shù)全局敏感性分析排名前10的參數(shù)

    模型校正結(jié)果為241個(gè)城市的標(biāo)準(zhǔn)參數(shù)組。采用線性回歸結(jié)果的斜率、R和nRMSE,通過(guò)比較每個(gè)城市的模擬產(chǎn)量和統(tǒng)計(jì)產(chǎn)量對(duì)經(jīng)過(guò)參數(shù)校正的模型進(jìn)行精度檢驗(yàn),結(jié)果如圖2所示。模擬產(chǎn)量統(tǒng)計(jì)產(chǎn)量的皮爾遜相關(guān)系數(shù)為0.82,擬合直線的斜率為0.81(通過(guò)顯著性檢驗(yàn),在0.01水平上顯著),散點(diǎn)集中在1:1線附近,說(shuō)明利用作物模型模擬出的產(chǎn)量與實(shí)際產(chǎn)量接近;2=0.67,擬合程度較高;nRMSE=17%,介于10%~20%之間,結(jié)果良好??傮w來(lái)說(shuō),模型校正的精度達(dá)到預(yù)期水平。經(jīng)驗(yàn)證后的各城市的率定參數(shù)K為0.315~1.785、HI0為0.144~0.816。

    圖2 模型驗(yàn)證結(jié)果

    3.2 區(qū)域脆弱性曲線

    圖3為中國(guó)5個(gè)玉米種植區(qū)的干旱脆弱性曲線擬合結(jié)果及對(duì)應(yīng)的災(zāi)損函數(shù)表達(dá)式。北方春播玉米種植區(qū)、黃淮海夏播玉米種植區(qū)、西南山地玉米種植區(qū)、南方丘陵玉米種植區(qū)、西北灌溉玉米種植區(qū)的擬和結(jié)果對(duì)應(yīng)的2依次為0.93、0.86、0.47、0.70、0.98,總體來(lái)說(shuō)位于北方的3個(gè)區(qū)域的擬合效果優(yōu)于南方的2個(gè)區(qū)域。由脆弱性曲線擬合結(jié)果圖可以看出,5個(gè)區(qū)域的擬合曲線均近似為“S”形。當(dāng)干旱強(qiáng)度指標(biāo)DHI達(dá)到0.2附近時(shí),產(chǎn)量損失率開始迅速增加;當(dāng)DHI達(dá)到0.6左右時(shí),產(chǎn)量損失率接近最大值。

    觀察圖中散點(diǎn)可知,西北灌溉玉米種植區(qū)各城市單元的干旱強(qiáng)度明顯高于其他地區(qū),西南山地區(qū)和南方丘陵區(qū)的干旱強(qiáng)度則明顯最低。表3對(duì)241個(gè)城市單元的年均干旱強(qiáng)度指數(shù)進(jìn)行了統(tǒng)計(jì)。結(jié)果顯示,西北灌溉玉米種植區(qū)的干旱水平最高,所有城市單元的DHI均在0.2以上,其中81.25%的城市高于0.4。其次是北方春播玉米種植區(qū),近半數(shù)城市單元的DHI達(dá)到0.2以上,20%左右超過(guò)0.4。再次是黃淮海夏播玉米種植區(qū)和南方丘陵玉米種植區(qū)。干旱水平最低的是西南山地玉米種植區(qū),該區(qū)域中所有城市單元的DHI均在0.2以下。

    圖3 不同玉米種植區(qū)區(qū)域干旱脆弱性曲線

    表3 不同干旱強(qiáng)度指標(biāo)下城市個(gè)數(shù)及占比

    4 討 論

    4.1 模型校正與驗(yàn)證

    作物模型校正和驗(yàn)證的精度通常會(huì)受到多種因素的影響,例如校正單元、樣本數(shù)量以及模型輸入數(shù)據(jù)的質(zhì)量等。在現(xiàn)有全國(guó)尺度基于模型模擬的脆弱性曲線研究中,模型校正通常以較大的自然區(qū)域作為基本校正單元[15,33-34,47]。根據(jù)地理學(xué)第一定律[48],如果基本校正單元過(guò)大,則難以保證同一區(qū)域內(nèi)的均質(zhì)化,同一組標(biāo)定參數(shù)不能真實(shí)反映區(qū)域內(nèi)的地理環(huán)境差異,從而導(dǎo)致作物模型的模擬結(jié)果產(chǎn)生較大的誤差。故本研究采用市級(jí)行政單元作為作物模型的基本校正單元,有效提高了同一校正單元內(nèi)地理環(huán)境的相似性,從而提高了模型的模擬運(yùn)算精度。

    除去基本校正單元的選擇,還有一些其他因素可能影響模型參數(shù)標(biāo)定的精度。研究中采用擴(kuò)展傅里葉幅度檢驗(yàn)法選取了2個(gè)敏感參數(shù)進(jìn)行標(biāo)定,但其他非敏感參數(shù)同樣對(duì)模型模擬結(jié)果存在或多或少的影響;同時(shí)由于缺乏實(shí)際的田間管理數(shù)據(jù)(如灌溉、施肥、農(nóng)藥、地表覆蓋等),研究中采用模型默認(rèn)值作為輸入,同樣可能導(dǎo)致模擬結(jié)果的誤差;此外,對(duì)不同玉米品種間作物參數(shù)差異的忽略也可能是模型校正不確定性的來(lái)源之一。

    4.2 干旱指標(biāo)構(gòu)建

    在已有的研究中,基于水分脅迫構(gòu)建干旱指標(biāo)的方法通常采用生育期內(nèi)每日水分脅迫指數(shù)的累加值來(lái)計(jì)算干旱強(qiáng)度指數(shù),常用的2種計(jì)算方法如式(7)[16,32]和式(8)[33,47,49]所示:

    式中DHI代表第年第個(gè)區(qū)域/格網(wǎng)的干旱強(qiáng)度指數(shù),CWSI為第年中第天的水分脅迫指數(shù),max(DHI)和min(DHI)分別代表第個(gè)區(qū)域/格網(wǎng)中所有站點(diǎn)/場(chǎng)景下歷年DHI的最大值和最小值。這2種方法均采用以水分脅迫指數(shù)累加值計(jì)算干旱強(qiáng)度指數(shù)DHI的方法(稱為累加法)。這種方法存在2個(gè)弊端:1)累加法只考慮到作物干旱的累積效應(yīng)而忽略了干旱持續(xù)時(shí)間的影響。當(dāng)作物因干旱提前死亡而導(dǎo)致生育期提前結(jié)束的情況下,模型模擬會(huì)提前結(jié)束,由累加法計(jì)算得到的干旱強(qiáng)度將會(huì)比實(shí)際值偏小,這將可能導(dǎo)致干旱脆弱性曲線擬合過(guò)程中出現(xiàn)DHI值偏小而對(duì)應(yīng)產(chǎn)量損失率接近1的異常點(diǎn),使得脆弱性曲線擬合結(jié)果出現(xiàn)偏差。2)由于累加法計(jì)算得到的DHI值通常大于1,為了使DHI的取值范圍落在0~1之間,上述公式均對(duì)累加結(jié)果進(jìn)行了標(biāo)準(zhǔn)化處理。但是這樣一來(lái),不僅干旱強(qiáng)度指數(shù)失去了其作為水分脅迫指數(shù)的真實(shí)物理意義,同時(shí)不同區(qū)域之間也會(huì)失去可比性。

    本研究中用水分脅迫指數(shù)的日均值(式(6))代替?zhèn)鹘y(tǒng)的累加值,提出了更加合理的干旱強(qiáng)度指數(shù)構(gòu)建方法。當(dāng)一個(gè)情景中不存在作物提前死亡的情況時(shí),使用日均法與累加法構(gòu)建的干旱強(qiáng)度指數(shù)是相同的,而當(dāng)作物因干旱死亡、提前結(jié)束生育期的情景下,使用日均值能夠有效解決干旱指數(shù)計(jì)算小于真實(shí)值的問(wèn)題,從而提高脆弱性曲線擬合的精度;同時(shí)也使得不同情景之間更具可比性,因此日均法的適用范圍更為廣泛。此外,由于水分脅迫指數(shù)的取值范圍本身就介于0~1之間,不需要進(jìn)行額外的標(biāo)準(zhǔn)化處理,故很容易對(duì)不同區(qū)域間干旱強(qiáng)度指數(shù)的大小進(jìn)行統(tǒng)計(jì)和比較。綜上,該文中在傳統(tǒng)累加法的基礎(chǔ)上改進(jìn)后得到的日均法更適用于脆弱性曲線研究中干旱強(qiáng)度指標(biāo)的構(gòu)建。

    4.3 作物在不同生長(zhǎng)階段對(duì)干旱的響應(yīng)

    該文僅研究了玉米在完整生育期中產(chǎn)量受水分脅迫的影響,在構(gòu)建脆弱性曲線的時(shí)候,各生長(zhǎng)階段的水分脅迫大小在干旱致在強(qiáng)度的計(jì)算中是等權(quán)重的。但事實(shí)上,由于作物在不同生長(zhǎng)階段對(duì)水分的需求有所差異,在相同的水分脅迫強(qiáng)度下產(chǎn)生的影響也應(yīng)有所不同。故在未來(lái)的研究中,應(yīng)當(dāng)進(jìn)一步完善對(duì)作物在不同生長(zhǎng)階段對(duì)干旱強(qiáng)度的響應(yīng)規(guī)律,在建立干旱強(qiáng)度指標(biāo)時(shí)對(duì)不同階段的水分脅迫賦予相應(yīng)的權(quán)重,進(jìn)一步增加作物干旱脆弱性曲線的理論和實(shí)用價(jià)值。

    5 結(jié) 論

    該文以中國(guó)5大玉米種植區(qū)域中的241個(gè)玉米種植城市為研究單元,在對(duì)AquaCrop作物模型完成參數(shù)標(biāo)定的基礎(chǔ)上,提出了基于模型模擬的脆弱性曲線構(gòu)建方法。主要結(jié)果如下:

    1)AquaCrop作物模型中,對(duì)玉米產(chǎn)量變化最為敏感的2個(gè)參數(shù)分別是作物冠層形成后到衰老之前的作物系數(shù)和參考收獲指數(shù)。對(duì)這2個(gè)敏感參數(shù)在241個(gè)玉米種植城市中進(jìn)行逐市標(biāo)定,最終得到對(duì)應(yīng)241個(gè)城市的241個(gè)標(biāo)準(zhǔn)參數(shù)組。驗(yàn)證的結(jié)果表明,模型校正的精度較高(2=0.67)。

    2)5個(gè)玉米種植區(qū)的玉米干旱脆弱性曲線擬合結(jié)果均接近“S”形,當(dāng)干旱強(qiáng)度指標(biāo)(drought hazard index, DHI)達(dá)到0.2附近時(shí),產(chǎn)量損失率開始迅速增加;當(dāng)DHI達(dá)到0.6左右時(shí),產(chǎn)量損失率接近最大值。北方春播玉米種植區(qū)、黃淮海夏播玉米種植區(qū)、西南山地玉米種植區(qū)、南方丘陵玉米種植區(qū)、西北灌溉玉米種植區(qū)擬合函數(shù)的決定系數(shù)依次為0.93、0.86、0.47、0.70、0.98,曲線擬合結(jié)果較好,在中國(guó)區(qū)域性玉米干旱脆弱性研究方面具有一定參考價(jià)值,在干旱風(fēng)險(xiǎn)評(píng)估領(lǐng)域具有一定的實(shí)際應(yīng)用價(jià)值。

    [1]Rosenzweig C, Elliott J, Deryng D, et al. Assessing agricultural risks of climate change in the 21st century in a global gridded crop model intercomparison[J]. Proceedings of the National Academy of Sciences, 2014, 111(9): 3268-3273.

    [2]Xie W, Xiong W, Pan J, et al. Decreases in global beer supply due to extreme drought and heat[J]. Nature Plants, 2018, 4(11):964-973.

    [3]IPCC. Climate change 2007: Impact, adaption, and vulnerability[M]. Cambridge: Cambridge University Press, 2007.

    [3]IPCC. Climate change 2014: Impact, adaption, and vulnerability[M]. Cambridge: Cambridge University Press, 2014.

    [4]FAO. Global report on food crises 2017[R/OL]. 2017-03-28 [2019-05-01]. http://www.fao.org

    [5]何斌,武建軍,呂愛(ài)鋒. 農(nóng)業(yè)干旱風(fēng)險(xiǎn)研究進(jìn)展[J]. 地理科學(xué)進(jìn)展,2010,29(5):557-564.

    He Bin, Wu Jianjun, Lü Aifeng. New advances in agricultural drought risk study[J]. Progress in Geography, 2010, 29(5): 557-564. (in Chinese with English abstract)

    [6]Dai A. Increasing drought under global warming in observations and models[J]. Nature Climate Change, 2012, 3(1): 52-58.

    [7]Chen H, Sun J. Changes in drought characteristics over China using the standardized precipitation evapotranspiration index[J]. Journal of Climate, 2015, 28(13): 5430-5447.

    [8]Yao N, Li Y, Lei T, et al. Drought evolution, severity and trends in mainland China over 1961-2013[J]. Science of the Total Environment, 2018, 616/617: 73-89.

    [10]國(guó)家統(tǒng)計(jì)局. 中國(guó)統(tǒng)計(jì)年鑒[M]. 北京:中國(guó)統(tǒng)計(jì)出版社, 2017.

    [11]Lobell D B, Schlenker W, Costa-Roberts J. Climate trends and global crop production since 1980[J]. Science, 2011, 333(6042): 616-620.

    [12]Leng G, Tang Q, Rayburg S. Climate change impacts on meteorological, agricultural and hydrological droughts in China[J]. Global & Planetary Change, 2015, 126(126): 23-34.

    [13]Birkmann J. Measuring vulnerability to promote disaster- resilient societies: Conceptual frameworks and definitions[M]. Measuring Vulnerability to Natural Hazards: Towards Disaster Resilient Societies. Tokyo: United Nations University Press, 2006.

    [14]Wisner B. At risk: Natural hazards, people’s vulnerability and disasters[M]. London: Psychology Press, 2004.

    [15]Wang Z, He F, Fang W, et al. Assessment of physical vulnerability to agricultural drought in China[J]. Natural Hazards, 2013, 67(2): 645-657.

    [16]Guo H, Zhang X, Lian F, et al. Drought risk assessment based on vulnerability surfaces: A case study of maize[J]. Sustainability, 2016, 8(8): 813.

    [17]Dilley M. Natural disaster hotspots: A global risk analysis[J]. Uwe Deichmann, 2005, 20(4): 1-145.

    [18]Wu H, Wilhite D A. An operational agricultural drought risk assessment model for Nebraska, USA[J]. Natural Hazards, 2004, 33(1): 1-21.

    [19]Smith D I. Flood damage estimation: A review of urban stage-damage curves and loss functions[J]. Water SA, 1994, 20(3): 231-238.

    [20]Dutta D, Herath S, Musiakec K. A mathematical model for flood loss estimation[J]. Journal of Hydrology, 2003, 277(1/2): 24-49.

    [21]Singhal A, Kiremidjian A S. Method for probabilistic evaluation of seismic structural damage[J]. Journal of Structural Engineering, 1996, 122(12): 1459-1467.

    [22]Colombi M, Borzi B, Crowley H, et al. Deriving vulnerability curves using Italian earthquake damage data[J]. Bulletin of Earthquake Engineering, 2008, 6(3): 485-504.

    [23]Orsini G. A model for buildings’ vulnerability assessment using the parameterless scale of seismic intensity (PSI)[J]. Earthquake Spectra, 2012, 15(3): 463-483.

    [24]Khanduri A C, Morrow G C. Vulnerability of buildings to windstorms and insurance loss estimation[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2003, 91(4): 455-467.

    [25]Lee K H, Rosowsky D V. Fragility assessment for roof sheathing failure in high wind regions[J]. Engineering Structures, 2005, 27(6): 857-868.

    [26]Bell R, Glade T. Quantitative risk analysis for landslides: Examples from Bildudalur, NW-Iceland[J]. Natural Hazards and Earth System Sciences, 2004, 4(1): 117-131.

    [27]Galli M, Guzzetti F. Landslide vulnerability criteria: A case study from Umbria, central Italy[J]. Environmental Management, 2007, 40(4): 649-664.

    [28]Keylock C J, Barbolini M. Snow avalanche impact pressure - vulnerability relations for use in risk assessment[J]. Canadian Geotechnical Journal, 2001, 38(2): 227-238.

    [29]Cappabianca F, Barbolini M, Natale L. Snow avalanche risk assessment and mapping: A new method based on a combination of statistical analysis, avalanche dynamics simulation and empirically-based vulnerability relations integrated in a GIS platform[J]. Cold Regions Science and Technology, 2008, 54(3): 193-205.

    [30]Hohl R, Schiesser H H, Aller D. Hailfall: The relationship between radar-derived hail kinetic energy and hail damage to buildings[J]. Atmospheric Research, 2002, 63(3): 177-207.

    [31]周瑤,王靜愛(ài). 自然災(zāi)害脆弱性曲線研究進(jìn)展[J]. 地球科學(xué)進(jìn)展,2012,27(4):435-442.

    Zhou Yao, Wang Jingai. A review on development of vulnerability curve of natural disaster[J]. Advance in Earth Sciences, 2012, 27(4): 435-442. (in Chinese with English abstract)

    [32]Yin Y, Zhang X, Lin D, et al. GEPIC-V-R model: A GIS-based tool for regional crop drought risk assessment[J]. Agricultural Water Management, 2014, 144: 107-119.

    [33]Yue Y, Li J, Ye X, et al. An EPIC model-based vulnerability assessment of wheat subject to drought[J]. Natural Hazards, 2015, 78(3): 1629-1652.

    [34]Jia H, Wang J, Cao C, et al. Maize drought disaster risk assessment of China based on EPIC model[J]. International Journal of Digital Earth, 2012, 5(6): 488-515.

    [35]董姝娜,龐澤源,張繼權(quán),等. 基于CERES-Maize模型的吉林西部玉米干旱脆弱性曲線研究[J]. 災(zāi)害學(xué),2014,29(3):115-119.

    Dong Shuna, Pang Zeyuan, Zhang Jiquan, et al. Research on vulnerability curve of drought disaster of maize on CERES-Maize model in western Jilin province[J]. Journal of Catastrophology, 2014, 29(3): 115-119. (in Chinese with English abstract)

    [36]佟屏亞. 中國(guó)玉米種植區(qū)劃[M]. 北京:中國(guó)農(nóng)業(yè)科技出版社,1992.

    [37]朱秀芳,李宜展,潘耀忠,等. AquaCrop作物模型研究和應(yīng)用進(jìn)展[J]. 中國(guó)農(nóng)學(xué)通報(bào),2014,30(8):270-278.

    Zhu Xiufang, Li Yizhan, Pan Yaozhong, et al. A review on the research and application of AquaCrop model[J]. Chinese Agricultural Science Bulletin, 2014, 30(8): 270-278. (in Chinese with English abstract)

    [38]Todorovic M, Albrizio R, Zivotic L, et al. Assessment of AquaCrop, CropSyst, and WOFOST models in the simulation of sunflower growth under different water regimes[J]. Agronomy Journal, 2009, 101(3): 509-521.

    [39]邢會(huì)敏,相詩(shī)堯,徐新剛,等. 基于EFAST方法的AquaCrop作物模型參數(shù)全局敏感性分析[J]. 中國(guó)農(nóng)業(yè)科學(xué),2017,50(1):64-76.

    Xing Huimin, Xiang Shiyao, Xu Xingang, et al. Global sensitivity analysis of AquaC rop crop model parameters based on EFAST method[J]. Scientia Agricultura Sinica, 2017, 50(1): 64-76. (in Chinese with English abstract)

    [40]Saltelli A. Sensitivity analysis: Could better methods be used?[J]. Journal of Geophysical Research Atmospheres, 1999, 104(D3): 3789-3793.

    [41]姜志偉,陳仲新,周清波,等. CERES-Wheat作物模型參數(shù)全局敏感性分析[J]. 農(nóng)業(yè)工程學(xué)報(bào),2011,27(1):236-242.

    Jiang Zhiwei, Chen Zhongxin, Zhou Qingbo, et al. Global sensitivity analysis of CERES-Wheat model parameters[J]. Transaction of the Chinese Society of Agriculture Engineering (Transactions of the CSAE), 2011, 27(1): 236-242. (in Chinese with English abstract)

    [42]DeJonge K C, Li J C A, Ahmadi M, et al. Global sensitivity and uncertainty analysis of a dynamic agroecosystem model under different irrigation treatments[J]. Ecological Modelling, 2012, 231: 113-125.

    [43]Wang J, Li X, Lu L, et al. Parameter sensitivity analysis of crop growth models based on the extended Fourier amplitude sensitivity test method[J]. Environmental Modelling & Software, 2013, 48(5): 171-182.

    [44]Vanuytrecht E, Raes D, Willems P. Global sensitivity analysis of yield output from the water productivity model AquaCrop[J]. Environmental Modelling & Software, 2014, 51(1): 323-332.

    [45]Jackson R D, Idso S B R J, Reginato R J, et al. Canopy temperature as a crop water stress indicator[J]. Water Resources Research, 1981, 17(4): 1133-1138.

    [46]Cui Y, Jiang S, Jin J, et al. Quantitative assessment of soybean drought loss sensitivity at different growth stages based on S-shaped damage curve[J]. Agricultural Water Management, 2019, 213: 821-832.

    [47]Jia H C, Pan D H, Li J, et al. Risk assessment of maize drought disaster in southwest China using the Environmental Policy Integrated Climate model[J]. Journal of Mountain Science, 2016, 13(3): 465-475.

    [48]Tobler W R. A computer movie simulatiing urban growth in the detroit region[J]. Economic Geography, 1970, 46(2): 234-240.

    [49]Wang Z, Fang W, Sun P, et al. Assessment on typical drought risk for wheat production in China based on natural vulnerability[J]. Arid Zone Research, 2010, 27(1): 23-35.

    Vulnerability of drought disaster of maize in China based on AquaCrop model

    Xu Kun1, Zhu Xiufang1,2,3※, Liu Ying1, Hou Chenyao1

    (1100875,; 2.100875,; 3.100875,)

    Drought disaster assessment has become increasingly significant in ensuring national food security and sustainable agricultural development. Vulnerability assessment plays a significant role in disaster research area and vulnerability curve is one of the common quantitative evaluation methods in the field of vulnerability research. In this paper, using the AquaCrop model that has been calibrated city by city, we simulated the response of maize yield to different water stress and then constructed drought vulnerability curves for 5 maize planting regions in China: the north spring maize planting region, the Huang-Huai-Hai summer maize planting region, the southwest mountain maize planting region, the south hilly maize planting region and the northwest irrigated maize planting region. In this research, firstly, 2 of 36 main crop parameters of maize were selected as sensitive parameters based on a global sensitivity analysis method, Extended Fourier Amplitude Sensitivity Test. Then, AquaCrop model was calibrated city by city in 241 maize-growing cities and used to simulate the maize yield under different irrigation scenarios. Finally, we built drought vulnerability curves of 5 main maize plating regions with an improved drought hazard index construction method, which used an average value of daily drought hazard indexes instead of the commom accumulate value, thus we raised comparability of drought hazard index between different maize planting regions and took extreme drought situation into account. The results showed that: 1) The 2 most sensitive parameters to maize yield in the Aquacrop model were the crop coefficient when canopy growth was complete but prior to senescence and the reference harvest index. We finally obtain 241 groups of parameters for the 241 maize planting cities after finishing model calibration and according to the result of validation, the accuracy of the model calibration was satisfactory (2=0.67). 2) All the 5 vulnerability curves followed an “S” shape. And we found that when the drought hazard index reached 0.2, the yield loss rate began to increase rapidly; and when it reached 0.6, the yield loss rate approached the maximum value. The2of the fitted functions in 5 maize planting regions were 0.93, 0.86, 0.47, 0.70, 0.98, respectively. The northwest irrigated maize planting region had the highest2and the southwest mountain maize planting region had the lowest. The drought situation was more serious in the northwest irrigated maize planting region, followed by the north spring maize planting region, the Huang-Huai-Hai summer maize planting region, the south hilly maize planting region and the southwest mountain maize planting region. The research enriched case studies of the AquaCrop model and vulnerability curve construction, quantitatively explored the spatial and temporal differences in drought effects on maize yield in China and enhanced the researches on yield loss prediction. It provides valulble information for the study of drought hazard vulnerability of maize in China and has a certain practical value in the field of drought risk assessment.

    drought; water stress; vulnerability curves; AquaCrop model; maize

    徐 昆,朱秀芳,劉 瑩,侯陳瑤. 采用AquaCrop作物生長(zhǎng)模型研究中國(guó)玉米干旱脆弱性[J]. 農(nóng)業(yè)工程學(xué)報(bào),2020,36(1):154-161.doi:10.11975/j.issn.1002-6819.2020.01.018 http://www.tcsae.org

    Xu Kun, Zhu Xiufang, Liu Ying, Hou Chenyao. Vulnerability of drought disaster of maize in China based on AquaCrop model[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(1): 154-161. (in Chinese with English abstract) doi:10.11975/j.issn.1002-6819.2020.01.018 http://www.tcsae.org

    2019-05-09

    2019-09-10

    國(guó)家重點(diǎn)研發(fā)計(jì)劃(2019YFAO606900);國(guó)家自然科學(xué)基金青年基金項(xiàng)目(41401479)

    徐 昆,博士生,專業(yè)方向?yàn)榈貓D學(xué)與地理信息系統(tǒng)。Email:xukun@mail.bnu.edu.cn

    朱秀芳,副教授,博士,主要研究方向?yàn)檫b感應(yīng)用。Email:zhuxiufang@bnu.edu.cn

    10.11975/j.issn.1002-6819.2020.01.018

    S423; S513

    A

    1002-6819(2020)-01-0154-08

    猜你喜歡
    種植區(qū)脆弱性校正
    C市主要草莓種植區(qū)土壤重金屬鎘、鉛現(xiàn)狀調(diào)查
    不同種植區(qū)隴東苜蓿營(yíng)養(yǎng)價(jià)值的比較研究
    劉光第《南旋記》校正
    草莓種植區(qū)土壤中典型有機(jī)氮化合物的分布及來(lái)源
    一類具有校正隔離率隨機(jī)SIQS模型的絕滅性與分布
    機(jī)內(nèi)校正
    煤礦電網(wǎng)脆弱性評(píng)估
    電子制作(2017年10期)2017-04-18 07:23:09
    殺毒軟件中指令虛擬機(jī)的脆弱性分析
    基于攻擊圖的工控系統(tǒng)脆弱性量化方法
    河北昌黎縣葡萄種植區(qū)農(nóng)業(yè)地球化學(xué)特征
    18禁裸乳无遮挡动漫免费视频| 国产亚洲欧美精品永久| 性色av乱码一区二区三区2| 女人久久www免费人成看片| 亚洲成人国产一区在线观看| 欧美老熟妇乱子伦牲交| 国产精品九九99| 日韩熟女老妇一区二区性免费视频| videosex国产| 欧美精品高潮呻吟av久久| 建设人人有责人人尽责人人享有的| 热re99久久精品国产66热6| 国产精品亚洲一级av第二区| 欧美不卡视频在线免费观看 | 久久天堂一区二区三区四区| av欧美777| 999精品在线视频| 中文字幕精品免费在线观看视频| 国产免费男女视频| 欧美日韩成人在线一区二区| 桃红色精品国产亚洲av| a级毛片在线看网站| 国产伦人伦偷精品视频| 国产精品久久久久久人妻精品电影| 大型av网站在线播放| 成人国语在线视频| 国产精品乱码一区二三区的特点 | 成人永久免费在线观看视频| 老熟女久久久| 日韩欧美一区视频在线观看| 久久 成人 亚洲| 精品熟女少妇八av免费久了| 18禁黄网站禁片午夜丰满| 亚洲国产中文字幕在线视频| 欧美性长视频在线观看| 国产高清激情床上av| 亚洲国产精品一区二区三区在线| svipshipincom国产片| 欧美日韩瑟瑟在线播放| 亚洲国产看品久久| 日韩大码丰满熟妇| 人人妻人人添人人爽欧美一区卜| 天天躁夜夜躁狠狠躁躁| 久久中文字幕一级| 欧美成人免费av一区二区三区 | 天天躁日日躁夜夜躁夜夜| ponron亚洲| 国产成人系列免费观看| 丁香欧美五月| 亚洲专区字幕在线| 精品少妇久久久久久888优播| 黄色成人免费大全| 免费黄频网站在线观看国产| 日韩熟女老妇一区二区性免费视频| av超薄肉色丝袜交足视频| 亚洲五月天丁香| 极品人妻少妇av视频| av国产精品久久久久影院| 午夜两性在线视频| 亚洲国产看品久久| 亚洲九九香蕉| 99香蕉大伊视频| 大片电影免费在线观看免费| 50天的宝宝边吃奶边哭怎么回事| 国产免费男女视频| 丰满饥渴人妻一区二区三| avwww免费| 亚洲精品国产精品久久久不卡| 亚洲aⅴ乱码一区二区在线播放 | 岛国在线观看网站| 777久久人妻少妇嫩草av网站| 亚洲三区欧美一区| 啦啦啦 在线观看视频| 国产欧美日韩综合在线一区二区| 亚洲专区字幕在线| 在线视频色国产色| av超薄肉色丝袜交足视频| 亚洲精品国产色婷婷电影| 亚洲专区国产一区二区| 久久中文看片网| 国产精品久久视频播放| 高潮久久久久久久久久久不卡| 婷婷成人精品国产| 视频区欧美日本亚洲| 国产麻豆69| 午夜91福利影院| 久久国产乱子伦精品免费另类| 一夜夜www| 国产精品免费大片| 美女视频免费永久观看网站| 大陆偷拍与自拍| 99久久精品国产亚洲精品| 咕卡用的链子| 国产高清激情床上av| 国产精品久久视频播放| 99精品在免费线老司机午夜| 国产欧美日韩一区二区三区在线| 精品乱码久久久久久99久播| 久久久国产精品麻豆| 一级毛片女人18水好多| 欧美人与性动交α欧美软件| 纯流量卡能插随身wifi吗| 精品国产美女av久久久久小说| 亚洲黑人精品在线| 精品无人区乱码1区二区| 日韩欧美在线二视频 | 1024香蕉在线观看| 中文字幕人妻丝袜制服| 国产国语露脸激情在线看| 在线播放国产精品三级| 午夜日韩欧美国产| 手机成人av网站| 国产男靠女视频免费网站| 90打野战视频偷拍视频| 精品一区二区三区四区五区乱码| 亚洲精品av麻豆狂野| 天堂俺去俺来也www色官网| 丁香欧美五月| 亚洲男人天堂网一区| 色94色欧美一区二区| 麻豆av在线久日| 精品久久久久久电影网| 国产日韩一区二区三区精品不卡| 成人18禁高潮啪啪吃奶动态图| 两人在一起打扑克的视频| 国产精品一区二区免费欧美| 欧美精品啪啪一区二区三区| 在线观看www视频免费| 午夜福利一区二区在线看| 中文字幕另类日韩欧美亚洲嫩草| 高清在线国产一区| 亚洲精品在线观看二区| 免费看十八禁软件| 午夜日韩欧美国产| 热99久久久久精品小说推荐| 91国产中文字幕| 搡老乐熟女国产| 丁香六月欧美| 91九色精品人成在线观看| 免费在线观看视频国产中文字幕亚洲| 精品福利观看| 久久精品亚洲熟妇少妇任你| 色综合婷婷激情| 色综合婷婷激情| 老司机靠b影院| 免费高清在线观看日韩| av天堂久久9| 老熟妇仑乱视频hdxx| 俄罗斯特黄特色一大片| 9色porny在线观看| 一级a爱片免费观看的视频| 亚洲五月天丁香| 亚洲精品久久午夜乱码| 国产精品久久久av美女十八| 国产欧美日韩一区二区精品| 99国产精品一区二区三区| 国产不卡一卡二| 黄片大片在线免费观看| 搡老熟女国产l中国老女人| 80岁老熟妇乱子伦牲交| 久久这里只有精品19| 曰老女人黄片| 久久久国产精品麻豆| 婷婷丁香在线五月| 91字幕亚洲| 51午夜福利影视在线观看| 国产日韩一区二区三区精品不卡| 涩涩av久久男人的天堂| 后天国语完整版免费观看| 国产精品电影一区二区三区 | 美女午夜性视频免费| 777久久人妻少妇嫩草av网站| 不卡一级毛片| 80岁老熟妇乱子伦牲交| 母亲3免费完整高清在线观看| 亚洲精品乱久久久久久| 久久精品国产综合久久久| 老熟女久久久| 亚洲七黄色美女视频| 中出人妻视频一区二区| 欧美乱色亚洲激情| 国产日韩欧美亚洲二区| 男女免费视频国产| 男女高潮啪啪啪动态图| 欧美乱码精品一区二区三区| 久久久久久久国产电影| 久久久水蜜桃国产精品网| 制服诱惑二区| 淫妇啪啪啪对白视频| av免费在线观看网站| 精品视频人人做人人爽| 免费观看人在逋| 在线观看舔阴道视频| 国产一区有黄有色的免费视频| 久久精品亚洲av国产电影网| 手机成人av网站| 亚洲精品美女久久av网站| 狂野欧美激情性xxxx| 日本黄色视频三级网站网址 | 午夜免费鲁丝| 国产真人三级小视频在线观看| 午夜精品在线福利| 久久久久国产精品人妻aⅴ院 | 香蕉久久夜色| 成年人免费黄色播放视频| 久久精品国产99精品国产亚洲性色 | 精品久久久久久,| 一区二区三区精品91| 韩国av一区二区三区四区| 国产蜜桃级精品一区二区三区 | 午夜福利欧美成人| 亚洲欧美激情综合另类| 妹子高潮喷水视频| 日韩三级视频一区二区三区| 69精品国产乱码久久久| 多毛熟女@视频| 老司机在亚洲福利影院| av在线播放免费不卡| 男女下面插进去视频免费观看| 大片电影免费在线观看免费| 亚洲av成人不卡在线观看播放网| 一区二区三区精品91| 中文字幕精品免费在线观看视频| 久久久国产精品麻豆| 国产片内射在线| 亚洲伊人色综图| 丝袜人妻中文字幕| 好看av亚洲va欧美ⅴa在| 亚洲精品一二三| 热99re8久久精品国产| 91九色精品人成在线观看| 午夜福利视频在线观看免费| 俄罗斯特黄特色一大片| 欧美丝袜亚洲另类 | 校园春色视频在线观看| 亚洲专区国产一区二区| 日韩视频一区二区在线观看| 99精品在免费线老司机午夜| 欧美黄色片欧美黄色片| 黑人猛操日本美女一级片| 女人高潮潮喷娇喘18禁视频| 欧美日韩视频精品一区| 久久国产精品影院| 亚洲九九香蕉| 男男h啪啪无遮挡| 午夜福利一区二区在线看| 操美女的视频在线观看| 男人舔女人的私密视频| 欧美人与性动交α欧美精品济南到| 757午夜福利合集在线观看| 19禁男女啪啪无遮挡网站| 色尼玛亚洲综合影院| 涩涩av久久男人的天堂| 美女 人体艺术 gogo| 一二三四在线观看免费中文在| 久久中文字幕人妻熟女| 在线观看免费高清a一片| 国产欧美日韩一区二区三| 十分钟在线观看高清视频www| 亚洲精品国产区一区二| 黄片大片在线免费观看| 国产精品国产高清国产av | 97人妻天天添夜夜摸| 露出奶头的视频| 久久国产精品人妻蜜桃| 少妇被粗大的猛进出69影院| 老司机福利观看| 老熟女久久久| 18禁黄网站禁片午夜丰满| 久久香蕉精品热| 午夜亚洲福利在线播放| 免费一级毛片在线播放高清视频 | 国产不卡一卡二| 日本vs欧美在线观看视频| 亚洲精品国产区一区二| 国产精品成人在线| 亚洲欧美日韩另类电影网站| 成在线人永久免费视频| 看片在线看免费视频| 一进一出好大好爽视频| 一级片'在线观看视频| 亚洲av成人不卡在线观看播放网| 国产精品美女特级片免费视频播放器 | 多毛熟女@视频| 无遮挡黄片免费观看| 日韩有码中文字幕| 飞空精品影院首页| 久久午夜亚洲精品久久| 成人特级黄色片久久久久久久| cao死你这个sao货| 欧美性长视频在线观看| 男女之事视频高清在线观看| 天天躁日日躁夜夜躁夜夜| 岛国毛片在线播放| 欧美最黄视频在线播放免费 | 国产真人三级小视频在线观看| 高清毛片免费观看视频网站 | 欧美丝袜亚洲另类 | 999精品在线视频| 精品久久久久久久毛片微露脸| 婷婷丁香在线五月| 日韩中文字幕欧美一区二区| 十八禁网站免费在线| 久久ye,这里只有精品| 国产不卡一卡二| 国产精品亚洲av一区麻豆| 日韩制服丝袜自拍偷拍| 18禁黄网站禁片午夜丰满| 不卡一级毛片| 后天国语完整版免费观看| 女同久久另类99精品国产91| 最新的欧美精品一区二区| 色老头精品视频在线观看| 香蕉丝袜av| 中文字幕高清在线视频| 国产成人精品在线电影| 亚洲午夜理论影院| 黄色视频,在线免费观看| 亚洲 欧美一区二区三区| a级片在线免费高清观看视频| 美女扒开内裤让男人捅视频| 国产不卡一卡二| 可以免费在线观看a视频的电影网站| 丝瓜视频免费看黄片| 一区在线观看完整版| 中文亚洲av片在线观看爽 | 国产成人欧美| 国产欧美日韩一区二区精品| 国产高清视频在线播放一区| 久9热在线精品视频| 热re99久久精品国产66热6| 最新的欧美精品一区二区| 国产精品1区2区在线观看. | 美女高潮喷水抽搐中文字幕| 黑人欧美特级aaaaaa片| 免费在线观看完整版高清| www.熟女人妻精品国产| 国产一卡二卡三卡精品| 91大片在线观看| 午夜福利视频在线观看免费| 最新在线观看一区二区三区| 久久久国产欧美日韩av| 十分钟在线观看高清视频www| 亚洲精品在线观看二区| 国产97色在线日韩免费| 久久久国产欧美日韩av| 91麻豆av在线| avwww免费| 欧美日韩亚洲国产一区二区在线观看 | 国产精品.久久久| 精品少妇久久久久久888优播| 欧美黑人欧美精品刺激| 国产99久久九九免费精品| www.999成人在线观看| 91大片在线观看| 18禁国产床啪视频网站| 精品久久久久久久久久免费视频 | 国产精品一区二区在线观看99| 国精品久久久久久国模美| 久久久水蜜桃国产精品网| 脱女人内裤的视频| 女人久久www免费人成看片| 久久人妻av系列| 国产精品免费一区二区三区在线 | avwww免费| 熟女少妇亚洲综合色aaa.| 巨乳人妻的诱惑在线观看| 可以免费在线观看a视频的电影网站| 色在线成人网| 大片电影免费在线观看免费| 国产精品久久久人人做人人爽| 午夜久久久在线观看| 一区二区日韩欧美中文字幕| 制服人妻中文乱码| 高清视频免费观看一区二区| 一级毛片精品| 亚洲人成77777在线视频| 国产成人精品久久二区二区免费| 欧美日韩亚洲综合一区二区三区_| 国产亚洲精品一区二区www | 亚洲欧美激情在线| 中文欧美无线码| 一区二区三区激情视频| 麻豆国产av国片精品| 国产精品综合久久久久久久免费 | 人人妻人人澡人人爽人人夜夜| 一区二区三区国产精品乱码| 国产精品自产拍在线观看55亚洲 | 中亚洲国语对白在线视频| 少妇 在线观看| 在线观看午夜福利视频| 欧美日韩黄片免| 久9热在线精品视频| 纯流量卡能插随身wifi吗| 无人区码免费观看不卡| 久久国产乱子伦精品免费另类| 可以免费在线观看a视频的电影网站| 一级a爱片免费观看的视频| 大片电影免费在线观看免费| 少妇被粗大的猛进出69影院| 国产成人av激情在线播放| 女人高潮潮喷娇喘18禁视频| 中文欧美无线码| 身体一侧抽搐| 国产极品粉嫩免费观看在线| 国产在线一区二区三区精| 色婷婷av一区二区三区视频| 精品久久久久久,| 91麻豆av在线| 精品一品国产午夜福利视频| 欧美丝袜亚洲另类 | 丝瓜视频免费看黄片| 一边摸一边抽搐一进一小说 | 人成视频在线观看免费观看| 在线观看免费高清a一片| 一二三四社区在线视频社区8| 亚洲精品国产色婷婷电影| 两个人免费观看高清视频| 欧美成人免费av一区二区三区 | 亚洲人成电影观看| 亚洲色图av天堂| 午夜亚洲福利在线播放| 国产成人精品久久二区二区免费| 亚洲熟妇熟女久久| 国产aⅴ精品一区二区三区波| 天天影视国产精品| 久久久久久亚洲精品国产蜜桃av| 黄色 视频免费看| 最近最新免费中文字幕在线| 午夜日韩欧美国产| 欧美激情 高清一区二区三区| 自拍欧美九色日韩亚洲蝌蚪91| 大型黄色视频在线免费观看| 亚洲五月婷婷丁香| 免费日韩欧美在线观看| 成人亚洲精品一区在线观看| 精品久久久久久久久久免费视频 | 在线观看一区二区三区激情| 国内久久婷婷六月综合欲色啪| 午夜福利,免费看| 美女扒开内裤让男人捅视频| 人人妻,人人澡人人爽秒播| 乱人伦中国视频| 十八禁人妻一区二区| 国产亚洲av高清不卡| 老鸭窝网址在线观看| 岛国在线观看网站| 亚洲九九香蕉| 少妇猛男粗大的猛烈进出视频| 高清视频免费观看一区二区| 伦理电影免费视频| 午夜福利欧美成人| 国产区一区二久久| www.自偷自拍.com| 欧美日韩一级在线毛片| 身体一侧抽搐| 中文欧美无线码| 天堂俺去俺来也www色官网| 亚洲一码二码三码区别大吗| 久久国产亚洲av麻豆专区| 久久精品国产清高在天天线| 欧美一级毛片孕妇| 国产男靠女视频免费网站| av福利片在线| 日韩中文字幕欧美一区二区| 亚洲五月天丁香| 97人妻天天添夜夜摸| 亚洲aⅴ乱码一区二区在线播放 | 午夜福利一区二区在线看| 午夜免费鲁丝| 日韩欧美免费精品| 久久香蕉精品热| 操出白浆在线播放| 天天影视国产精品| 黑人欧美特级aaaaaa片| 成在线人永久免费视频| 一本大道久久a久久精品| 精品国产一区二区三区久久久樱花| 丝袜美足系列| 久久香蕉国产精品| 亚洲在线自拍视频| 19禁男女啪啪无遮挡网站| 久久国产精品人妻蜜桃| 国产欧美日韩精品亚洲av| 亚洲欧洲精品一区二区精品久久久| 国产不卡一卡二| 一级片'在线观看视频| 免费黄频网站在线观看国产| av片东京热男人的天堂| 99久久人妻综合| 精品一区二区三区四区五区乱码| 老司机靠b影院| 99riav亚洲国产免费| 国产精品免费视频内射| 亚洲五月天丁香| 欧美精品亚洲一区二区| 久久久久久久精品吃奶| 免费日韩欧美在线观看| а√天堂www在线а√下载 | 99精国产麻豆久久婷婷| av电影中文网址| 国内毛片毛片毛片毛片毛片| netflix在线观看网站| 国产精品98久久久久久宅男小说| 亚洲av日韩精品久久久久久密| 正在播放国产对白刺激| 欧美午夜高清在线| 国产高清国产精品国产三级| 亚洲国产看品久久| 成人三级做爰电影| 99re在线观看精品视频| 亚洲国产中文字幕在线视频| 精品国产一区二区三区四区第35| 熟女少妇亚洲综合色aaa.| av欧美777| 99久久综合精品五月天人人| 男女午夜视频在线观看| 中出人妻视频一区二区| 午夜福利欧美成人| 777米奇影视久久| 亚洲国产精品sss在线观看 | 成在线人永久免费视频| 久久精品国产综合久久久| 天天添夜夜摸| 97人妻天天添夜夜摸| 国产精品美女特级片免费视频播放器 | 咕卡用的链子| 十八禁人妻一区二区| 亚洲少妇的诱惑av| 国产aⅴ精品一区二区三区波| 黑丝袜美女国产一区| 变态另类成人亚洲欧美熟女 | 一区在线观看完整版| 国产又爽黄色视频| 国精品久久久久久国模美| 大型av网站在线播放| 女人被躁到高潮嗷嗷叫费观| 亚洲精品国产一区二区精华液| 捣出白浆h1v1| 下体分泌物呈黄色| 热99re8久久精品国产| 80岁老熟妇乱子伦牲交| 王馨瑶露胸无遮挡在线观看| 最新的欧美精品一区二区| 免费看a级黄色片| 国产人伦9x9x在线观看| 国内毛片毛片毛片毛片毛片| 精品一区二区三区视频在线观看免费 | 成人黄色视频免费在线看| 天天躁狠狠躁夜夜躁狠狠躁| 久久人人爽av亚洲精品天堂| 狠狠狠狠99中文字幕| 少妇粗大呻吟视频| 色综合欧美亚洲国产小说| 19禁男女啪啪无遮挡网站| 欧美日韩精品网址| 极品人妻少妇av视频| 欧美激情 高清一区二区三区| 国产在视频线精品| 久久午夜亚洲精品久久| 国产不卡一卡二| 一区二区三区国产精品乱码| 久久久久国内视频| 在线视频色国产色| 亚洲 欧美一区二区三区| 久久亚洲真实| av中文乱码字幕在线| 国产在线精品亚洲第一网站| 男女免费视频国产| 大型黄色视频在线免费观看| 正在播放国产对白刺激| 免费观看人在逋| 12—13女人毛片做爰片一| 午夜成年电影在线免费观看| 黑丝袜美女国产一区| 妹子高潮喷水视频| 熟女少妇亚洲综合色aaa.| 大香蕉久久网| 欧美日韩精品网址| 99久久精品国产亚洲精品| 99精品欧美一区二区三区四区| 成年人午夜在线观看视频| 午夜老司机福利片| 18禁国产床啪视频网站| 国产精品乱码一区二三区的特点 | 两性夫妻黄色片| 国产不卡一卡二| 五月开心婷婷网| 日本黄色视频三级网站网址 | 欧美日韩亚洲国产一区二区在线观看 | 亚洲全国av大片| 99精品久久久久人妻精品| 亚洲成人免费电影在线观看| 免费观看a级毛片全部| av线在线观看网站| x7x7x7水蜜桃| 巨乳人妻的诱惑在线观看| 国产成人免费无遮挡视频| 色精品久久人妻99蜜桃| 精品久久久久久,| 日日摸夜夜添夜夜添小说| 国产黄色免费在线视频| 日韩免费av在线播放| 十分钟在线观看高清视频www| 人人妻,人人澡人人爽秒播| 国产欧美日韩精品亚洲av| 精品高清国产在线一区| 99热只有精品国产| 美女福利国产在线| 亚洲av欧美aⅴ国产| 国产高清视频在线播放一区| 久久久国产成人精品二区 | 宅男免费午夜| 王馨瑶露胸无遮挡在线观看| 亚洲国产欧美一区二区综合|