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

    基于Copula函數(shù)的陜西蘋果晚霜凍特征分析

    2022-03-10 02:26:26茹曉雅李美榮王景紅蘇寶峰何建強(qiáng)
    關(guān)鍵詞:霜凍歷時(shí)花期

    姜 元,茹曉雅,羅 琦,李美榮,王 釗,王景紅,馮 浩,張 東,蘇寶峰,于 強(qiáng),何建強(qiáng)

    ·農(nóng)業(yè)生物環(huán)境與能源工程·

    基于Copula函數(shù)的陜西蘋果晚霜凍特征分析

    姜 元1,2,茹曉雅1,2,羅 琦1,2,李美榮3,王 釗3,王景紅3,馮 浩2,4,張 東5,蘇寶峰6,于 強(qiáng)4,何建強(qiáng)1,2※

    (1. 西北農(nóng)林科技大學(xué)旱區(qū)農(nóng)業(yè)水土工程教育部重點(diǎn)實(shí)驗(yàn)室,楊凌 712100;2. 西北農(nóng)林科技大學(xué)中國旱區(qū)節(jié)水農(nóng)業(yè)研究院,楊凌 712100;3. 陜西省農(nóng)業(yè)遙感與經(jīng)濟(jì)作物氣象服務(wù)中心,西安 710015;4. 中國科學(xué)院水利部水土保持研究所黃土高原土壤侵蝕與旱地農(nóng)業(yè)國家重點(diǎn)實(shí)驗(yàn)室,楊凌 712100;5. 西北農(nóng)林科技大學(xué)園藝學(xué)院,楊凌 712100;6. 西北農(nóng)林科技大學(xué)農(nóng)業(yè)農(nóng)村部農(nóng)業(yè)物聯(lián)網(wǎng)重點(diǎn)實(shí)驗(yàn)室,楊凌 712100)

    為探討利用Copula函數(shù)對(duì)蘋果晚霜凍進(jìn)行特征分析的適用性,該研究首先基于陜西省蘋果主產(chǎn)區(qū)7個(gè)氣象站點(diǎn)1971-2018年的逐日最低氣溫(min)數(shù)據(jù)集,提取出晚霜凍事件的歷時(shí)和強(qiáng)度兩個(gè)特征變量。然后,基于6種不同的Copula函數(shù)構(gòu)建晚霜凍特征變量的聯(lián)合分布,并進(jìn)行擬合優(yōu)度評(píng)價(jià)。最后,利用優(yōu)選的Copula函數(shù)分析晚霜凍發(fā)生的概率及重現(xiàn)期。結(jié)果表明:陜西蘋果產(chǎn)區(qū)各站點(diǎn)1971-2018年受晚霜凍的影響在空間分布上由東南向西北方向加重,各站點(diǎn)晚霜凍的歷時(shí)和強(qiáng)度之間均具有顯著的正相關(guān)關(guān)系。當(dāng)晚霜凍強(qiáng)度和歷時(shí)增大時(shí),其聯(lián)合累積概率也相應(yīng)增大,且增大趨勢(shì)變緩。各站點(diǎn)聯(lián)合重現(xiàn)期代表的“或”事件比同現(xiàn)重現(xiàn)期所代表的“且”事件更容易發(fā)生。當(dāng)單變量重現(xiàn)期取值較小時(shí),可將聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期視為單變量重現(xiàn)期的兩種極端情況,對(duì)其實(shí)際范圍進(jìn)行估計(jì)??傮w而言,陜西蘋果產(chǎn)區(qū)各站點(diǎn)發(fā)生長歷時(shí)且高強(qiáng)度晚霜凍事件的概率較小,但是位于延安果區(qū)的站點(diǎn)相較于其他果區(qū)站點(diǎn)更容易發(fā)生高強(qiáng)度或長歷時(shí)的晚霜凍事件,以及高強(qiáng)度長歷時(shí)同時(shí)發(fā)生的晚霜凍事件,需要重點(diǎn)加以關(guān)注。該研究可為陜西蘋果產(chǎn)區(qū)應(yīng)對(duì)晚霜凍災(zāi)害提供理論依據(jù)。

    氣象;災(zāi)害;晚霜凍;Copula函數(shù);頻率分析;重現(xiàn)期;蘋果

    0 引 言

    蘋果在中國的栽培面積和產(chǎn)量均位居世界首位[1-2]。蘋果的廣泛栽培使得蘋果生產(chǎn)更易受到各種氣象災(zāi)害的影響,晚霜凍是中國北方蘋果產(chǎn)區(qū)面臨的首要極端氣象災(zāi)害[3-5],會(huì)對(duì)蘋果產(chǎn)量和質(zhì)量有嚴(yán)重影響[4,6],給當(dāng)?shù)靥O果產(chǎn)業(yè)造成重大經(jīng)濟(jì)損失[4]。目前,全球氣候變暖的形勢(shì)日趨嚴(yán)峻[7],冬季增溫幅度也更為明顯[8],暖冬導(dǎo)致果樹春季萌芽及開花提前[9-10],使其遭受晚霜凍的風(fēng)險(xiǎn)日益增加。陜西地區(qū)的大陸性季風(fēng)氣候特征十分明顯,全年約有40%的強(qiáng)降溫天氣發(fā)生在3-4月份。然而,該時(shí)期正處于蘋果開花期,頻繁的冷空氣過程導(dǎo)致晚霜凍災(zāi)害風(fēng)險(xiǎn)較高[11-13]。

    目前,關(guān)于蘋果晚霜凍的研究多集中在建立蘋果晚霜凍指數(shù)等級(jí)并探究晚霜凍事件時(shí)空變化特征和成災(zāi)因子。例如,李健等[14]基于氣象站最低氣溫構(gòu)建了蘋果果樹花期凍害預(yù)警指標(biāo),并分析了主果區(qū)花期凍害的時(shí)空變化特征。屈振江等[13]構(gòu)建了凍害風(fēng)險(xiǎn)災(zāi)損率對(duì)蘋果花期凍害風(fēng)險(xiǎn)進(jìn)行了評(píng)估。柏秦鳳等[15]基于花期凍害風(fēng)險(xiǎn)指數(shù)對(duì)蘋果花期凍害進(jìn)行降尺度風(fēng)險(xiǎn)區(qū)劃。王景紅等[12]利用人工氣候箱模擬低溫并結(jié)合歷史災(zāi)情調(diào)查對(duì)蘋果花期霜凍指標(biāo)進(jìn)行了修訂優(yōu)化。張琪等[3]基于物候模型探究了蘋果霜凍發(fā)生特征。邱美娟等[1]結(jié)合晚霜凍氣象指標(biāo),對(duì)蘋果花期晚霜凍氣候風(fēng)險(xiǎn)進(jìn)行了評(píng)估。王明昌等[11]采用線性傾向法和風(fēng)險(xiǎn)指數(shù)法分析了陜西省禮泉和旬邑兩縣蘋果花期凍害風(fēng)險(xiǎn)情況。

    上述研究大多基于晚霜凍指標(biāo)進(jìn)行蘋果晚霜凍災(zāi)害風(fēng)險(xiǎn)的分析,但是評(píng)估霜凍等氣象災(zāi)害事件風(fēng)險(xiǎn)的較佳方法是通過詳細(xì)了解以持續(xù)時(shí)間和低溫強(qiáng)度為特征的霜凍事件的發(fā)生頻率[16]。李美榮等[17]建立了蘋果花期嚴(yán)重凍害最低氣溫的概率分布模型及重現(xiàn)期預(yù)測(cè),但是只使用了單變量頻率分析方法且僅考慮了影響霜凍的一個(gè)變量,無法捕捉到極端氣象災(zāi)害問題的復(fù)雜性。

    另一種評(píng)估氣象災(zāi)害事件風(fēng)險(xiǎn)更為有效的的方法是建立多元聯(lián)合分布[18]。然而,許多傳統(tǒng)多元聯(lián)合分布要求各變量的邊緣分布屬于同一類型[19]。Copula函數(shù)作為一種構(gòu)造靈活的模型不限制各變量邊緣分布的選擇,能夠通過邊緣分布和相關(guān)性結(jié)構(gòu)兩部分來構(gòu)建多維聯(lián)合分布,形式靈活多樣,被廣泛應(yīng)用于降雨頻率[20]、干旱頻率[21]和洪水頻率[22]等水文和災(zāi)害事件分析。Chatrabgoun等[16]首次采用Copula函數(shù)對(duì)葡萄園霜凍進(jìn)行了概率評(píng)估,通過對(duì)霜凍持續(xù)時(shí)間和嚴(yán)重程度這2個(gè)特征變量建立聯(lián)合分布,探究了位于伊朗馬來爾地區(qū)葡萄園霜凍現(xiàn)象的風(fēng)險(xiǎn)和影響,表明Copula函數(shù)可以作為構(gòu)建霜凍多維特征變量聯(lián)合分布函數(shù)的有效工具。然而,將Copula函數(shù)應(yīng)用于蘋果晚霜凍特征分析的研究目前還鮮有報(bào)道。

    本研究以陜西省蘋果主產(chǎn)區(qū)7個(gè)地面氣象觀測(cè)站1971-2018年的氣象數(shù)據(jù)集為基礎(chǔ)數(shù)據(jù),采用Copula函數(shù)對(duì)陜西蘋果產(chǎn)區(qū)的晚霜凍災(zāi)害進(jìn)行概率評(píng)估,研究目標(biāo)包括:1)逐個(gè)氣象站提取晚霜凍事件發(fā)生的連續(xù)天數(shù)和日最低溫度最小值的絕對(duì)值分別作為表征晚霜凍歷時(shí)和強(qiáng)度的2個(gè)特征變量;2)對(duì)2個(gè)特征變量的單變量邊緣分布進(jìn)行擬合,并分析變量間相關(guān)性;3)基于6種Copula函數(shù),建立2個(gè)特征變量的多元聯(lián)合分布模型并進(jìn)行擬合優(yōu)度檢驗(yàn);4)分析計(jì)算聯(lián)合累積概率和重現(xiàn)期,從而對(duì)陜西蘋果產(chǎn)區(qū)各站點(diǎn)的晚霜凍災(zāi)害風(fēng)險(xiǎn)進(jìn)行評(píng)估。本研究的結(jié)果有望在站點(diǎn)尺度上為陜西蘋果晚霜凍災(zāi)害防御和蘋果產(chǎn)業(yè)生態(tài)管理提供科學(xué)依據(jù)。

    1 材料和方法

    1.1 研究區(qū)概況與數(shù)據(jù)來源

    陜西蘋果主產(chǎn)區(qū)位于黃土高原地區(qū),是符合蘋果生長氣象指標(biāo)的優(yōu)質(zhì)產(chǎn)區(qū)之一。根據(jù)氣候特點(diǎn)及物候差異,將陜西蘋果主產(chǎn)區(qū)劃分為延安果區(qū)、渭北東部果區(qū)、渭北西部果區(qū)和關(guān)中西部果區(qū)[23](圖1)。本研究使用的1971-2018年陜西省氣象站點(diǎn)的逐日最低氣溫(min)數(shù)據(jù)來源于中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)(http://cdc.cma.gov.cn/)。使用該氣象數(shù)據(jù)前需對(duì)各個(gè)站點(diǎn)進(jìn)行檢查,刪除數(shù)據(jù)序列長度過短或缺失的站點(diǎn),并對(duì)數(shù)據(jù)異常值進(jìn)行處理。最終選擇研究區(qū)域內(nèi)的延安市寶塔區(qū)、洛川、長武、旬邑、銅川、白水和禮泉等7個(gè)氣象站作為研究站點(diǎn)(圖1)。由于陜西蘋果花期物候觀測(cè)數(shù)據(jù)不夠完整,所以本文使用的蘋果花期物候數(shù)據(jù)主要來自王潤紅等[24]基于物候模型對(duì)陜西省蘋果花期進(jìn)行模擬所得的數(shù)據(jù)。

    1.2 研究方法

    1.2.1 蘋果晚霜凍事件識(shí)別

    劉映寧等[26]根據(jù)蘋果晚霜凍氣象指標(biāo),結(jié)合災(zāi)情調(diào)查資料,發(fā)現(xiàn)蘋果花期內(nèi)不同程度的低溫會(huì)對(duì)蘋果生長發(fā)育產(chǎn)生不同的影響:日最低溫度低于?4 ℃時(shí),出現(xiàn)嚴(yán)重晚霜凍,中心花受凍率高達(dá)80%,減產(chǎn)30%以上;日最低溫度在?2至?4 ℃時(shí),出現(xiàn)中度晚霜凍,中心花受凍率達(dá)60%至80%,對(duì)產(chǎn)量、品質(zhì)、商品率產(chǎn)生嚴(yán)重影響;日最低溫度在0至?2 ℃時(shí),出現(xiàn)輕度晚霜凍,中心花受凍率達(dá)30%至60%,對(duì)產(chǎn)量、品質(zhì)、商品率產(chǎn)生明顯影響。此外,相關(guān)文獻(xiàn)也給出了蘋果花期凍害的臨界溫度。李美榮等[26]將日最低氣溫min為0、?2、?4 ℃作為蘋果在花期內(nèi)不同時(shí)期受凍的臨界溫度;劉映寧等[26]提出把日最低氣溫分別低于0、?2、?4 ℃作為陜西蘋果花期凍害農(nóng)業(yè)保險(xiǎn)的3個(gè)等級(jí);王景紅等[27]則采用0、?2、?4 ℃作為陜西蘋果不同等級(jí)花期凍害的臨界指標(biāo)。因此,本研究將日最低氣溫min為0 作為發(fā)生蘋果晚霜凍事件的臨界溫度,每個(gè)晚霜凍事件中日最低溫度最小值的絕對(duì)值被作為該期間的晚霜凍強(qiáng)度。蘋果晚霜凍不僅與低溫強(qiáng)度有關(guān),也與低溫的持續(xù)時(shí)間相關(guān),強(qiáng)度越大或持續(xù)時(shí)間越長受凍率越高[12],因此將晚霜凍發(fā)生的連續(xù)天數(shù)定義為晚霜凍歷時(shí)。

    圖1 陜西蘋果產(chǎn)區(qū)亞區(qū)劃分及氣象站點(diǎn)分布

    1.2.2 構(gòu)建晚霜凍歷時(shí)和晚霜凍強(qiáng)度的邊緣分布函數(shù)

    在基于Copula函數(shù)構(gòu)建晚霜凍歷時(shí)和強(qiáng)度聯(lián)合分布之前,首先要確定晚霜凍歷時(shí)和強(qiáng)度這2個(gè)特征變量的邊緣分布函數(shù),同時(shí)要考慮它們之間的相依性。本文選用正態(tài)分布、指數(shù)分布、伽馬分布、威爾布分布、柯西分布、邏輯斯諦分布、對(duì)數(shù)正態(tài)分布7種概率分布函數(shù)構(gòu)建這2個(gè)特征變量邊緣函數(shù)。使用R語言的“fitdistrplus”包,利用極大似然法估計(jì)相關(guān)參數(shù),對(duì)晚霜凍歷時(shí)和強(qiáng)度這2個(gè)特征變量的邊緣分布函數(shù)進(jìn)行擬合,并通過Kolmogorov-Smirnov(K-S)檢驗(yàn)選取最優(yōu)的單變量邊緣分布函數(shù)。采取3種常用的相關(guān)系數(shù)來進(jìn)行2個(gè)晚霜凍特征變量間的相關(guān)性度量,包括皮爾遜積矩相關(guān)系數(shù)(Pearson product-moment correlation coefficient)、斯皮爾曼秩相關(guān)系數(shù)(Spearman’s rank correlation coefficient)、肯德爾秩相關(guān)系數(shù)(Kendall rank correlation coefficient)。

    以陜西省禮泉站為例,因?yàn)槭褂弥鹑兆畹蜌鉁貋碜R(shí)別晚霜凍事件,所以晚霜凍歷時(shí)數(shù)據(jù)均為整數(shù)天,導(dǎo)致晚霜凍歷時(shí)存在較多的相同值,因此晚霜凍歷時(shí)的經(jīng)驗(yàn)分布呈條帶狀(圖2)。在對(duì)晚霜凍歷時(shí)邊緣分布函數(shù)進(jìn)行擬合過程中,這種條帶狀的離散分布會(huì)導(dǎo)致樣本個(gè)數(shù)明顯少于實(shí)際值,進(jìn)而影響擬合效果。參考韓會(huì)明等[28-29]引入干旱歷時(shí)離散變量連續(xù)化處理方法,本研究對(duì)整數(shù)的晚霜凍歷時(shí)加上[?0.5, 0.5]的均勻隨機(jī)分布變量,從而使離散化的數(shù)據(jù)序列連續(xù)化作為處理后的經(jīng)驗(yàn)分布。此外,Michele等[30]證明,對(duì)樣本添加均勻分布的擾動(dòng)不會(huì)改變?cè)嫉慕y(tǒng)計(jì)量信息。

    圖2 禮泉站晚霜凍歷時(shí)邊緣分布擬合情況

    1.2.3 構(gòu)建晚霜凍歷時(shí)和晚霜凍強(qiáng)度的Copula聯(lián)合分布函數(shù)

    Copula函數(shù)可用于構(gòu)造邊緣分布不同的多個(gè)變量的聯(lián)合分布函數(shù)。根據(jù)Sklar定理[31],設(shè)=(<)和=(<)分別為特征變量晚霜凍歷時(shí)和強(qiáng)度的邊緣分布函數(shù),其聯(lián)合分布函數(shù)為(,)(式(1))。

    (,)=(<,<)=[(),()](1)

    式中為特征變量取值,為Copula函數(shù),為邊緣分布函數(shù)。

    確定單個(gè)變量的邊緣分布函數(shù)后,要優(yōu)選一種Copula函數(shù)連接單變量分布函數(shù)。本文選用了6種常用的二維Copula聯(lián)合分布函數(shù)(表1),包括Archimedean Copula函數(shù)中的Gumbel Copula,F(xiàn)rank Copula,Clayton Copula以及Joe Copula,和橢圓Copula函數(shù)中最常見的Normal Copula和T Copula[16,20-21]。其中,參數(shù)估計(jì)均選用極大似然法進(jìn)行估計(jì),擬合優(yōu)度的檢驗(yàn)選用赤池信息準(zhǔn)則(Akaike Information Criterion, AIC;式(2))和貝葉斯信息準(zhǔn)則(Bayesian Information Criterion, BIC;式(3)),以赤池信息量AIC值和貝葉斯信息量BIC值最小作為擬合最優(yōu)的檢驗(yàn)標(biāo)準(zhǔn)。

    AIC=?2ln()+2(2)

    BIC=?2ln()+ln()·(3)

    式中為參數(shù)的數(shù)量;為似然函數(shù);為數(shù)據(jù)數(shù)量。

    表1 本研究中所采用的6種二維Copula函數(shù)

    注:為特征變量取值,為邊緣分布函數(shù),為關(guān)聯(lián)參數(shù),為生成函數(shù)自變量。

    Notes:are characteristic variable values;,are the marginal distribution functions;denotes to the association parameters;is the generator function argument.

    1.2.4 基于晚霜凍歷時(shí)和強(qiáng)度的晚霜凍事件重現(xiàn)期確定

    重現(xiàn)期作為水文分析中評(píng)估風(fēng)險(xiǎn)的一個(gè)通用標(biāo)準(zhǔn),被廣泛用于確定重大災(zāi)害事件發(fā)生的概率,極具現(xiàn)實(shí)意義和應(yīng)用價(jià)值。晚霜凍是一種對(duì)蘋果生長極具威脅性的氣象災(zāi)害,采用重現(xiàn)期對(duì)晚霜凍事件進(jìn)行評(píng)估,有助于了解晚霜凍事件的發(fā)生規(guī)模與發(fā)生頻率之間的關(guān)系,對(duì)晚霜凍事件的發(fā)生進(jìn)行預(yù)測(cè),有助于陜西蘋果產(chǎn)區(qū)更好地防御晚霜凍災(zāi)害。一般情況下,重現(xiàn)期的分析僅限于單變量,所謂單變量重現(xiàn)期(式(4))是指變量超過某一特定值出現(xiàn)一次的間隔時(shí)間,然而由于極端氣象災(zāi)害問題的復(fù)雜性,導(dǎo)致單變量重現(xiàn)期分析結(jié)果具有一定的不準(zhǔn)確性[17]。晚霜凍事件是晚霜凍歷時(shí)和強(qiáng)度兩個(gè)因素相互作用的結(jié)果,所以在晚霜凍事件特征分析中應(yīng)該重點(diǎn)考慮重現(xiàn)期的聯(lián)合屬性和條件屬性。對(duì)于多特征變量的聯(lián)合分布而言,聯(lián)合重現(xiàn)期(Joint return period;式(5))代表多個(gè)特征變量中某一個(gè)特征變量大于給定閾值出現(xiàn)一次的“或”事件的間隔時(shí)間,而同現(xiàn)重現(xiàn)期(Co-occurrence return period;式(6))代表多個(gè)特征變量中每一個(gè)特征變量均大于給定閾值出現(xiàn)一次的“且”事件的間隔時(shí)間。

    傳統(tǒng)基于單變量重現(xiàn)期的計(jì)算公式為

    式中F()為變量的邊緣分布函數(shù);單變量同理。

    二維聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期的計(jì)算公式分別為[32]:

    2 結(jié)果與分析

    2.1 陜西省蘋果主產(chǎn)區(qū)晚霜凍事件特征分析

    本研究利用陜西省蘋果主產(chǎn)區(qū)各個(gè)氣象站點(diǎn)數(shù)據(jù),提取出1971-2018年間該區(qū)域蘋果花期發(fā)生晚霜凍事件的晚霜凍歷時(shí)和晚霜凍強(qiáng)度(表2),以及每次晚霜凍發(fā)生的起始和終止日期。該區(qū)域不同果區(qū)的站點(diǎn)發(fā)生晚霜凍的頻率存在明顯差異,其中位于延安和渭北西部果區(qū)的4個(gè)站點(diǎn)共發(fā)生晚霜凍頻率高,共計(jì)266次,且晚霜凍的歷時(shí)和強(qiáng)度值也相應(yīng)較大;渭北東部和關(guān)中西部的站點(diǎn)發(fā)生晚霜凍共89次,頻次相對(duì)較少,且歷時(shí)和強(qiáng)度值也較小。在陜西蘋果產(chǎn)區(qū)各氣象站中,延安站不僅嚴(yán)重晚霜凍和總晚霜凍發(fā)生次數(shù)最多,而且晚霜凍強(qiáng)度平均值也最大;旬邑站的晚霜凍歷時(shí)和強(qiáng)度的最大值均為最大。

    基于晚霜凍歷時(shí)的起止時(shí)間進(jìn)一步分析晚霜凍歷時(shí)過程(圖3),可見總體上陜西省蘋果主產(chǎn)區(qū)晚霜凍事件歷時(shí)多分布于3月31日至4月16日。關(guān)中西部果區(qū)的站點(diǎn)(圖3a)明顯發(fā)生晚霜凍的次數(shù)少、歷時(shí)短、強(qiáng)度低;渭北西部果區(qū)(圖3d和圖3e)和延安果區(qū)(圖3f和圖3g)的站點(diǎn)發(fā)生的晚霜凍次數(shù)多、歷時(shí)長、強(qiáng)度高。陜西蘋果主產(chǎn)區(qū)各站點(diǎn)在1980年晚霜凍最多發(fā)生達(dá)17次;在2013年發(fā)生晚霜凍共12次,為近10年晚霜凍發(fā)生次數(shù)最多的年份。

    表2 陜西省蘋果主產(chǎn)區(qū)1971-2018年晚霜凍事件統(tǒng)計(jì)

    圖3 陜西蘋果產(chǎn)區(qū)7個(gè)氣象站的晚霜凍歷時(shí)分布

    2.2 晚霜凍歷時(shí)和強(qiáng)度邊緣分布函數(shù)的擬合

    本研究采用正態(tài)分布、指數(shù)分布、伽馬分布、威布爾分布、柯西分布、邏輯斯諦(Logistic)分布、對(duì)數(shù)正態(tài)分布等7種概率分布,對(duì)陜西蘋果產(chǎn)區(qū)晚霜凍歷時(shí)和強(qiáng)度的邊緣分布函數(shù)進(jìn)行擬合,得到晚霜凍歷時(shí)和強(qiáng)度邊緣分布的K-S檢驗(yàn)結(jié)果以及變量之間的相關(guān)系數(shù)(表3)。選取擬合度最高的概率分布類型(顯著性水平=0.05),最終確定各站點(diǎn)晚霜凍歷時(shí)邊緣分布函數(shù)為對(duì)數(shù)正態(tài)分布。對(duì)于晚霜凍強(qiáng)度這一特征變量的邊緣分布擬合,各站點(diǎn)情況不一,其中禮泉站、白水站和旬邑站為Logistic分布,銅川站、洛川站和延安站為指數(shù)分布,長武站為伽馬分布。相關(guān)性檢驗(yàn)結(jié)果表明各站點(diǎn)Pearson、Spearman和Kendall這3個(gè)相關(guān)系數(shù)均大于0.4且都通過了顯著性水平=0.05的顯著性檢驗(yàn),說明各站點(diǎn)的晚霜凍歷時(shí)和強(qiáng)度具有一定的相關(guān)性,可基于Copula函數(shù)建立聯(lián)合分布函數(shù)進(jìn)行進(jìn)一步分析。

    表3 單變量邊緣分布的K-S檢驗(yàn)和變量間的相關(guān)系數(shù)

    注:*表示通過了顯著性水平= 0.05的顯著性檢驗(yàn),**表示通過了顯著性水平= 0.01的顯著性檢驗(yàn)。

    Notes: * is statistically significance level of α= 0.05; ** is statistically extreme significant level of= 0.01.

    2.3 基于Copula函數(shù)晚霜凍事件聯(lián)合分布函數(shù)的建立

    本研究基于擬合的晚霜凍歷時(shí)和強(qiáng)度的邊緣分布函數(shù),建立了晚霜凍歷時(shí)和強(qiáng)度之間的6種Copula分布函數(shù),利用極大似然法對(duì)其中未知參數(shù)進(jìn)行了估計(jì)。此外,還計(jì)算了相應(yīng)的AIC、BIC值(式(2)和式(3)),對(duì)Copula分布函數(shù)擬合優(yōu)度進(jìn)行評(píng)價(jià)(圖4)。根據(jù)AIC和BIC準(zhǔn)則判斷可知,Normal Copula函數(shù)對(duì)禮泉站和旬邑站晚霜凍事件的擬合效果最好;Clayton Copula函數(shù)對(duì)白水站和銅川站晚霜凍事件的擬合效果最好;對(duì)洛川站晚霜凍事件擬合最優(yōu)的是Frank Copula函數(shù);對(duì)延安站和長武站晚霜凍事件擬合優(yōu)度最高的是Joe Copula函數(shù)。

    圖4 陜西蘋果產(chǎn)區(qū)晚霜凍歷時(shí)和強(qiáng)度的6種Copula分布函數(shù)擬合優(yōu)度檢驗(yàn)

    2.4 晚霜凍歷時(shí)和強(qiáng)度的聯(lián)合概率

    進(jìn)一步進(jìn)行本研究中各站點(diǎn)晚霜凍特征變量二維聯(lián)合概率分布分析(圖5),可知當(dāng)晚霜凍強(qiáng)度和歷時(shí)值增大時(shí),其聯(lián)合累積概率也相應(yīng)增大??傮w而言,各站點(diǎn)所在區(qū)域發(fā)生短歷時(shí)低強(qiáng)度、短歷時(shí)高強(qiáng)度、長歷時(shí)低強(qiáng)度的晚霜凍事件的概率較大,而同時(shí)滿足長歷時(shí)和高強(qiáng)度的晚霜凍事件的發(fā)生概率較小。具體地,在禮泉站(圖5a),晚霜凍強(qiáng)度低于3 ℃、晚霜凍歷時(shí)在0.5至2 d時(shí),其聯(lián)合累積概率隨特征變量的增大而迅速增大;在白水站(圖5b),晚霜凍強(qiáng)度高于3 ℃、晚霜凍歷時(shí)高于3 d時(shí),其聯(lián)合累積概率增大趨勢(shì)明顯變緩;在銅川站(圖5c),晚霜凍強(qiáng)度不超過4℃、晚霜凍歷時(shí)在0.5至3 d時(shí),其聯(lián)合累積概率隨特征變量的增大而迅速增大;在旬邑站(圖5d),晚霜凍強(qiáng)度不超過2 ℃、晚霜凍歷時(shí)不超過4 d時(shí),其聯(lián)合累積概率隨特征變量的增大而迅速增大;在長武站(圖5e),晚霜凍強(qiáng)度超過1 ℃、晚霜凍歷時(shí)超過2 d時(shí),其聯(lián)合累積概率增大趨勢(shì)明顯變緩;在洛川站(圖5f),晚霜凍強(qiáng)度不超過2 ℃、晚霜凍歷時(shí)不超過3 d時(shí),其聯(lián)合累積概率隨特征變量的增大而迅速增大;在延安站(圖5g),晚霜凍強(qiáng)度超過1 ℃、晚霜凍歷時(shí)超過2 d時(shí),其聯(lián)合累積概率增大趨勢(shì)明顯變緩。

    圖5 陜西蘋果產(chǎn)區(qū)7個(gè)氣象站晚霜凍歷時(shí)和強(qiáng)度聯(lián)合概率分布

    2.5 晚霜凍歷時(shí)和強(qiáng)度的重現(xiàn)期

    以等值線圖的形式對(duì)晚霜凍事件的聯(lián)合重現(xiàn)期進(jìn)行描繪(圖6)。根據(jù)晚霜凍兩個(gè)特征變量值垂直相交的交點(diǎn)所在等值線可查出聯(lián)合重現(xiàn)期,可發(fā)現(xiàn)隨著晚霜凍特征變量不斷增大,兩變量的聯(lián)合重現(xiàn)期也隨之增大。根據(jù)等線圖的疏密情況可明顯看到,不同站點(diǎn)間晚霜凍事件的聯(lián)合重現(xiàn)期存在一定的差異,各站點(diǎn)中禮泉站最不容易發(fā)生高強(qiáng)度或長歷時(shí)的晚霜凍事件(圖6a),而白水站和銅川站較不易發(fā)生高強(qiáng)度或長歷時(shí)的晚霜凍事件(圖6b和6c)。對(duì)于其他站點(diǎn),若霜凍強(qiáng)度達(dá)到8 ℃或者晚霜凍歷時(shí)達(dá)到8 d時(shí),在旬邑站其聯(lián)合重現(xiàn)期約為150 a(圖6d),在長武站其聯(lián)合重現(xiàn)期約為280 a(圖6e),在洛川站其聯(lián)合重現(xiàn)期約為150 a(圖6f),在延安站其聯(lián)合重現(xiàn)期約為50 a(圖6g)。這說明相較于旬邑站、長武站和洛川站,延安站更容易受到高強(qiáng)度或長歷時(shí)晚霜凍事件的威脅,也進(jìn)一步表明聯(lián)合重現(xiàn)期等值線圖可以表征不同區(qū)域間的晚霜凍事件發(fā)生概率,能夠服務(wù)于減災(zāi)預(yù)防工作。

    通過單變量重現(xiàn)期(式(4))以及構(gòu)建的單變量邊緣分布的逆函數(shù),分別計(jì)算當(dāng)單變量重現(xiàn)期為3、5、10、20、50、100 a時(shí),各站點(diǎn)晚霜凍歷時(shí)和強(qiáng)度值,再計(jì)算聯(lián)合重現(xiàn)期(式(5))和同現(xiàn)重現(xiàn)期(式(6))對(duì)應(yīng)值(表4)。結(jié)果發(fā)現(xiàn)聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期均隨著晚霜凍特征變量的增大而增大,這與等值線圖(圖6)所示的聯(lián)合重現(xiàn)期變化規(guī)律相一致。在單變量取值同等增幅條件下,同現(xiàn)重現(xiàn)期的增幅要明顯高于聯(lián)合重現(xiàn)期,這說明該區(qū)域各站點(diǎn)聯(lián)合重現(xiàn)期代表的“或”事件比同現(xiàn)重現(xiàn)期所代表的“且”事件更容易發(fā)生。然而,不同站點(diǎn)的重現(xiàn)期變化存在一定的差異,例如,當(dāng)單變量重現(xiàn)期增幅相同時(shí),延安站對(duì)應(yīng)的同現(xiàn)重現(xiàn)期增幅明顯大于旬邑站,這說明延安站更易受到長歷時(shí)且高強(qiáng)度晚霜凍事件的威脅。各站的單變量重現(xiàn)期、聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期具有相同的大小順序,單變量重現(xiàn)期總介于聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期之間,因此可以通過計(jì)算聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期來估計(jì)單變量重現(xiàn)期的范圍,即可得到不同晚霜凍特征變量所代表的晚霜凍事件發(fā)生的頻率,但是僅當(dāng)單變量重現(xiàn)期較小時(shí)估計(jì)的范圍才比較準(zhǔn)確。

    圖6 陜西蘋果產(chǎn)區(qū)7個(gè)氣象站晚霜凍歷時(shí)和強(qiáng)度聯(lián)合重現(xiàn)期等值線

    表4 陜西蘋果產(chǎn)區(qū)7個(gè)氣象站晚霜凍歷時(shí)和強(qiáng)度的單變量重現(xiàn)期和聯(lián)合分布重現(xiàn)期

    續(xù)表

    3 討 論

    在對(duì)極端氣象災(zāi)害事件進(jìn)行特征分析時(shí),首先要明確氣象災(zāi)害指標(biāo)以及事件相關(guān)的特征變量,要能夠?qū)O端氣象災(zāi)害事件進(jìn)行綜合性評(píng)價(jià)。以往對(duì)晚霜凍災(zāi)害的研究重點(diǎn)是探究低溫強(qiáng)度的影響,然而張曉煜等[33]提出低溫和低溫持續(xù)時(shí)間是對(duì)霜凍危害程度起決定作用的關(guān)鍵因子;王景紅等[12]表示蘋果晚霜凍不僅與低溫強(qiáng)度有關(guān),也與低溫持續(xù)時(shí)間密切相關(guān),溫度越低或持續(xù)時(shí)間越長,蘋果受凍率越高;袁佰順等[6]發(fā)現(xiàn)0 以下溫度持續(xù)時(shí)間和低溫強(qiáng)度都需作為判斷霜凍受災(zāi)和受災(zāi)程度的因素。因此,本研究選取晚霜凍歷時(shí)和晚霜凍強(qiáng)度作為特征變量,基于Copula函數(shù)建立特征變量的多元聯(lián)合分布,克服了單變量頻率分析方法不能捕捉到極端氣象災(zāi)害復(fù)雜性的不足,為陜西蘋果產(chǎn)區(qū)晚霜凍特征研究提供了新的思路和方法。

    本研究通過對(duì)陜西蘋果產(chǎn)區(qū)晚霜凍事件的特征分析可知,在空間上陜西蘋果產(chǎn)區(qū)受晚霜凍的影響總體上呈現(xiàn)東南向西北加重的趨勢(shì),這與前人對(duì)陜西蘋果晚霜凍的探究結(jié)果基本一致[13-14,17,25-26,34]。該空間分布趨勢(shì)是由于延安果區(qū)和渭北果區(qū)所處的經(jīng)緯度、高海拔及地勢(shì)起伏大所導(dǎo)致的。李美榮等[17]對(duì)旬邑站進(jìn)行單變量重現(xiàn)期計(jì)算,得到10 a重現(xiàn)期對(duì)應(yīng)的晚霜凍強(qiáng)度為5.2 ℃,而本研究為3.81 ℃,這可能是由于研究所選的時(shí)間尺度不同、識(shí)別晚霜凍所用的氣象指標(biāo)不同,以及擬合分布類型不同所導(dǎo)致的。此外,本研究采用Copula函數(shù)通過對(duì)晚霜凍歷時(shí)和強(qiáng)度這2個(gè)特征變量建立多元聯(lián)合分布評(píng)估陜西蘋果晚霜凍風(fēng)險(xiǎn),相比于李美榮等[17]僅考慮低溫強(qiáng)度而使用單變量頻率的分析方法,更能反映晚霜凍事件受晚霜凍歷時(shí)和強(qiáng)度兩個(gè)特征變量共同影響的復(fù)雜性,從而為精準(zhǔn)評(píng)估陜西蘋果晚霜凍風(fēng)險(xiǎn)提供更為全面和科學(xué)的信息。

    準(zhǔn)確的蘋果物候數(shù)據(jù)是評(píng)估蘋果花期凍害風(fēng)險(xiǎn)的重要依據(jù)。然而,目前國內(nèi)關(guān)于蘋果物候觀測(cè)數(shù)據(jù)的積累十分有限,存在時(shí)間序列短、缺測(cè)嚴(yán)重等問題[11],并且現(xiàn)有的驅(qū)動(dòng)物候模型模擬的氣象數(shù)據(jù)也是有限的,從中國氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)所獲得的氣象數(shù)據(jù)基本只能到縣域尺度,缺乏空間密集的氣象觀測(cè)數(shù)據(jù),難以對(duì)蘋果晚霜凍事件進(jìn)行區(qū)域精準(zhǔn)識(shí)別,所以未來研究要加強(qiáng)對(duì)陜西蘋果產(chǎn)區(qū)蘋果花期的多站點(diǎn)和時(shí)序觀測(cè),采用質(zhì)量較高的氣象數(shù)據(jù)產(chǎn)品來驅(qū)動(dòng)物候模型進(jìn)行模擬。蘋果屬于多年生作物,不同地區(qū)、不同樹齡、不同品種的蘋果對(duì)氣候響應(yīng)有所不同,然而目前晚霜凍指標(biāo)沒有進(jìn)行針對(duì)性的細(xì)化,使得目前對(duì)晚霜凍事件的識(shí)別工作具有一定的局限性。本文僅對(duì)晚霜凍歷時(shí)和強(qiáng)度這2個(gè)特征變量進(jìn)行聯(lián)合,后續(xù)研究可以增加對(duì)晚霜凍特征量的選取,但是多元聯(lián)合分布函數(shù)也會(huì)隨著變量維數(shù)增多變得愈發(fā)復(fù)雜。所以,蘋果花期模擬、晚霜凍指標(biāo)細(xì)化、晚霜凍事件精準(zhǔn)識(shí)別、基于Copula函數(shù)的更高維特征變量分析等都是今后進(jìn)行蘋果晚霜凍災(zāi)害研究需要關(guān)注的重點(diǎn)方向。

    4 結(jié) 論

    本研究根據(jù)日最低氣溫閾值提取1971–2018年間陜西蘋果產(chǎn)區(qū)各氣象站晚霜凍事件的歷時(shí)和強(qiáng)度,基于Copula函數(shù)建立了晚霜凍歷時(shí)和強(qiáng)度的聯(lián)合分布,并分析了研究區(qū)內(nèi)晚霜凍事件的重現(xiàn)期概況,得到以下主要結(jié)論:

    1)位于延安和渭北西部果區(qū)的站點(diǎn)共發(fā)生晚霜凍266次,是主要的晚霜凍災(zāi)害頻發(fā)區(qū),并且歷時(shí)長、強(qiáng)度高;渭北東部和關(guān)中西部的站點(diǎn)發(fā)生晚霜凍共89次,頻次相對(duì)較少,而且歷時(shí)短、強(qiáng)度低。陜西蘋果產(chǎn)區(qū)各站點(diǎn)受晚霜凍的影響在空間分布上呈由東南向西北加重的趨勢(shì)。

    2)陜西蘋果產(chǎn)區(qū)各站點(diǎn)晚霜凍歷時(shí)的最優(yōu)邊緣分布均為對(duì)數(shù)正態(tài)分布,晚霜凍強(qiáng)度的最優(yōu)邊緣分布隨站點(diǎn)而各異。各站點(diǎn)晚霜凍的歷時(shí)和強(qiáng)度之間的相關(guān)系數(shù)均大于0.4,具有顯著的正相關(guān)關(guān)系。

    3)本研究采用各站點(diǎn)擬合優(yōu)度最高的Copula 函數(shù)建立陜西蘋果產(chǎn)區(qū)各站點(diǎn)晚霜凍歷時(shí)和強(qiáng)度的聯(lián)合分布函數(shù),分析晚霜凍特征變量的聯(lián)合累積概率可知,同時(shí)滿足長歷時(shí)和高強(qiáng)度的晚霜凍事件的發(fā)生概率較小。

    4)不同站點(diǎn)間晚霜凍事件的重現(xiàn)期存在一定的差異,若霜凍強(qiáng)度達(dá)到8℃或者晚霜凍歷時(shí)達(dá)到8 d時(shí),延安站聯(lián)合重現(xiàn)期約為50 a,相較于其他站點(diǎn)百年以上的聯(lián)合重現(xiàn)期,延安站更容易受到高強(qiáng)度或長歷時(shí),以及高強(qiáng)度且長歷時(shí)晚霜凍事件的威脅。陜西蘋果產(chǎn)區(qū)各站點(diǎn)聯(lián)合重現(xiàn)期代表的“或”事件比同現(xiàn)重現(xiàn)期所代表的“且”事件更容易發(fā)生。當(dāng)單變量重現(xiàn)期的值較小時(shí),可將聯(lián)合重現(xiàn)期和同現(xiàn)重現(xiàn)期可看作單變量重現(xiàn)期的兩種極端情況,對(duì)單變量重現(xiàn)期的實(shí)際范圍進(jìn)行估計(jì)。

    總之,本文基于Copula函數(shù)建立的聯(lián)合分布模型能夠同時(shí)描述晚霜凍事件的歷時(shí)和強(qiáng)度之間的關(guān)系,可以利用重現(xiàn)期評(píng)估陜西蘋果產(chǎn)區(qū)站點(diǎn)尺度未來晚霜凍發(fā)生的風(fēng)險(xiǎn),從而為研究該區(qū)域各站點(diǎn)晚霜凍事件的發(fā)生規(guī)律,提高該區(qū)域各站點(diǎn)晚霜凍的預(yù)報(bào)能力提供新的科學(xué)方法和依據(jù)。

    [1] 邱美娟,劉布春,劉園,等. 中國北方主產(chǎn)地蘋果始花期模擬及晚霜凍風(fēng)險(xiǎn)評(píng)估[J]. 農(nóng)業(yè)工程學(xué)報(bào),2020,36(21):154-163.

    Qiu Meijuan, Liu Buchun, Liu Yuan, et al. Simulation of first flowering date for apple and risk assessment of late frost in main producing areas of northern China[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(21): 154-163. (in Chinese with English abstract)

    [2] 萬煒,師紀(jì)博,劉忠,等. 棲霞市蘋果園氮磷養(yǎng)分平衡及環(huán)境風(fēng)險(xiǎn)評(píng)價(jià)[J]. 農(nóng)業(yè)工程學(xué)報(bào),2020,36(4):211-219.

    Wan Wei, Shi Jibo, Liu Zhong, et al. Nitrogen and phosphorus nutrient balance and environmental risk assessment of apple orchard in Qixia city[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2020, 36(4): 211-219. (in Chinese with English abstract)

    [3] 張琪,李榮平,張曉月,等. 綏中蘋果霜凍災(zāi)害特征及其對(duì)氣候變暖的適應(yīng)分析[J]. 北方園藝,2016(13):30-33.

    Zhang Qi, Li Rongping, Zhang Xiaoyue, et al. Analysis of apple frost disaster characteristics and adapt to the climate warming in Suizhong[J]. Northern Horticulture, 2016(13): 30-33. (in Chinese with English abstract)

    [4] 馬延慶,劉長民,朱海利,等. 陜西咸陽渭北旱塬地區(qū)優(yōu)質(zhì)蘋果基地生態(tài)氣候特征分析[J]. 干旱地區(qū)農(nóng)業(yè)研究,2008(1):146-153.

    Ma Yanqing, Liu Changmin, Zhu Haili, et al. The analysis of ecosystem characteristic of high quality apple producing base on Weibei arid plain region in Xianyang of Shaanxi Province[J]. Agricultural Research in the Arid Areas, 2008(1): 146-153. (in Chinese with English abstract)

    [5] 李化龍,趙西社. 陜西黃土高原果業(yè)氣候生態(tài)條件研究及應(yīng)用[M]. 北京:氣象出版社,2010.

    [6] 袁佰順,許彥平,姚曉紅,等. 基于春季晚霜凍害的天水蘋果氣候風(fēng)險(xiǎn)區(qū)劃[J]. 干旱區(qū)資源與環(huán)境,2015,29(2):185-189.

    Yuan Baishun, Xu Yanping, Yao Xiaohong, et al. The compartment of later frost climate risk for apple in Tianshui[J]. Journal of Arid Land Resources and Environment, 2015, 29(2): 185-189. (in Chinese with English abstract)

    [7] Organization W M. WMO Provisional Statement on the State of the Global Climate in 2019[EB/OL]. http://www.cma.gov.cn/en2014/climate/featutes/202003/P020200312816904145935.pdf, 2022-6-23.

    [8] Cook B I, Wolkovich E M, Parmesan C, Divergent responses to spring and winter warming drive community level flowering trends[J]. 2012, 109(23): 9000-9005.

    [9] Ge Q, Wang H, Rutishauser T, et al. Phenological response to climate change in China: A meta-analysis[J]. Global Change Biology, 2015, 21(1): 265-274.

    [10] 郭佳,張寶林,高聚林,等. 氣候變化對(duì)中國農(nóng)業(yè)氣候資源及農(nóng)業(yè)生產(chǎn)影響的研究進(jìn)展[J]. 北方農(nóng)業(yè)學(xué)報(bào),2019,47(1):105-113.

    Guo Jia, Zhang Baolin, Gao Julin, et al. Advances on the impacts of climate change on agro-climatic resources and agricultural production in China[J]. Journal of Northern Agriculture, 2019, 47(1): 105-113. (in Chinese with English abstract)

    [11] 王明昌,劉布春,劉園,等. 陜西蘋果主產(chǎn)縣花期凍害風(fēng)險(xiǎn)評(píng)估[J]. 中國農(nóng)業(yè)氣象,2020,41(6):381-392.

    Wang Mingchang, Liu Buchun, Liu Yuan, et al. Assessment on the freezing injury risk during apple flowering in Liquan and Xunyi[J]. Chinese Journal of Agrometeorolog, 2020, 41(6): 381-392. (in Chinese with English abstract)

    [12] 王景紅,劉璐,高峰,等. 陜西富士系蘋果花期霜凍災(zāi)害氣象指標(biāo)的修訂[J]. 中國農(nóng)業(yè)氣象,2015,36(1):50-56.

    Wang Jinghong, Liu Lu, Gao Feng, et al. Revision on meteorological indices of florescence frost disaster for Fuji apple in Shaanxi Province[J]. Chinese Journal of Agrometeorology, 2015, 36(1): 50-56. (in Chinese with English abstract)

    [13] 屈振江,劉瑞芳,郭兆夏,等. 陜西省蘋果花期凍害風(fēng)險(xiǎn)評(píng)估及預(yù)測(cè)技術(shù)研究[J]. 自然災(zāi)害學(xué)報(bào),2013,22(1):219-225.

    Qu Zhenjiang, Liu Ruifang, Guo Zhaoxia, et al. Study of risk assessment and prediction of apple blooming freezing injury in Shaanxi Province[J]. Journal of Natural Disasters, 2013, 22(1): 219-225. (in Chinese with English abstract)

    [14] 李健,劉映寧,李美榮,等. 陜西果樹花期低溫凍害特征及防御對(duì)策[J]. 氣象科技,2008(3):318-322.

    Li Jian, Liu Yingning, Li Meirong, et al. Low temperature and freezing injury to fruit trees at bloom stage in Shaanxi and countermeasures[J]. Meteorological Science and Technology, 2008(3): 318-322. (in Chinese with English abstract)

    [15] 柏秦鳳,王景紅,梁軼,等. 基于縣域單元的降尺度蘋果花期凍害風(fēng)險(xiǎn)區(qū)劃[J]. 中國農(nóng)學(xué)通報(bào),2013,29(16):153-158.

    Bai Qinfeng, Wang Jinghong, Liang Yi, et al. The apple flowering frost damage risk zoning down-scaling based on county level[J]. Chinese Agricultural Science Bulletin, 2013, 29(16): 153-158. (in Chinese with English abstract)

    [16] Chatrabgoun O, Karimi R, Daneshkhah A, et al. Copula-based probabilistic assessment of intensity and duration of cold episodes: A case study of Malayer vineyard region[J]. Agricultural and Forest Meteorology, 2020, 295: 108150.

    [17] 李美榮,李星敏,柏秦鳳,等. 蘋果極端氣象災(zāi)害氣溫極值的分布及重現(xiàn)期預(yù)測(cè)[J]. 干旱地區(qū)農(nóng)業(yè)研究,2012,30(3):257-261.

    Li Meirong, Li Xingmin, Bai Qinfeng, et al. Probability distribution and recurrence interval prediction of the extreme value of air temperature in extreme meteorological disasters of apples[J]. Agricultural Research in the Arid Areas, 2012, 30(3): 257-261. (in Chinese with English abstract)

    [18] Pouteau R, Rambal S, Ratte J P, et al. Downscaling MODIS-derived maps using GIS and boosted regression trees: The case of frost occurrence over the arid Andean highlands of Bolivia[J]. Remote Sensing of Environment, 2011, 115(1): 117-129.

    [19] Kurowicka D, Cooke R M. Uncertainty Analysis with High Dimensional Dependence Modelling[M]. Chichester: John Wiley & Sons, 2006.

    [20] 郭愛軍,暢建霞,王義民,等. 近50年涇河流域降雨-徑流關(guān)系變化及驅(qū)動(dòng)因素定量分析[J].農(nóng)業(yè)工程學(xué)報(bào),2015,31(14):165-171.

    Guo Aijun, Chang Jianxia, Wang Yimin, et al. Variation characteristics of rainfall-runoff relationship and driving factors analysis in Jinghe river basin in nearly 50 years [J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2015, 31(14): 165-171. (in Chinese with English abstract)

    [21] 王曉峰,張園,馮曉明,等. 基于游程理論和Copula函數(shù)的干旱特征分析及應(yīng)用[J]. 農(nóng)業(yè)工程學(xué)報(bào),2017,33(10):206-214.

    Wang Xiaofeng, Zhang Yuan, Feng Xiaoming, et al. Analysis and application of drought characteristics based on run theory and Copula function[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2017, 33(10): 206-214. (in Chinese with English abstract)

    [22] 藺文慧,宋松柏. 基于混合Copula函數(shù)的洪水聯(lián)合概率計(jì)算方法[J]. 水力發(fā)電學(xué)報(bào),2021,40(12):40-51.

    Lin Wenhui, Song Songbai. Study on calculation method of flood joint probability based on mixed Copula function[J]. Journal of Hydroelectric Engineering, 2021, 40(12): 40-51. (in Chinese with English abstract)

    [23] 鄔定榮,霍治國,王培娟,等. 陜西蘋果花期機(jī)理性預(yù)報(bào)模型的適用性評(píng)價(jià)[J]. 應(yīng)用氣象學(xué)報(bào),2019,30(5):555-564.

    Wu Dingrong, Huo Zhiguo, Wang Peijuan, et al. The applicability of mechanism phenology model to simulating apple flowering date in Shaanxi provice[J]. Journal of Applied Meteorological Science, 2019, 30(5): 555-564. (in Chinese with English abstract)

    [24] 王潤紅,茹曉雅,蔣騰聰,等. 基于物候模型研究未來氣候情景下陜西蘋果花期的可能變化[J]. 中國農(nóng)業(yè)氣象,2021,42(9):729-745.

    Wang Runhong, Ru Xiaoya, Jiang Tengcong, et al. Based on the phenological model to study the possible changes of apple flowering dates under future climate scenarios in Shaanxi province[J]. Chinese Journal of Agrometeorology, 2021, 42(9): 729-745. (in Chinese with English abstract)

    [25] 劉映寧,賀文麗,李艷莉,等. 陜西果區(qū)蘋果花期凍害農(nóng)業(yè)保險(xiǎn)風(fēng)險(xiǎn)指數(shù)的設(shè)計(jì)[J]. 中國農(nóng)業(yè)氣象,2010,31(1):125-129,136.

    Liu Yingning, He Wenli, Li Yanli, et al. Chinese Journal of Agrometeorolog[J]. A study on the risk index design of agricultural insurance on apple florescence freezing injury in Shaanxi fruit zone, 2010, 31(1): 125-129,136. (in Chinese with English abstract)

    [26] 李美榮,朱琳,杜繼穩(wěn). 陜西蘋果花期霜凍災(zāi)害分析[J]. 果樹學(xué)報(bào),2008(5):666-670.

    Li Meirong, Zhu Lin, Du Jiwen. Analysis of frost damage during apple bloom period in Shaanxi Province[J]. Journal of Natural Disasters, 2008(5): 666-670. (in Chinese with English abstract)

    [27] 王景紅. 陜西主要果樹氣候適宜性與氣象災(zāi)害風(fēng)險(xiǎn)區(qū)劃圖集[M]. 西安:陜西科學(xué)技術(shù)出版社,2012:85-95.

    [28] 韓會(huì)明,劉喆玥,劉成林,等. 基于Copula函數(shù)的贛江流域氣象干旱特征分析[J]. 水電能源科學(xué),2020,38(8):9-13.

    Han Huiming, Liu Zheyue, Liu Chenglin, et al. Analysis of meteorological drought characteristics in Ganjiang River basin based on Copula function[J]. Water Resources and Power, 2020, 38(8): 9-13. (in Chinese with English abstract)

    [29] 張卓群,馮冬發(fā),侯宇恒. 基于Copula函數(shù)的黃河流域干旱特征研究[J]. 干旱區(qū)資源與環(huán)境,2022,36(1):66-72.

    Zhang Zhuoqun, Feng Dongfa, Hou Yuheng. Study on drought characteristics in Yellow River basin based on Copula function[J]. Journal of Arid Land Resources and Environment, 2022, 36(1): 66-72. (in Chinese with English abstract)

    [30] Michele C D, Salvadori G, Vezzoli R, et al., Multivariate assessment of droughts: Frequency analysis and dynamic return period[J]. Water Resources Research, 2013, 49(10): 6985-6994.

    [31] Sklar A, Fonctions de repartition an dimensions et leurs marges[J]. Publ. inst. statist. univ. Paris, 1959, 8: 229-231.

    [32] Salvadori G. Bivariate return periods via 2-Copulas[J]. Statistical Methodology, 2004, 1(1): 129-144.

    [33] 張曉煜,馬玉平,蘇占勝,等. 寧夏主要作物霜凍試驗(yàn)研究[J]. 干旱區(qū)資源與環(huán)境,2001(2):50-54.

    Zhang Xiaoyu, Ma Yuping, Su Zhansheng, et al. Experiments on frost injuries of main crops in Ningxia Province[J]. Journal of Arid Land Resources and Environment, 2001(2): 50-54. (in Chinese with English abstract)

    [34] 王景紅,柏秦鳳,梁軼,等. 2013年陜西蘋果花期凍害氣象條件分析及受凍指標(biāo)研究[J]. 果樹學(xué)報(bào),2015,32(1):100-107,174.

    Wang Jinghong, Bai Qinfeng, Liang Yi, et al. Study on freezing index and weather conditions causing Shaanxi apple florescence freezing injury in 2013[J]. Journal of Natural Disasters, 2015, 32(1): 100-107, 174. (in Chinese with English abstract)

    Analysis of the characteristics of apple later frost risks in Shaanxi Province based on Copula functions

    Jiang Yuan1,2, Ru Xiaoya1,2, Luo Qi1,2, Li Meirong3, Wang Zhao3, Wang Jinghong3, Feng Hao2,4, Zhang Dong5, Su Baofeng6, Yu Qiang4, He Jianqiang1,2※

    (1.712100,; 2.712100,; 3.710015,; 4.,712100,; 5.,712100,; 6.,,712100,)

    Late frost is one of the most destructive meteorological disasters in the Loess Plateau of China. A great threat can be posed to the sustainable production of apples, leading to great economic losses in the apple industry. Thus, it is a high demand to explore the occurrence of late frost events for the prevention of apple late frost disasters in the decision-making of the local apple industry. In this study, the late frost return periods were investigated on the duration and severity of late frost events using the Coupla functions. The reliability of the model was then verified to analyze the characteristics of apple late frost. The study area was taken as the apple-producing areas of Shaanxi Province in western China. The meteorological datasets were collected from the seven weather stations from 1971-2018. The daily minimum temperature (min) of 0℃ was taken as the critical temperature for the occurrence of apple late frost events, in order to extract the two characteristic variables of duration and severity of late frost events. These characteristic variables of late frost events were first fitted by seven common distribution functions, respectively. Kolmogorov Smirnov (K-S) test was then carried out to verify the model. The joint distributions of late frost characteristic variables were constructed to evaluate the goodness-of-fit using six Copula functions. The occurrence probability and return period of late frost events were analyzed with the optimized Copula functions. The results showed that the severity of late frost risks generally increased from the southeast to northwest in the study area from 1971 to 2018. The optimal marginal distribution of late frost duration was in the log-normal distribution, while there was a great difference in the optimal marginal distributions of late frost severity. A significant positive correlation was found between the duration variables and the severity of late frost at each station. The joint cumulative probability increased significantly, as the severity and duration of late frost increased, but the increasing trend was much slower than before. A much more significant increase was observed in the co-occurrence return period under the same increase of univariate value, compared with the joint. The univariate return period was always between the joint and co-occurrence return period. Once the univariate return period was small enough, the optimal range of the actual univariate return period can be estimated, according to the joint and the co-occurrence return period. In general, a low probability was found in the late frost events with long duration and high severity at the weather stations in apple-producing areas in Shaanxi Province of China. However, the stations in the Yan’an area were more susceptible to the late frost events with high severity or long duration, as well as the late frost events with both high severity and long duration. Thus, more attention can be paid to late frost risks in the Yan’an area. This finding can provide a theoretical base to deal with the late frost disaster in apple production.

    meteorological; disasters; late frost; Copula function; frequency analysis; return period; apple

    10.11975/j.issn.1002-6819.2022.23.018

    S166

    A

    1002-6819(2022)-23-0170-11

    姜元,茹曉雅,羅琦,等. 基于Copula函數(shù)的陜西蘋果晚霜凍特征分析[J]. 農(nóng)業(yè)工程學(xué)報(bào),2022,38(23):170-180.doi:10.11975/j.issn.1002-6819.2022.23.018 http://www.tcsae.org

    Jiang Yuan, Ru Xiaoya, Luo Qi, et al. Analysis of the characteristics of apple later frost risks in Shaanxi Province based on Copula functions[J]. Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2022, 38(23): 170-180. (in Chinese with English abstract) doi:10.11975/j.issn.1002-6819.2022.23.018 http://www.tcsae.org

    2022-08-21

    2022-11-26

    國家重點(diǎn)研發(fā)計(jì)劃(2021YFD1900700);陜西省重點(diǎn)研發(fā)計(jì)劃重點(diǎn)產(chǎn)業(yè)創(chuàng)新鏈(群)-農(nóng)業(yè)領(lǐng)域項(xiàng)目(2019ZDLNY07-03);陜西省氣象局秦嶺和黃土高原生態(tài)環(huán)境氣象重點(diǎn)實(shí)驗(yàn)室2019年開放研究基金課題(2019Z-5);西北農(nóng)林科技大學(xué)人才專項(xiàng)資金(千人計(jì)劃項(xiàng)目);高等學(xué)校學(xué)科創(chuàng)新引智計(jì)劃(111計(jì)劃)(B12007)

    姜元,博士生,研究方向?yàn)檗r(nóng)業(yè)生態(tài)系統(tǒng)模擬。Email:yuanjiang59@nwafu.edu.cn

    何建強(qiáng),教授,博士生導(dǎo)師,研究方向?yàn)檗r(nóng)業(yè)生態(tài)系統(tǒng)模擬。Email:jianqiang_he@nwafu.edu.cn

    猜你喜歡
    霜凍歷時(shí)花期
    近50年海原縣霜凍發(fā)生特征及其對(duì)農(nóng)業(yè)的影響
    大豆:花期結(jié)莢期巧管理
    量詞“只”的形成及其歷時(shí)演變
    常用詞“怠”“惰”“懶”的歷時(shí)演變
    農(nóng)作物防御霜凍六法
    作物遭受霜凍該如何補(bǔ)救
    對(duì)《紅樓夢(mèng)》中“不好死了”與“……好的”的歷時(shí)考察
    古今字“兌”“說”“悅”“敚”歷時(shí)考察
    Current status and challenges in sentinel node navigation surgery for early gastric cancer
    農(nóng)作物的殺手——霜凍
    精品99又大又爽又粗少妇毛片| 免费观看无遮挡的男女| 精品少妇黑人巨大在线播放| 成人国语在线视频| 黄色 视频免费看| 亚洲四区av| 亚洲,一卡二卡三卡| 一二三四中文在线观看免费高清| 欧美日韩视频高清一区二区三区二| 97在线视频观看| 亚洲欧洲日产国产| 亚洲国产欧美在线一区| 国产精品国产三级国产av玫瑰| 国产一区二区在线观看av| 又大又黄又爽视频免费| 91午夜精品亚洲一区二区三区| 国产男女超爽视频在线观看| 十八禁网站网址无遮挡| 欧美人与善性xxx| 尾随美女入室| 国产一区二区激情短视频 | 欧美激情极品国产一区二区三区 | 五月玫瑰六月丁香| 国产精品国产三级国产av玫瑰| 男人操女人黄网站| 久久久亚洲精品成人影院| 国产免费现黄频在线看| 99国产综合亚洲精品| 日本-黄色视频高清免费观看| 亚洲国产精品999| 国产成人午夜福利电影在线观看| 高清视频免费观看一区二区| 精品少妇黑人巨大在线播放| 伦理电影大哥的女人| a级片在线免费高清观看视频| 欧美人与善性xxx| 五月开心婷婷网| 黄色配什么色好看| 欧美激情国产日韩精品一区| 欧美激情 高清一区二区三区| kizo精华| 久久狼人影院| 久久这里有精品视频免费| 久热这里只有精品99| 国产又色又爽无遮挡免| 亚洲在久久综合| 国产在线视频一区二区| kizo精华| 人人妻人人爽人人添夜夜欢视频| 国产伦理片在线播放av一区| 国产探花极品一区二区| 日韩,欧美,国产一区二区三区| 久久av网站| 久久女婷五月综合色啪小说| 伊人亚洲综合成人网| 最新的欧美精品一区二区| 成人黄色视频免费在线看| 少妇人妻 视频| 大片电影免费在线观看免费| 久久热在线av| 在线观看国产h片| 国产日韩一区二区三区精品不卡| 精品午夜福利在线看| 有码 亚洲区| 国产精品国产三级国产av玫瑰| 成人手机av| 69精品国产乱码久久久| 亚洲国产av新网站| 啦啦啦在线观看免费高清www| 啦啦啦在线观看免费高清www| 一本色道久久久久久精品综合| 亚洲一级一片aⅴ在线观看| 精品久久国产蜜桃| 日本av手机在线免费观看| 亚洲成人手机| 99久久综合免费| 国产1区2区3区精品| 欧美日韩国产mv在线观看视频| 精品国产国语对白av| 欧美少妇被猛烈插入视频| 狂野欧美激情性xxxx在线观看| 黄色怎么调成土黄色| 中文字幕av电影在线播放| 在线精品无人区一区二区三| 亚洲精品一二三| 看十八女毛片水多多多| 亚洲国产精品一区二区三区在线| 黑丝袜美女国产一区| 永久网站在线| 成人毛片60女人毛片免费| 亚洲av电影在线进入| 伦精品一区二区三区| 亚洲精华国产精华液的使用体验| 一个人免费看片子| 欧美日韩亚洲高清精品| 久久久久久人妻| 一级爰片在线观看| 少妇被粗大的猛进出69影院 | 久久久久精品人妻al黑| 色网站视频免费| 国产av一区二区精品久久| 日韩视频在线欧美| 侵犯人妻中文字幕一二三四区| 丝袜在线中文字幕| 天天影视国产精品| 久久av网站| 久久精品国产亚洲av涩爱| 色5月婷婷丁香| 日本黄色日本黄色录像| 中文字幕人妻丝袜制服| 午夜日本视频在线| 国产国语露脸激情在线看| 日韩三级伦理在线观看| 制服人妻中文乱码| 91精品国产国语对白视频| 成年人午夜在线观看视频| 人妻少妇偷人精品九色| 99热这里只有是精品在线观看| 丝袜脚勾引网站| 国产成人aa在线观看| 美女视频免费永久观看网站| 涩涩av久久男人的天堂| 在线观看免费视频网站a站| 色94色欧美一区二区| 伊人亚洲综合成人网| 哪个播放器可以免费观看大片| 精品一区二区免费观看| 女人精品久久久久毛片| 免费av中文字幕在线| 欧美性感艳星| 啦啦啦中文免费视频观看日本| 久久人人爽av亚洲精品天堂| av一本久久久久| 91精品国产国语对白视频| 美国免费a级毛片| 美女国产高潮福利片在线看| 久久99热6这里只有精品| 香蕉精品网在线| 蜜桃国产av成人99| 日日啪夜夜爽| 一区二区日韩欧美中文字幕 | 精品福利永久在线观看| 亚洲一区二区三区欧美精品| 制服诱惑二区| 少妇人妻 视频| 精品午夜福利在线看| 青青草视频在线视频观看| 国产在视频线精品| 欧美成人午夜精品| 蜜桃国产av成人99| 涩涩av久久男人的天堂| 国产成人一区二区在线| 黑丝袜美女国产一区| 欧美最新免费一区二区三区| 国产女主播在线喷水免费视频网站| 国产永久视频网站| 街头女战士在线观看网站| 国产精品无大码| 日韩大片免费观看网站| 看免费av毛片| 亚洲综合精品二区| 香蕉精品网在线| 久久狼人影院| 久久鲁丝午夜福利片| av不卡在线播放| 中文字幕av电影在线播放| 最近的中文字幕免费完整| 成年av动漫网址| 久久精品aⅴ一区二区三区四区 | 中文字幕免费在线视频6| 国产精品久久久av美女十八| h视频一区二区三区| 岛国毛片在线播放| 午夜福利网站1000一区二区三区| 草草在线视频免费看| 青青草视频在线视频观看| 国产精品熟女久久久久浪| 亚洲国产欧美日韩在线播放| 日韩一区二区视频免费看| 老司机影院毛片| a级毛片在线看网站| 在线观看三级黄色| 9191精品国产免费久久| 性色av一级| 两个人免费观看高清视频| 久久精品人人爽人人爽视色| 精品福利永久在线观看| 日本vs欧美在线观看视频| 中国国产av一级| 午夜激情av网站| av黄色大香蕉| 99re6热这里在线精品视频| 99国产综合亚洲精品| 视频在线观看一区二区三区| 精品一区二区免费观看| 国产精品女同一区二区软件| 最近手机中文字幕大全| 亚洲精品成人av观看孕妇| 亚洲精品一区蜜桃| 另类亚洲欧美激情| 免费高清在线观看日韩| 亚洲一区二区三区欧美精品| 成人国语在线视频| 黑人欧美特级aaaaaa片| 美国免费a级毛片| 性色av一级| 亚洲国产看品久久| 成年人免费黄色播放视频| 天天躁夜夜躁狠狠久久av| 草草在线视频免费看| 精品第一国产精品| 一级毛片我不卡| 婷婷色av中文字幕| 日韩精品免费视频一区二区三区 | av网站免费在线观看视频| 两性夫妻黄色片 | 精品一区在线观看国产| 性色avwww在线观看| 热99国产精品久久久久久7| 日本欧美国产在线视频| 国产精品久久久久成人av| 久久久久精品性色| 男女午夜视频在线观看 | 日韩三级伦理在线观看| 亚洲,一卡二卡三卡| 国产片特级美女逼逼视频| 99热这里只有是精品在线观看| 在线观看一区二区三区激情| 看十八女毛片水多多多| 国产 一区精品| 少妇的逼好多水| 久久久久精品性色| 一级爰片在线观看| 五月开心婷婷网| 只有这里有精品99| 国产男女内射视频| 久久午夜福利片| 51国产日韩欧美| 在线观看人妻少妇| 久久鲁丝午夜福利片| 日产精品乱码卡一卡2卡三| 国产精品免费大片| 亚洲综合精品二区| 久久人人97超碰香蕉20202| 免费看不卡的av| 制服丝袜香蕉在线| 亚洲国产av影院在线观看| 久久午夜综合久久蜜桃| 国产成人aa在线观看| 老女人水多毛片| 免费高清在线观看日韩| 亚洲精品久久久久久婷婷小说| 国产精品欧美亚洲77777| 人妻人人澡人人爽人人| 久久久久久人人人人人| 亚洲精品一二三| av在线观看视频网站免费| 日本猛色少妇xxxxx猛交久久| 午夜91福利影院| 九九爱精品视频在线观看| 成人亚洲精品一区在线观看| 我的女老师完整版在线观看| 亚洲四区av| 日韩电影二区| 22中文网久久字幕| 国产欧美另类精品又又久久亚洲欧美| 午夜福利影视在线免费观看| 91成人精品电影| 免费高清在线观看日韩| 国产精品国产三级国产av玫瑰| 中文字幕精品免费在线观看视频 | 少妇人妻久久综合中文| 成人无遮挡网站| 国产免费又黄又爽又色| 高清黄色对白视频在线免费看| 欧美精品国产亚洲| 水蜜桃什么品种好| 午夜免费男女啪啪视频观看| 亚洲欧美精品自产自拍| 深夜精品福利| 中文天堂在线官网| 青春草亚洲视频在线观看| 婷婷色av中文字幕| 国产一区二区激情短视频 | 久热这里只有精品99| 一级片免费观看大全| 有码 亚洲区| 一级黄片播放器| 亚洲国产精品999| 国产黄色免费在线视频| 久久久久人妻精品一区果冻| 人妻系列 视频| 伦精品一区二区三区| av免费观看日本| 久久韩国三级中文字幕| 久久精品久久久久久久性| 黄色配什么色好看| 久久精品国产自在天天线| 18禁在线无遮挡免费观看视频| 黑人巨大精品欧美一区二区蜜桃 | 国产欧美日韩综合在线一区二区| 亚洲一区二区三区欧美精品| 极品人妻少妇av视频| 各种免费的搞黄视频| 国产精品久久久久久精品电影小说| 国产日韩欧美亚洲二区| 99视频精品全部免费 在线| 久久精品久久久久久噜噜老黄| 中国国产av一级| 国产精品99久久99久久久不卡 | 五月开心婷婷网| 久久精品国产自在天天线| 夫妻午夜视频| 天堂8中文在线网| 欧美3d第一页| 亚洲,欧美精品.| 亚洲av国产av综合av卡| 欧美精品国产亚洲| 欧美国产精品va在线观看不卡| 三上悠亚av全集在线观看| 久久国产精品大桥未久av| 亚洲av男天堂| 水蜜桃什么品种好| 精品人妻熟女毛片av久久网站| 最新中文字幕久久久久| 国国产精品蜜臀av免费| 三级国产精品片| 夜夜爽夜夜爽视频| 人体艺术视频欧美日本| 亚洲欧美成人综合另类久久久| 日韩成人av中文字幕在线观看| 国产成人aa在线观看| 男女无遮挡免费网站观看| 免费女性裸体啪啪无遮挡网站| 精品一区二区三区四区五区乱码 | 亚洲人成网站在线观看播放| 日韩,欧美,国产一区二区三区| 欧美精品一区二区免费开放| 久久99热这里只频精品6学生| 三级国产精品片| 伦精品一区二区三区| av在线观看视频网站免费| 成人影院久久| 51国产日韩欧美| 97超碰精品成人国产| 国产成人精品无人区| 日本wwww免费看| 国精品久久久久久国模美| 99热6这里只有精品| 久久av网站| 成人黄色视频免费在线看| 美女xxoo啪啪120秒动态图| 九色亚洲精品在线播放| 80岁老熟妇乱子伦牲交| 精品福利永久在线观看| 自拍欧美九色日韩亚洲蝌蚪91| 亚洲欧美中文字幕日韩二区| 丝袜喷水一区| 香蕉丝袜av| 色吧在线观看| 高清欧美精品videossex| 伊人久久国产一区二区| 90打野战视频偷拍视频| 色婷婷av一区二区三区视频| www.熟女人妻精品国产 | 久久精品熟女亚洲av麻豆精品| 中文字幕人妻丝袜制服| 一本色道久久久久久精品综合| 一级毛片我不卡| videos熟女内射| 亚洲欧美一区二区三区黑人 | 国产精品 国内视频| 黑人高潮一二区| 美女内射精品一级片tv| 极品少妇高潮喷水抽搐| 99精国产麻豆久久婷婷| 色吧在线观看| 亚洲精品日本国产第一区| 国产成人午夜福利电影在线观看| 欧美日韩成人在线一区二区| 精品第一国产精品| 性色av一级| 国产免费现黄频在线看| 制服诱惑二区| 黄色怎么调成土黄色| 亚洲av综合色区一区| 亚洲国产日韩一区二区| 少妇人妻 视频| 亚洲精品av麻豆狂野| 韩国精品一区二区三区 | 免费不卡的大黄色大毛片视频在线观看| 少妇猛男粗大的猛烈进出视频| 精品视频人人做人人爽| 久久精品国产自在天天线| 日韩 亚洲 欧美在线| 久久久欧美国产精品| 哪个播放器可以免费观看大片| av国产精品久久久久影院| 免费在线观看完整版高清| 两性夫妻黄色片 | 在线观看国产h片| 日本91视频免费播放| 1024视频免费在线观看| 午夜福利视频精品| 亚洲欧美色中文字幕在线| 日韩一本色道免费dvd| 青春草亚洲视频在线观看| 男人舔女人的私密视频| 亚洲欧美精品自产自拍| 激情五月婷婷亚洲| 高清在线视频一区二区三区| 纵有疾风起免费观看全集完整版| 成人18禁高潮啪啪吃奶动态图| 18禁观看日本| 中文字幕另类日韩欧美亚洲嫩草| 伦理电影免费视频| 极品少妇高潮喷水抽搐| 人体艺术视频欧美日本| 又大又黄又爽视频免费| 国语对白做爰xxxⅹ性视频网站| 18+在线观看网站| 国产 一区精品| 毛片一级片免费看久久久久| 十八禁高潮呻吟视频| 亚洲精品自拍成人| 国产日韩一区二区三区精品不卡| 熟妇人妻不卡中文字幕| 国产av精品麻豆| 男女边吃奶边做爰视频| 超碰97精品在线观看| 国产欧美亚洲国产| 成人国语在线视频| 亚洲精品日本国产第一区| 亚洲欧美成人精品一区二区| 街头女战士在线观看网站| 亚洲高清免费不卡视频| 午夜影院在线不卡| av在线观看视频网站免费| 日本爱情动作片www.在线观看| 人妻系列 视频| 亚洲中文av在线| 国产激情久久老熟女| 亚洲精品456在线播放app| 精品少妇黑人巨大在线播放| 22中文网久久字幕| 国产男人的电影天堂91| 麻豆乱淫一区二区| 不卡视频在线观看欧美| 亚洲国产av影院在线观看| 亚洲av免费高清在线观看| 亚洲,欧美精品.| 亚洲熟女精品中文字幕| 国语对白做爰xxxⅹ性视频网站| 免费人成在线观看视频色| 婷婷色麻豆天堂久久| 久久国产亚洲av麻豆专区| 久久久久久久亚洲中文字幕| 美女大奶头黄色视频| 看免费av毛片| 国产日韩欧美亚洲二区| freevideosex欧美| 久久午夜福利片| 亚洲成av片中文字幕在线观看 | 大香蕉久久成人网| 亚洲精华国产精华液的使用体验| 免费高清在线观看视频在线观看| 亚洲精品av麻豆狂野| 免费少妇av软件| 黑丝袜美女国产一区| www.熟女人妻精品国产 | 亚洲经典国产精华液单| 99热网站在线观看| 国产日韩欧美视频二区| 全区人妻精品视频| av天堂久久9| 中文字幕亚洲精品专区| 人成视频在线观看免费观看| 午夜免费鲁丝| videos熟女内射| 国产黄色免费在线视频| 夜夜骑夜夜射夜夜干| 尾随美女入室| 亚洲欧美成人综合另类久久久| 亚洲av免费高清在线观看| www.av在线官网国产| 免费av中文字幕在线| 国产乱人偷精品视频| 这个男人来自地球电影免费观看 | 岛国毛片在线播放| av一本久久久久| 久久精品熟女亚洲av麻豆精品| 99热网站在线观看| 午夜精品国产一区二区电影| 欧美日韩一区二区视频在线观看视频在线| 久久久久视频综合| 自线自在国产av| 日本av免费视频播放| 欧美国产精品va在线观看不卡| 国产白丝娇喘喷水9色精品| 少妇的逼好多水| 精品亚洲乱码少妇综合久久| 国产一区二区在线观看日韩| www日本在线高清视频| 精品人妻熟女毛片av久久网站| 少妇的逼好多水| 一级爰片在线观看| 26uuu在线亚洲综合色| 麻豆乱淫一区二区| 一级片免费观看大全| 国产成人免费无遮挡视频| 少妇被粗大的猛进出69影院 | 女人精品久久久久毛片| 美女国产高潮福利片在线看| 插逼视频在线观看| 黄色怎么调成土黄色| 男人添女人高潮全过程视频| 又大又黄又爽视频免费| 97精品久久久久久久久久精品| 日本猛色少妇xxxxx猛交久久| 国产精品国产三级国产专区5o| 天天操日日干夜夜撸| av天堂久久9| av在线老鸭窝| 久久久久视频综合| 午夜影院在线不卡| 五月玫瑰六月丁香| 91成人精品电影| 最近最新中文字幕免费大全7| 日本91视频免费播放| xxx大片免费视频| 精品亚洲乱码少妇综合久久| 久久精品国产自在天天线| 欧美丝袜亚洲另类| 精品熟女少妇av免费看| 日韩成人av中文字幕在线观看| 欧美丝袜亚洲另类| 99re6热这里在线精品视频| 国产一区二区激情短视频 | 日韩伦理黄色片| 九色亚洲精品在线播放| 最近最新中文字幕大全免费视频 | 有码 亚洲区| 亚洲精品久久成人aⅴ小说| 美女主播在线视频| 熟女人妻精品中文字幕| 欧美精品人与动牲交sv欧美| 婷婷色综合www| 亚洲成人av在线免费| 嫩草影院入口| 亚洲成色77777| 亚洲国产最新在线播放| 90打野战视频偷拍视频| 久久韩国三级中文字幕| 精品一品国产午夜福利视频| 高清av免费在线| 久热这里只有精品99| 91久久精品国产一区二区三区| 日韩视频在线欧美| xxx大片免费视频| 精品一区二区三卡| 久久人人爽人人爽人人片va| 男女无遮挡免费网站观看| 一级a做视频免费观看| 久热这里只有精品99| √禁漫天堂资源中文www| 亚洲第一av免费看| 各种免费的搞黄视频| 男女无遮挡免费网站观看| 多毛熟女@视频| 日韩欧美一区视频在线观看| 精品一区在线观看国产| 亚洲精品一二三| 国产精品女同一区二区软件| 另类亚洲欧美激情| 黑人欧美特级aaaaaa片| 亚洲丝袜综合中文字幕| 少妇人妻 视频| 久久精品国产a三级三级三级| 亚洲成国产人片在线观看| 精品国产一区二区三区四区第35| 黑人猛操日本美女一级片| 午夜视频国产福利| 成人漫画全彩无遮挡| 亚洲久久久国产精品| 亚洲精品av麻豆狂野| 亚洲成人av在线免费| 我要看黄色一级片免费的| 欧美精品一区二区大全| 在线免费观看不下载黄p国产| 国产精品国产三级专区第一集| 999精品在线视频| 免费在线观看黄色视频的| 高清毛片免费看| 99香蕉大伊视频| 国产日韩欧美亚洲二区| 又粗又硬又长又爽又黄的视频| 欧美老熟妇乱子伦牲交| 免费少妇av软件| 一本色道久久久久久精品综合| 夜夜爽夜夜爽视频| 自拍欧美九色日韩亚洲蝌蚪91| 午夜视频国产福利| 91久久精品国产一区二区三区| 99re6热这里在线精品视频| 久久久久久久精品精品| www.av在线官网国产| 亚洲精品av麻豆狂野| 国产av精品麻豆| 精品一区二区三区视频在线| 午夜福利网站1000一区二区三区| 99热6这里只有精品| 母亲3免费完整高清在线观看 | 美女国产视频在线观看| 在线观看免费视频网站a站| 欧美少妇被猛烈插入视频|