基于LAI和VTCI及Copula函數(shù)的冬小麥單產(chǎn)估測(cè)

王鵬新 陳 弛 張樹譽(yù) 張 悅 李紅梅

(1.中國農(nóng)業(yè)大學(xué)信息與電氣工程學(xué)院， 北京 100083; 2.陜西省氣象局， 西安 710014)

0 引言

冬小麥?zhǔn)俏覈饕Z食作物之一，對(duì)冬小麥進(jìn)行大規(guī)模的長(zhǎng)勢(shì)監(jiān)測(cè)和產(chǎn)量估測(cè)，不僅能為地方政府制定生產(chǎn)計(jì)劃提供科學(xué)依據(jù)，也能為糧食安全提供保障。近年來，遙感技術(shù)以其快速、簡(jiǎn)便、宏觀、無損及客觀等優(yōu)點(diǎn)被廣泛應(yīng)用于農(nóng)業(yè)生產(chǎn)的各個(gè)環(huán)節(jié)[1]。國內(nèi)外學(xué)者通常使用可以表征作物生長(zhǎng)信息的植被指數(shù)(Vegetation index，VI)來建立與作物產(chǎn)量間的經(jīng)驗(yàn)?zāi)Ｐ停M(jìn)而進(jìn)行作物估產(chǎn)。以往的研究中常用的植被指數(shù)有歸一化植被指數(shù)(Normalized difference vegetation index，NDVI)[2]、增強(qiáng)型植被指數(shù)(Enhanced vegetation index，EVI)[3]、比值植被指數(shù)(Ratio vegetation index，RVI)[4]等。利用遙感技術(shù)獲得的葉面積指數(shù)(Leaf area index，LAI)可以反映作物光合作用、呼吸作用、蒸騰作用及碳氮循環(huán)等生物物理過程[5]，常作為作物長(zhǎng)勢(shì)監(jiān)測(cè)和產(chǎn)量估測(cè)的重要參數(shù)[6]。然而作物單產(chǎn)不僅與植被指數(shù)及其相關(guān)的長(zhǎng)勢(shì)參數(shù)有關(guān)，還與土壤水分有密切聯(lián)系。在NDVI和地表溫度(Land surface temperature，LST)的散點(diǎn)圖呈三角形區(qū)域分布的基礎(chǔ)上，文獻(xiàn)[7]提出條件植被溫度指數(shù)(Vegetation temperature condition index，VTCI)的干旱監(jiān)測(cè)方法，彌補(bǔ)了單一遙感干旱監(jiān)測(cè)指數(shù)的不足，并在干旱監(jiān)測(cè)、預(yù)測(cè)和影響評(píng)估等領(lǐng)域得到了廣泛的應(yīng)用。在此基礎(chǔ)上，綜合考慮作物長(zhǎng)勢(shì)與干旱脅迫可以更加準(zhǔn)確地描述作物生長(zhǎng)狀態(tài)，從而提高作物產(chǎn)量估測(cè)精度。

Copula函數(shù)是構(gòu)造多元聯(lián)合分布和隨機(jī)變量間相關(guān)結(jié)構(gòu)的重要工具。近年來，Copula函數(shù)以其在多變量分析中的便利性引起了學(xué)者們的廣泛關(guān)注，它能夠?qū)⒆兞康南嚓P(guān)性結(jié)構(gòu)與邊緣分布分開處理，考慮變量間的非線性與非正態(tài)關(guān)系，沒有任何限制[8-10]。VERGNI等[11]將Copula理論應(yīng)用于農(nóng)業(yè)干旱領(lǐng)域，利用Student Copula獲得了向日葵相對(duì)發(fā)病率和相對(duì)嚴(yán)重程度的聯(lián)合概率以及雙變量返回期，為干旱規(guī)劃和管理提供了有用的信息。張迎等[12]綜合降雨和徑流兩種要素，基于Archimedean Copula構(gòu)建了一種能夠綜合表征氣象干旱和水文干旱的新型綜合干旱指標(biāo)，并用其表征渭河流域的干旱演變特征。針對(duì)Copula函數(shù)構(gòu)建多元隨機(jī)變量間聯(lián)合分布時(shí)參數(shù)求解困難、不同Copula函數(shù)的參數(shù)不盡相同的問題，王鵬新等[13-14]分別利用主成分分析法(Principal component analysis，PCA)和核主成分分析法(Kernel PCA，KPCA)提取VTCI的主成分因子，基于Copula函數(shù)評(píng)估冬小麥主要生育期干旱對(duì)其產(chǎn)量的影響，驗(yàn)證了Copula函數(shù)對(duì)基于VTCI的干旱影響評(píng)估研究具有較好的適用性。然而作物的長(zhǎng)勢(shì)和產(chǎn)量是多個(gè)指標(biāo)共同作用的結(jié)果，僅靠單一的指標(biāo)進(jìn)行作物長(zhǎng)勢(shì)監(jiān)測(cè)研究難以反映不同指標(biāo)在作物生長(zhǎng)與產(chǎn)量估測(cè)過程中的相互影響[15]。因此，本文以陜西省關(guān)中平原為研究區(qū)域，選取VTCI和LAI為長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)，將PCA應(yīng)用于提取單個(gè)、兩個(gè)特征變量的主成分?jǐn)?shù)據(jù)，再結(jié)合無需求解參數(shù)的乘積Copula函數(shù)構(gòu)建綜合長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)，并建立綜合長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)與冬小麥單產(chǎn)之間的線性回歸模型，以期獲得更準(zhǔn)確的冬小麥單產(chǎn)估測(cè)方法。

1 材料與方法

1.1 研究區(qū)域概況

關(guān)中平原位于陜西省中部的渭河流域，西起寶雞大散關(guān)，東至渭南潼關(guān)，北到陜北黃土高原，南至秦嶺，坐標(biāo)為106°22′～110°24′E，33°57′～35°39′N。行政區(qū)域包括西安市、銅川市、寶雞市、咸陽市、渭南市和楊凌國家農(nóng)業(yè)高新技術(shù)產(chǎn)業(yè)示范區(qū)，其中楊凌國家農(nóng)業(yè)高新技術(shù)產(chǎn)業(yè)示范區(qū)由于面積較小，故將其劃分到咸陽市境內(nèi)。關(guān)中平原是陜西省的農(nóng)業(yè)基地，也是全國糧食主要生產(chǎn)基地之一[16]。該區(qū)域年平均氣溫在6～13℃之間，年平均降雨量在550～700 mm之間，但分布極不均勻，11月到次年5月易發(fā)生干旱[17]。受全球氣候變暖影響，農(nóng)業(yè)干旱災(zāi)害不斷加劇，由干旱造成的冬小麥減產(chǎn)現(xiàn)象時(shí)有發(fā)生[18-19]。

1.2 試驗(yàn)數(shù)據(jù)

1.2.1時(shí)間序列VTCI的生成

選取關(guān)中平原2012—2017年每年3—5月的Aqua-MODIS遙感數(shù)據(jù)的日地表反射率產(chǎn)品(MYD09GA)和日地表溫度產(chǎn)品(MYD11A1)，獲取日LST和日NDVI，通過最大值合成法分別逐像素取10 d內(nèi)所包含的多日LST和NDVI的最大值作為該像素的LST和NDVI值，即獲得旬LST和旬NDVI最大值合成產(chǎn)品。將多年某一旬的LST和NDVI最大值合成產(chǎn)品再一次運(yùn)用最大值合成技術(shù)，生成多年旬尺度LST和NDVI最大值合成產(chǎn)品。基于最小值合成技術(shù)將多年某一旬LST的最大值合成產(chǎn)品逐像素取最小值，得到多年旬尺度LST最大-最小值合成產(chǎn)品；通過多年旬尺度NDVI和LST最大值合成產(chǎn)品確定研究區(qū)域VTCI的熱邊界，通過多年旬尺度NDVI最大值合成產(chǎn)品和多年旬尺度LST最大-最小值合成產(chǎn)品確定VTCI的冷邊界。VTCI計(jì)算公式為[20-21]

(1)

其中

LSTmax(NDVIi)=a+bNDVIi

(2)

LSTmin(NDVIi)=a′+b′NDVIi

(3)

式中LST(NDVIi)——研究區(qū)域內(nèi)某一像素的NDVI值為NDVIi時(shí)的地表溫度

LSTmax(NDVIi)——研究區(qū)域內(nèi)當(dāng)NDVI值等于NDVIi時(shí)的所有像素地表溫度的最大值，稱作熱邊界

LSTmin(NDVIi)——研究區(qū)域內(nèi)當(dāng)NDVI值等于NDVIi時(shí)的所有像素地表溫度的最小值，稱作冷邊界

a、b、a′、b′——待定系數(shù)，由研究區(qū)域NDVI和LST的散點(diǎn)圖近似獲得

VTCI的取值范圍為[0,1]，其中，VTCI值越小，表示旱情越嚴(yán)重，值越大，則相反。根據(jù)越冬后冬小麥的生長(zhǎng)情況，將越冬后的生育時(shí)期劃分為返青期(3月上旬—中旬)、拔節(jié)期(3月下旬—4月中旬)、抽穗-灌漿期(4月下旬—5月上旬)和乳熟期(5月中旬—下旬)4個(gè)生育時(shí)期[22-23]。取某一生育時(shí)期內(nèi)多旬VTCI的平均值作為該生育時(shí)期的VTCI值，疊加關(guān)中平原的行政矢量圖，取各縣(區(qū))所包含像素的VTCI的平均值作為該地區(qū)該生育時(shí)期的VTCI值。

1.2.2時(shí)間序列LAI的生成

選取2012—2017年冬小麥主要生育時(shí)期的MCD15A3H產(chǎn)品，與MOD15A2(Terra-MODIS)和MYD15A2(Aqua-MODIS)產(chǎn)品相比，MCD15A3H具有較高的時(shí)間分辨率(4 d)和空間分辨率(500 m)，對(duì)于監(jiān)測(cè)農(nóng)作物的生長(zhǎng)和物候特性更為有利。原始葉面積指數(shù)產(chǎn)品由于云和大氣等因素的影響，會(huì)導(dǎo)致LAI數(shù)據(jù)出現(xiàn)驟降現(xiàn)象，從而降低準(zhǔn)確性。為了解決這個(gè)問題，應(yīng)用上包絡(luò)線Savitzky-Golay(S-G)濾波對(duì)原始時(shí)間序列LAI進(jìn)行平滑處理[24]，經(jīng)上包絡(luò)線S-G濾波平滑處理后的葉面積指數(shù)數(shù)據(jù)更加符合冬小麥的生長(zhǎng)情況。

為了使LAI和VTCI具有相同的取值范圍，將S-G濾波后的LAI進(jìn)行歸一化處理；逐像素取所在旬LAI的最大值作為該旬的LAI，取各生育時(shí)期所包含的多旬LAI的平均值作為該生育時(shí)期的LAI值，通過疊加研究區(qū)域的行政矢量圖，取各縣(區(qū))所包含像素的LAI的平均值作為該地區(qū)該生育時(shí)期的LAI值。

1.2.3冬小麥單產(chǎn)數(shù)據(jù)來源

《陜西統(tǒng)計(jì)年鑒》記錄了冬小麥的產(chǎn)量數(shù)據(jù)。因此，關(guān)中平原各縣(區(qū))2012—2017年的單產(chǎn)數(shù)據(jù)均來源于其所在市(西安市、咸陽市、寶雞市、渭南市)的統(tǒng)計(jì)年鑒。

1.3 PCA-Copula法

主成分分析是一種多元統(tǒng)計(jì)方法，可通過降維技術(shù)將多個(gè)變量轉(zhuǎn)換為少數(shù)幾個(gè)主成分，這些主成分能夠反映原始變量的大部分信息，通常表示為原始變量的線性組合。由于主成分因子是新的互相獨(dú)立的變量，因此在建立Copula函數(shù)時(shí)無需求解參數(shù)，計(jì)算簡(jiǎn)便。

1.3.1主成分分析

首先根據(jù)樣本數(shù)據(jù)構(gòu)建矩陣An×p，計(jì)算協(xié)方差矩陣Rp×p，進(jìn)而求出矩陣的特征值λi(i=1, 2,…,p)及其對(duì)應(yīng)的單位特征向量ei，各主成分的貢獻(xiàn)率Di用來反映其信息量[13]，計(jì)算公式為

(4)

最終選擇幾個(gè)主成分，即PC1,PC1,…,PCm中的m通常由主成分的累計(jì)貢獻(xiàn)率G(m)來確定，計(jì)算式為

(5)

根據(jù)經(jīng)驗(yàn)，當(dāng)G(m)大于85%時(shí)，認(rèn)為能足夠反映原來的信息，對(duì)應(yīng)的m就是選取的前m個(gè)主成分。第r個(gè)主成分的線性表示為

PCr=a1rX1+a2rX2+…+aprXp
(r=1,2,…,m)

(6)

式中PCr——第r個(gè)主成分因子

apr——所選指標(biāo)在第r個(gè)主成分線性組合中的系數(shù)

p——原變量所選指標(biāo)的個(gè)數(shù)

Xp——所選指標(biāo)(如VTCI、LAI)的值

1.3.2Copula函數(shù)

Copula函數(shù)是定義在[0, 1]上均勻分布的隨機(jī)向量的聯(lián)合分布函數(shù)。根據(jù)Sklar定理[25]，聯(lián)合分布函數(shù)與相關(guān)結(jié)構(gòu)函數(shù)之間存在一一對(duì)應(yīng)的關(guān)系，Copula函數(shù)形式可表示為

F(x1,x2,…,xj)=Cθ(F1(x1),F2(x2),…,Fj(xj))=
C(u1,u2,…,uj)

(7)

其中

uj=F(xj)

式中F——隨機(jī)變量分布函數(shù)

j——樣本容量

C——Copula函數(shù)

uj——隨機(jī)變量xj的最優(yōu)邊緣分布函數(shù)

由于采用PCA方法提取樣本數(shù)據(jù)的m個(gè)主成分因子間互無相關(guān)性，故可以通過相應(yīng)的乘積Copula獲取綜合長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)，即

(8)

式中I——基于前m個(gè)主成分確定的綜合長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)

1.3.3邊緣分布函數(shù)確定與優(yōu)度評(píng)價(jià)

由給定樣本集合求解隨機(jī)變量邊緣分布函數(shù)的方法包括參數(shù)法和非參數(shù)法。其中參數(shù)法是假定總體服從某種已知的分布，利用極大似然法估計(jì)分布函數(shù)的參數(shù)，這種方法依賴于實(shí)現(xiàn)對(duì)總體分布的假設(shè)。常用的分布線型包括皮爾遜Ⅲ型分布、正態(tài)分布、廣義極值分布、對(duì)數(shù)分布、指數(shù)分布等[26]。首先通過Kolmogorov-Smirnov(K-S)方法檢驗(yàn)樣本是否服從指定的分布，原假設(shè)為樣本服從指定分布，備擇假設(shè)為不服從指定分布。當(dāng)輸出檢驗(yàn)值h=1時(shí)，在顯著性水平下拒絕原假設(shè)；h=0時(shí)，則在顯著性水平下接受原假設(shè)。若該樣本某一主成分的多個(gè)指定分布均通過K-S檢驗(yàn)，則通過各邊緣分布的理論頻率與經(jīng)驗(yàn)頻率之間的均方根誤差(Root mean square error，RMSE)及赤池信息準(zhǔn)則值(Akaike information criterion，AIC)評(píng)估擬合效果，RMSE和AIC越小，表明擬合效果越好[27]。一維隨機(jī)分布的經(jīng)驗(yàn)頻率[28]計(jì)算公式為

(9)

式中N——樣本數(shù)

s——從最小樣本開始觀察的順序值

Qs——第s項(xiàng)的累計(jì)經(jīng)驗(yàn)頻率RMSE計(jì)算公式為

(10)

式中Qsir——第r個(gè)主成分第i個(gè)數(shù)據(jù)的經(jīng)驗(yàn)頻率

Qtir——第r個(gè)主成分第i個(gè)數(shù)據(jù)的理論頻率

n——第r個(gè)主成分因子中數(shù)據(jù)的個(gè)數(shù)，取n為144

AIC計(jì)算公式為

VAIC=nlnVRMSE+2σ

(11)

式中σ——所選模型中參數(shù)的個(gè)數(shù)

當(dāng)樣本總體分布未知或不符合常見分布的線型時(shí)，則不易用參數(shù)法對(duì)邊緣分布作出估計(jì)。非參數(shù)法則避開了線型選擇的問題，經(jīng)驗(yàn)分布與核密度估計(jì)均屬于非參數(shù)檢驗(yàn)方法，它們可以從樣本數(shù)據(jù)本身出發(fā)研究數(shù)據(jù)分布的特征，估計(jì)未知的分布函數(shù)，不受限于事先對(duì)總體分布做出假設(shè)。

2 結(jié)果與分析

2.1 冬小麥4個(gè)生育時(shí)期各長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)的主成分分析

由于VTCI和經(jīng)歸一化處理后的LAI的取值范圍均為[0,1]，因此直接通過PCA方法對(duì)原始數(shù)據(jù)進(jìn)行降維處理，根據(jù)主成分相應(yīng)累計(jì)貢獻(xiàn)率達(dá)到85%以上的原則選取主成分(表1)。VTCI第1主成分的貢獻(xiàn)率為64.78%，第2主成分貢獻(xiàn)率為25.61%，前2個(gè)主成分的累計(jì)貢獻(xiàn)率為90.39%；LAI的第1主成分貢獻(xiàn)率達(dá)到91.68%，第2主成分貢獻(xiàn)率為7.10%，前2個(gè)主成分的累計(jì)貢獻(xiàn)率為98.78%；VTCI和LAI雙變量數(shù)據(jù)、共8項(xiàng)指標(biāo)經(jīng)PCA處理后的前3個(gè)主成分的貢獻(xiàn)率分別為59.44%、14.41%和12.27%，其累計(jì)貢獻(xiàn)率為86.12%。

表1 冬小麥4個(gè)生育時(shí)期不同長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)的主成分貢獻(xiàn)率及其線性表達(dá)式Tab.1 Principal component contribution rates and linear expressions of different growth monitoring indices at four growth stages of winter wheat

2.2 主成分因子邊緣分布函數(shù)的確定

基于目前常用的分布線型，選取正態(tài)分布、廣義極值分布和對(duì)數(shù)分布分別擬合冬小麥4個(gè)生育時(shí)期各長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)的主成分因子，應(yīng)用極大似然法獲得對(duì)應(yīng)的邊緣分布函數(shù)參數(shù)，利用K-S檢驗(yàn)對(duì)各主成分的邊緣分布函數(shù)進(jìn)行擬合優(yōu)度評(píng)價(jià)(表2)?？梢钥闯?，各長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)主成分因子的分布擬合情況分為2種，指定的3種分布均未通過K-S檢驗(yàn)與僅通過正態(tài)分布和廣義極值分布的K-S檢驗(yàn)。以VTCI的第1、2主成分為例，闡述邊緣分布函數(shù)的選擇過程。VTCI第1主成分的正態(tài)分布、廣義極值分布和對(duì)數(shù)分布均未通過K-S檢驗(yàn)，因此選用非參數(shù)法擬合該主成分的邊緣分布。由圖1可見，VTCI第1主成分因子經(jīng)驗(yàn)分布和核分布估計(jì)曲線的走勢(shì)基本吻合，因此優(yōu)選更加光滑的核分布估計(jì)結(jié)果作為該主成分的邊緣分布。VTCI第2主成分的正態(tài)分布和廣義極值分布均能通過K-S檢驗(yàn)，依據(jù)邊緣分布的理論頻率與經(jīng)驗(yàn)頻率之間RMSE及AIC確定最優(yōu)邊緣分布(式(9))。經(jīng)計(jì)算，廣義極值分布的擬合結(jié)果(VRMSE=0.013，VAIC=-616.592)優(yōu)于正態(tài)分布的擬合結(jié)果(VRMSE=0.016，VAIC=-592.383)，因此VTCI的第2主成分優(yōu)選廣義極值分布。同理可以得出，LAI的第1主成分優(yōu)選廣義極值分布，第2主成分選擇核分布；VTCI和LAI的第1、2、3主成分均優(yōu)選廣義極值分布(表2)。

表2 冬小麥4個(gè)生育時(shí)期不同長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)主成分因子的K-S檢驗(yàn)值(h)及其優(yōu)選分布Tab.2 K-S test value (h) and optimal distribution of principal components of different growth monitoring indices at four growth stages of winter wheat

2.3 綜合長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)與冬小麥單產(chǎn)間的線性回歸分析

以VTCI為例，闡述綜合長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)的構(gòu)建過程。從表1中VTCI第1主成分PCV1的表達(dá)式看，返青期VTCI的系數(shù)最大，其他生育時(shí)期的系數(shù)比較小，說明返青期VTCI對(duì)PCV1的影響最大。從第2主成分PCV2的表達(dá)式看，V1的系數(shù)為負(fù)數(shù)，說明隨著返青期VTCI的減小，PCV2增大，可能是由于返青期的干旱的補(bǔ)償效應(yīng)造成的，即在冬小麥返青期發(fā)生一定程度的水分脅迫不會(huì)對(duì)作物生長(zhǎng)造成太大的影響，反而能夠提高作物在生育后期的水分利用效率，使產(chǎn)量增加[29]；PCV2在拔節(jié)期VTCI、抽穗-灌漿期VTCI、乳熟期VTCI上有均相近的正系數(shù)，說明這3個(gè)生育時(shí)期對(duì)PCV2的重要性都相似?；诖?，分別以各主成分的貢獻(xiàn)率在其累計(jì)貢獻(xiàn)率中所占的比重將主成分綜合，進(jìn)一步獲取前2個(gè)主成分的綜合線性表達(dá)式為

PCV=0.717PCV1+0.283PCV2

(12)

將VTCI前2個(gè)主成分的線性表達(dá)式(表1)代入式(12)得到綜合評(píng)價(jià)值為

PCV=0.664V1+0.246V2+0.218V3+0.213V4

(13)

由式(13)可見，返青期VTCI的系數(shù)最大，乳熟期VTCI的系數(shù)最小。但拔節(jié)期是干旱對(duì)冬小麥生長(zhǎng)影響的關(guān)鍵因素，其次為抽穗-灌漿期，返青期和乳熟期的影響相對(duì)較小[22-23]。該模型的構(gòu)建沒有考慮前2個(gè)主成分各自的分布特征，因此僅采用主成分分析得到的綜合評(píng)價(jià)值PCV存在不足，而將主成分分析與Copula函數(shù)結(jié)合的方法不受各個(gè)單因子變量邊緣分布的影響。因此，根據(jù)各主成分對(duì)應(yīng)的最優(yōu)邊緣分布特征建立聯(lián)合分布(式(8))，進(jìn)而構(gòu)建不同的綜合長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)。

由于表1中選擇各長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)的主成分個(gè)數(shù)均基于經(jīng)驗(yàn)方法，沒有理論支持該方法反映出的信息量最優(yōu)，因此利用PCA-Copula方法計(jì)算基于VTCI和LAI的包含不同主成分個(gè)數(shù)的綜合長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)，進(jìn)而與冬小麥單產(chǎn)進(jìn)行線性回歸分析(圖2、3、4)。結(jié)果表明，綜合VTCI包含第1個(gè)主成分時(shí)(記為綜合VTCI-1)與冬小麥單產(chǎn)的相關(guān)性較低(R2=0.066，P=0.002)(圖2)；綜合VTCI包含前2個(gè)主成分(記為綜合VTCI-2)與冬小麥的相關(guān)性達(dá)到極顯著水平(R2=0.246，P<0.001)，此時(shí)綜合VTCI-2中包含的主成分也符合主成分個(gè)數(shù)的選取原則。

綜合LAI包含第1個(gè)主成分(記為綜合LAI-1)與冬小麥單產(chǎn)的相關(guān)性(R2=0.567，P<0.001)優(yōu)于綜合LAI包含前2個(gè)主成分(記為綜合LAI-2)與冬小麥單產(chǎn)的相關(guān)性(R2=0.234，P<0.001)(圖3)。這是由于建立LAI的第2主成分PCL2表達(dá)式(表1)中抽穗-灌漿期LAI、乳熟期LAI系數(shù)均為負(fù)值，即PCL2在抽穗-灌漿期LAI、乳熟期LAI均有中等程度的負(fù)載荷。LAI與作物的長(zhǎng)勢(shì)及產(chǎn)量呈正相關(guān)，因此表達(dá)式中負(fù)系數(shù)與實(shí)際情況不符，加入后會(huì)帶來噪聲，使精度下降。分析表明，利用PCA-Copula方法進(jìn)行作物長(zhǎng)勢(shì)監(jiān)測(cè)和產(chǎn)量估測(cè)時(shí)，主成分的選取并非越多越好，不僅需分析所采用方法的理論基礎(chǔ)，而且需分析建立的各主成分表達(dá)式中各長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)的系數(shù)與作物長(zhǎng)勢(shì)間的關(guān)系。

綜合VTCI和LAI包含第1個(gè)主成分(記為綜合G-1)與冬小麥單產(chǎn)之間的回歸結(jié)果(R2=0.524，P<0.001)優(yōu)于綜合VTCI和LAI包含前2個(gè)主成分(記為綜合G-2)與冬小麥單產(chǎn)之間的回歸結(jié)果(R2=0.263，P<0.001)，以及綜合VTCI和LAI包含前3個(gè)主成分(記為綜合G-3)與冬小麥單產(chǎn)之間的回歸結(jié)果(R2=0.117，P<0.001)(圖4)。隨著第2、3主成分的加入，綜合VTCI和LAI與冬小麥單產(chǎn)之間的相關(guān)性逐漸降低。這是由于建立VTCI和LAI的第2主成分PCVL2表達(dá)式(表1)中拔節(jié)期VTCI、抽穗-灌漿期VTCI、乳熟期VTCI系數(shù)均為負(fù)值，其中PCVL2在乳熟期VTCI有中等程度的負(fù)載荷；建立VTCI和LAI的第3主成分PCVL3表達(dá)式(表1)中的4個(gè)負(fù)系數(shù)均與LAI相關(guān)。VTCI已被證明能夠準(zhǔn)確地反映作物水分脅迫信息，LAI與作物的長(zhǎng)勢(shì)和產(chǎn)量密切相關(guān)，PCVL2、PCVL3表達(dá)式中負(fù)系數(shù)指標(biāo)的加入會(huì)導(dǎo)致綜合值變小，精度下降，不符合實(shí)際情況。

由此可見，盡管沒有滿足主成分的累計(jì)貢獻(xiàn)率不低于85%的要求，但包含綜合長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)第1個(gè)主成分的綜合LAI-1、綜合G-1相比于綜合VTCI-2與冬小麥產(chǎn)量之間的關(guān)系更為緊密，以PCA-Copula方法為基礎(chǔ)計(jì)算作物綜合長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)時(shí)，僅僅考慮主成分分析的理論基礎(chǔ)不夠全面，更需考慮各長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)對(duì)作物長(zhǎng)勢(shì)監(jiān)測(cè)結(jié)果產(chǎn)生的實(shí)際影響，從而建立與作物長(zhǎng)勢(shì)和產(chǎn)量最密切的綜合長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)。

2.4 冬小麥單產(chǎn)估測(cè)的精度評(píng)價(jià)

基于包含最優(yōu)主成分個(gè)數(shù)的各綜合長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)分別建立估產(chǎn)模型，計(jì)算關(guān)中平原24個(gè)縣(區(qū))的估測(cè)單產(chǎn)。利用綜合VTCI-2建立的冬小麥估產(chǎn)模型估測(cè)單產(chǎn)與實(shí)際單產(chǎn)的均方根誤差為733.64 kg/hm2，相對(duì)誤差為0.05%～57.76%，平均相對(duì)誤差為13.99%；利用LAI-1建立的冬小麥估產(chǎn)模型估測(cè)單產(chǎn)與實(shí)際單產(chǎn)的均方根誤差為556.28 kg/hm2，相對(duì)誤差為0.18%～73.21%，平均相對(duì)誤差為10.31%；利用綜合G-1建立的冬小麥估產(chǎn)模型估測(cè)單產(chǎn)與實(shí)際單產(chǎn)的均方根誤差為582.81 kg/hm2，相對(duì)誤差為0.03%～54.79%，平均相對(duì)誤差為11.03%。相比之下，基于綜合VTCI-2建立的估產(chǎn)模型精度較低；基于綜合LAI-1和綜合G-1建立的估產(chǎn)模型總體估產(chǎn)誤差較小，其中綜合LAI-1建立的估產(chǎn)模型計(jì)算的估測(cè)單產(chǎn)存在部分年份或部分地區(qū)的冬小麥估測(cè)結(jié)果與實(shí)際偏差較大的情況。

為了進(jìn)一步驗(yàn)證估測(cè)模型的精度，將基于綜合LAI-1、綜合G-1建立的冬小麥估產(chǎn)模型分別用于分析2012—2017年關(guān)中平原24個(gè)縣(區(qū))估測(cè)單產(chǎn)與實(shí)際單產(chǎn)的均方根誤差和平均相對(duì)誤差(表3)。基于綜合LAI-1構(gòu)建的冬小麥單產(chǎn)回歸模型各縣(區(qū))估測(cè)單產(chǎn)與實(shí)際單產(chǎn)均方根誤差為126.74(三原縣)～1 148.87 kg/hm2(眉縣)，估測(cè)單產(chǎn)與實(shí)際單產(chǎn)的平均相對(duì)誤差以三原縣最低，為2.35%，以蒲城縣最高，為26.40%。基于綜合G-1構(gòu)建冬小麥單產(chǎn)回歸模型各縣(區(qū))估測(cè)單產(chǎn)與實(shí)際單產(chǎn)均方根誤差為125.91(大荔縣)～1 061.63 kg/hm2(閻良區(qū))，估測(cè)單產(chǎn)與實(shí)際單產(chǎn)的平均相對(duì)誤差以大荔縣最低，為2.56%，以合陽縣最高，為23.62%。上述結(jié)果表明，基于綜合LAI-1建立估產(chǎn)模型的決定系數(shù)(R2)雖略高于基于綜合G-1建立的估產(chǎn)模型，但具體到各縣(區(qū))的估產(chǎn)時(shí)，誤差范圍大于綜合G-1所建立的估產(chǎn)模型，說明利用單一LAI構(gòu)建的綜合LAI得到的冬小麥估產(chǎn)結(jié)果容易存在偏差。相對(duì)于單一指標(biāo)構(gòu)建的綜合VTCI或綜合LAI，基于雙變量構(gòu)建的綜合干旱指標(biāo)(綜合VTCI和LAI)得到的冬小麥估產(chǎn)結(jié)果更加全面可靠，可用于研究區(qū)域冬小麥單產(chǎn)的估測(cè)。

表3 關(guān)中平原各縣(區(qū))2012—2017年冬小麥估測(cè)單產(chǎn)與實(shí)際單產(chǎn)對(duì)比Tab.3 Comparison of estimated and actual winter wheat yields for each county (district) in Guanzhong Plain from 2012 to 2017

3 討論

已有研究表明Copula函數(shù)對(duì)于基于VTCI的干旱影響評(píng)估具有良好的適用性[13-14]，LAI可以反映作物在不同時(shí)間階段生長(zhǎng)發(fā)育的動(dòng)態(tài)特征和健康狀況，是表征作物長(zhǎng)勢(shì)與進(jìn)行產(chǎn)量預(yù)報(bào)的重要參數(shù)。鑒于此，本文基于冬小麥越冬后主要生育期VTCI、LAI為長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)，嘗試?yán)肞CA-Copula法分別建立冬小麥主要生育期單變量(綜合VTCI或綜合LAI)、雙變量(綜合VTCI和LAI)的估產(chǎn)模型。由主成分分析提取出的少數(shù)的幾個(gè)主成分因子幾乎保留了原始變量的全部信息，在大多實(shí)際研究工作中，選取主成分的個(gè)數(shù)有多種方法，常使用的原則是主成分的特征值大于1或累計(jì)貢獻(xiàn)率大于85%。它的優(yōu)勢(shì)在于計(jì)算簡(jiǎn)單，且對(duì)多數(shù)情況適用。但此方法沒有理論支持，屬于經(jīng)驗(yàn)性方法?；诖耍疚氖紫雀鶕?jù)主成分的累計(jì)貢獻(xiàn)率不小于85%的選取原則，在此基礎(chǔ)上嘗試令結(jié)合Copula獲得的綜合長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)包含不同的主成分個(gè)數(shù)，對(duì)比它們與冬小麥單產(chǎn)之間的擬合效果。

綜合LAI、綜合VTCI和LAI均包含第1個(gè)主成分時(shí)與單產(chǎn)間的擬合效果更好，這是由于LAI的第2個(gè)主成分表達(dá)式、VTCI和LAI的后2個(gè)主成分表達(dá)式中存在載荷量較大的負(fù)系數(shù)。VTCI和LAI均與冬小麥的長(zhǎng)勢(shì)與產(chǎn)量密切相關(guān)，兩者的值越大，作物長(zhǎng)勢(shì)越好、產(chǎn)量超高，因此負(fù)系數(shù)的出現(xiàn)不符合實(shí)際情況，加入后易導(dǎo)致精度降低。這一規(guī)律表明，開展多指標(biāo)的主成分分析時(shí)，不僅要考慮主成分的累計(jì)貢獻(xiàn)率，而且更要考慮各指標(biāo)的系數(shù)與作物長(zhǎng)勢(shì)間的關(guān)系。

利用綜合LAI構(gòu)建的估產(chǎn)模型雖然模型的擬合精度較高，但只考慮單一指標(biāo)對(duì)冬小麥產(chǎn)量進(jìn)行估測(cè)結(jié)果存在偏差。具體到每一個(gè)縣時(shí)，利用雙變量指標(biāo)綜合G-1構(gòu)建的冬小麥單產(chǎn)模型得到的估測(cè)單產(chǎn)誤差范圍比基于單變量構(gòu)建的模型更集中，結(jié)果更加可靠。研究成果表明，使用BP神經(jīng)網(wǎng)絡(luò)和IPSO-BP神經(jīng)網(wǎng)絡(luò)構(gòu)建雙變量(VTCI和LAI)冬小麥估產(chǎn)模型的R2分別為0.310、0.342[30]。相比之下，利用綜合G-1構(gòu)建估產(chǎn)模型的R2更高，進(jìn)一步說明基于PCA-Copula方法構(gòu)建的雙變量估產(chǎn)模型精度較高，可用于研究區(qū)域冬小麥單產(chǎn)估測(cè)。

4 結(jié)論

(1)采用PCA-Copula法分別對(duì)關(guān)中平原2012—2017年24個(gè)縣(區(qū))冬小麥主要生育時(shí)期單變量(VTCI或LAI)、雙變量(VTCI和LAI)進(jìn)行處理，獲得綜合VTCI、綜合LAI、綜合VTCI和LAI。結(jié)果表明，綜合VTCI包含前2個(gè)主成分時(shí)與冬小麥單產(chǎn)的擬合精度優(yōu)于包含第1個(gè)主成分的情況；綜合LAI包含第1個(gè)主成分時(shí)與單產(chǎn)的擬合精度高于綜合LAI包含前2個(gè)主成分；綜合VTCI和LAI包含第1個(gè)主成分時(shí)比包含前2個(gè)、3個(gè)主成分時(shí)與單產(chǎn)間的擬合效果好。

(2)綜合LAI-1與冬小麥單產(chǎn)間的回歸模型擬合效果最好(R2=0.567，P<0.001)，綜合G-1與冬小麥單產(chǎn)間的回歸模型略次之(R2=0.524，P<0.001)，綜合VTCI-2構(gòu)建的估產(chǎn)模型擬合效果相對(duì)較差(R2=0.246，P<0.001)。具體到各縣(區(qū))的單產(chǎn)估測(cè)時(shí)，盡管綜合LAI-1與單產(chǎn)間的相關(guān)性高于綜合G-1與單產(chǎn)間的相關(guān)性，但基于綜合LAI-1的估測(cè)單產(chǎn)相比于基于綜合G-1的估測(cè)單產(chǎn)與實(shí)際單產(chǎn)之間的誤差范圍更大，說明基于單一LAI構(gòu)建的綜合長(zhǎng)勢(shì)監(jiān)測(cè)指標(biāo)反映的冬小麥長(zhǎng)勢(shì)與產(chǎn)量信息不夠全面。因此，基于雙變量綜合VTCI和LAI構(gòu)建的冬小麥估產(chǎn)模型更加全面可靠。