基于貝葉斯統(tǒng)計的烤煙葉片臨界氮稀釋曲線構(gòu)建

摘要： 【目的】比較基于貝葉斯統(tǒng)計構(gòu)建的烤煙葉片臨界氮稀釋曲線與傳統(tǒng)方法的區(qū)別，以確定其是否可用于氮營養(yǎng)診斷，旨在簡化和準確進行烤煙氮營養(yǎng)診斷。【方法】基于兩年、3 個氮肥用量田間試驗，分析了不同氮肥施用量下移栽烤煙不同生長天數(shù)的葉片干物質(zhì)積累量、氮含量。利用各取樣時期最大葉片干物質(zhì)重，將所有處理樣本分為生長受氮限制組和不受氮限制組。然后，分別應(yīng)用貝葉斯統(tǒng)計和傳統(tǒng)兩步法建立葉片臨界氮稀釋曲線，使用 R 語言中rjags 軟件包實現(xiàn)馬爾可夫鏈蒙特卡洛 （MCMC） 算法，分析臨界氮稀釋曲線模型參數(shù)A1和A2 的后驗分布特征，分析曲線對受氮限制和不受氮限制的區(qū)分度，擬合計算葉片臨界氮濃度和氮營養(yǎng)指數(shù)，并比較其相對于實際觀測值的差異。【結(jié)果】施用氮肥顯著增加了烤煙葉片干物質(zhì)積累量，不同氮肥處理間差異顯著，烤煙葉片氮濃度隨烤煙生長進程而降低。參數(shù)A1 和A2 的95% 后驗分布分別為2.58～2.94 和0.13～0.33；擬合的烤煙葉片臨界氮稀釋曲線不確定性水平 （95% 可信區(qū)間的寬度） 隨著葉片干物重的增加先減少后增加，曲線的不確定性水平為0.16%～0.70%。基于貝葉斯統(tǒng)計的烤煙葉片臨界氮稀釋曲線對氮限制組和非氮限制組的區(qū)分度為71%，優(yōu)于傳統(tǒng)兩步法構(gòu)建曲線的區(qū)分度；兩種方法擬合的葉片臨界氮濃度與實際臨界氮濃度相關(guān)度基本一致；基于貝葉斯統(tǒng)計構(gòu)建的曲線擬合計算的氮營養(yǎng)指數(shù)略大于兩步法，通過兩種方法擬合計算的氮營養(yǎng)指數(shù)呈高度線性相關(guān)，決定系數(shù)R2=0.96，標準化均方根誤差 n-RMSE 為 6%，穩(wěn)定度較高。【結(jié)論】基于貝葉斯統(tǒng)計構(gòu)建的烤煙葉片臨界氮稀釋曲線為Nc=2.74×LDM?0.22。貝葉斯曲線在擬合臨界氮濃度的效果上，與兩步法曲線差異不大，但較兩步法曲線可以更好的區(qū)分生長受氮限制組與不受氮限制組，且計算的氮營養(yǎng)指數(shù)與兩步法曲線的線性變異決定系數(shù)R2 高達0.96，因此，貝葉斯曲線可更簡單、準確地用于評價烤煙氮營養(yǎng)狀況。

關(guān)鍵詞： 烤煙；臨界氮稀釋曲線；貝葉斯統(tǒng)計；不確定性；臨界氮濃度；氮營養(yǎng)指數(shù)

氮作為植物必需的營養(yǎng)元素，影響作物產(chǎn)量與品質(zhì)[1]，烤煙的產(chǎn)量與葉片氮含量呈正相關(guān)，而品質(zhì)與葉片氮含量呈負相關(guān)[2?3]，因此，協(xié)調(diào)產(chǎn)量與品質(zhì)的矛盾，需要將葉片氮含量控制在適宜的范圍。作物的氮素狀況可以用氮營養(yǎng)指數(shù)來判斷，氮營養(yǎng)指數(shù)為實際氮濃度與臨界氮濃度的比值[4]，根據(jù)臨界氮濃度構(gòu)建的臨界氮稀釋曲線是農(nóng)業(yè)生產(chǎn)中常用的氮營養(yǎng)狀況診斷工具[5]。自從Justes 等[6]在1994 年提出臨界氮稀釋曲線概念以來，研究者提出了多種方法來估計臨界氮稀釋曲線參數(shù)[7?8]，并建立了小麥[9?11]、水稻[12?13]、馬鈴薯[14]、番茄[15?16]、玉米[17?19]和冬油菜[20]等作物基于整株或不同組織器官[21]的臨界氮稀釋曲線來評估植物氮營養(yǎng)狀況。常規(guī)的臨界氮稀釋曲線構(gòu)建首先要根據(jù)觀測數(shù)據(jù)確定臨界氮濃度，然后采用方程 （Nc=A1×DM?A2） 擬合臨界氮稀釋曲線[5]，由于計算方法未考慮臨界氮濃度的不確定性，擬合參數(shù)的誤差較大[22?23]。

Makowski 等[24]首次應(yīng)用貝葉斯統(tǒng)計理論和馬爾科夫鏈蒙特卡洛 （MCMC） 不確定性分析方法構(gòu)建作物的臨界氮稀釋曲線。該方法不需要首先確定臨界氮濃度，可一步直接擬合原始觀測數(shù)據(jù)，因此該方法可以在觀測數(shù)據(jù)較少的情況下實現(xiàn)臨界氮稀釋曲線的構(gòu)建。隨后的一些研究應(yīng)用貝葉斯方法建立了番茄[22]、馬鈴薯[23]、小麥[25]和玉米[26]等作物的臨界氮稀釋曲線，并對臨界氮稀釋曲線的不確定性進行了分析。隨著植物生物量的增加，臨界氮稀釋曲線不確定性減小，且參數(shù) A1 （干物重為1 t/hm2時的臨界氮濃度值） 的不確定性明顯小于參數(shù) A2 （稀釋曲線的斜率，隨著生物量增加作物體內(nèi)氮濃度下降的速率）[ 2 5 ]；數(shù)據(jù)集的特征也影響構(gòu)建的臨界氮稀釋曲線，高氮條件下的觀測數(shù)據(jù)偏少和作物生育前期的觀測數(shù)據(jù)偏少，均可導(dǎo)致臨界氮稀釋曲線的不確定性增加[27]，增加試驗次數(shù)可以降低臨界氮稀釋曲線的不確定性[27]，增加觀測數(shù)據(jù)采集次數(shù)比增加氮肥處理能更有效地降低參數(shù)估計的不確定性[28]。

與其他農(nóng)作物相比，烤煙群體密度較小[29]，單位面積的干物質(zhì)積累量偏低，導(dǎo)致兩步法構(gòu)建的臨界氮稀釋曲線變異較大，應(yīng)用效果較差[30]，且地上部臨界氮稀釋曲線的應(yīng)用效果明顯差于葉片臨界氮稀釋曲線[30]。針對這一問題，本研究利用2 年的田間試驗數(shù)據(jù)，應(yīng)用貝葉斯統(tǒng)計和不確定性分析方法構(gòu)建烤煙葉片臨界氮稀釋曲線，并對擬合的臨界氮濃度與氮營養(yǎng)指數(shù)的可靠性進行了驗證，以期為烤煙氮肥施用提供較為可靠且簡單可行的診斷方法。

1 材料與方法

1.1 試驗地點

試驗于2022—2023 年在山東省濰坊市臨朐縣進行，該地區(qū)屬溫帶季風氣候，整個烤煙生育季內(nèi) （5—9 月），多年平均氣溫為23.92℃，降水量為432.58 mm，蒸發(fā)量為930.88 mm，日照時數(shù)為937.54 h。2022 年度烤煙生育期間降水量和日均溫分別為931.7 mm 和24.7°C，2023 年度烤煙生育期間降水量和日均溫分別為411.4 mm 和26.5°C。

供試土壤為褐土，兩年試驗地點的土壤基礎(chǔ)理化性質(zhì)見表1。

1.2 試驗設(shè)計

供試烤煙品種為“中川208”。

2022 年進行了試驗1，設(shè)6 個氮肥施用量處理，分別為N1 （0 kg/hm2）、N2 （30 kg/hm2）、N3 （60kg/hm2）、N4 （90 kg/hm2）、N5 （120 kg/hm2） 和N6 （150kg/hm2）；移栽時間為2022?05?04。

2023 年進行了試驗2 和試驗3，均設(shè)5 個氮肥用量處理，分別為N1 （0 kg/hm2）、N2 （45 kg/hm2）、N3 （90 kg/hm2）、N4 （135 kg/hm2）、N5 （180 kg/hm2），試驗2 與試驗3 移栽時間分別為2023?04?30 和2024?05?10。

3 個試驗均采用隨機區(qū)組排列，重復(fù)3 次，小區(qū)面積為96 m2，行距為1.2 m，株距為0.45 m。種植密度為1.8 萬棵/hm2。除氮肥外，所有處理均底施P2O5 90 kg /hm2、K2O 270 kg /hm2，供試化肥為硝酸銨鈣 （N 15.5%）、鈣鎂磷肥 （P2O5 12%）、硫酸鉀 （K2O50%）。肥料全部作為基肥，于移栽前5 天一次性條施于煙壟的壟底。灌溉、病蟲害防治等農(nóng)藝操作同當?shù)爻Ｒ?guī)生產(chǎn)。

1.3 取樣與樣品測定

試驗1 于移栽后47、61、71、83、93 天進行取樣，試驗2 于移栽后43、58、65、72、80、88 天進行取樣，試驗3 于移栽后48、53、60、67、77、85 天進行取樣。每次取樣時，每個小區(qū)選擇長勢一致的煙株3 株，采集所有葉片，105℃ 殺青30 min，65℃ 烘干至恒重，記錄葉片干物重，研磨過0.42 mm篩，經(jīng)H2SO4–H2O2 法聯(lián)合消煮后，采用流動分析儀測定葉片氮濃度。

1.4 數(shù)據(jù)分析與處理

1.4.1 兩步法擬合烤煙葉片臨界氮稀釋曲線 參考Justes 等[6]提出的計算方法，分兩步擬合烤煙葉片臨界氮稀釋曲線，第一步為計算臨界氮濃度，將每次取樣獲得的各處理葉片的葉片干物重進行方差分析，并用LSD 法進行多重比較。將葉片干物重無顯著差異的各處理定義為不受氮限制組，葉片干物量顯著小于最大干物重的各處理定義為受氮限制組。將受氮限制組各處理的葉片干物重與氮濃度進行線性擬合，對不受氮限制組的干物重取平均值，用平均值做垂直于橫軸的垂線，兩條直線的交點即為葉片臨界氮濃度。第二步為擬合臨界氮稀釋曲線，將葉片臨界氮濃度按冪函數(shù)擬合得到葉片臨界氮稀釋曲線，曲線方程為：

1.4.2 應(yīng)用貝葉斯統(tǒng)計擬合烤煙葉片臨界氮稀釋曲線 參考Makowski 等[24]提出的計算方法，采用貝葉斯統(tǒng)計理論擬合烤煙葉片臨界氮稀釋曲線方程，使用 R 包 rjags 實現(xiàn)馬爾可夫鏈蒙特卡洛 （MCMC） 算法來估計模型參數(shù)的后驗分布，其中擬合曲線時所用的參數(shù)A1、A2 的先驗區(qū)間分別為2～6、0～0.5。首先進行5 萬次迭代兩次，在驗證算法的收斂性后，再次運行 10 萬次迭代，得到參數(shù)A1、A2 的后驗分布結(jié)果樣本，計算參數(shù)A1、A2 的平均值、中位數(shù)和95% 可信區(qū)間。

1.4.3 烤煙氮營養(yǎng)指數(shù)計算 氮營養(yǎng)指數(shù)通常被用來定量植物體內(nèi)的氮素狀況，其計算公式為：

1.4.4 模型驗證 為了驗證模型的擬合度和精度，本研究通過計算均方根誤差 （RMSE） 和標準化均方根誤差 （n-RMSE） 來校驗?zāi)Ｐ停?/p>

2 結(jié)果與分析

2.1 施氮量對烤煙葉片干物質(zhì)積累的影響

如表2 所示，烤煙葉片干物質(zhì)積累量隨著生育進程呈逐漸上升趨勢，不同試驗、施氮水平和取樣時期，烤煙葉片干物質(zhì)積累量在0.20～4.05 t/hm2；各取樣時期烤煙葉片干物質(zhì)積累量均隨著施氮量的增加表現(xiàn)為先增加后保持穩(wěn)定的趨勢。參考何仲秋等[30]提出的方法，將不同施氮處理區(qū)分為受氮限制組 （虛線框內(nèi)） 和不受氮限制組 （虛線框外）。根據(jù)Justes[6]臨界氮濃度的定義：干物質(zhì)積累量較小，所有處理均為不受氮限制組，試驗1 移栽后47 和61天，試驗2 移栽后43 天，試驗3 移栽后48 天的所有處理均為不受氮限制組。根據(jù)多重比較的結(jié)果，試驗1 中，N1、N2 處理為受氮限制組，N4、N5、N6 處理為不受氮限制組；N3 處理在移栽后93 天時為受氮限制組，在移栽后71、83 天時為不受氮限制組。試驗2 中，N1、N2 處理為受氮限制組，N5 處理為不受氮限制組；N3 處理在移栽后65、72 天時為受氮限制組，在移栽后58、80、88 天時為不受氮限制組，N4 處理在移栽后65 天時為受氮限制組。試驗3 中，N1、N2 處理為生長受氮限制組，N4、N5 處理為不受氮限制組；N3 處理在移栽后53、60、77 天為受氮限制組，在移栽后67、85 天時為不受氮限制組。

2.2 基于貝葉斯統(tǒng)計構(gòu)建的烤煙葉片臨界氮稀釋曲線

隨著葉片干物重的增加，葉片氮濃度被稀釋，根據(jù)Makowski 等[ 2 4 ]提出的理論模型，使用 R 包rjags 實現(xiàn)的馬爾可夫鏈蒙特卡洛算法 （MCMC） 來估計方程 （1） 中參數(shù)A1、A2 的后驗分布，得到參數(shù)A1、A2 后驗分布頻率 （圖1），參數(shù)A1 后驗分布平均值、中值分別2.74 和2.73，分布范圍為2.58～2.94，參數(shù)A2 后驗分布平均值、中值分別為0.22 和0.21，分布范圍為0.13～0.33，其中參數(shù)A2 的變異度顯著大于參數(shù)A1 的變異度。

綜合3 個試驗數(shù)據(jù)，基于貝葉斯統(tǒng)計構(gòu)建的葉片臨界氮稀釋曲線 （簡稱貝葉斯曲線，見圖2）。圖中實線為根據(jù)參數(shù)A1 和A2 后驗分布的平均值 （2.74、0.22），帶入方程 （1） 構(gòu)建的貝葉斯曲線；兩條虛線間為貝葉斯曲線正負95% 置信區(qū)間寬度，定義為貝葉斯曲線的不確定性水平[24]。生物量水平?jīng)Q定了貝葉斯曲線的不確定性，貝葉斯曲線的不確定性水平隨著葉片干物重的增加先減少后增加，在葉干物重為1.8t/hm2 時降到最低，在葉片干物重大于0.5 t/hm2 時，貝葉斯曲線的不確定性水平為0.16%～0.70%。

2.3 基于貝葉斯統(tǒng)計與兩步法構(gòu)建的烤煙葉片臨界氮稀釋曲線模型的比較

根據(jù)兩步法構(gòu)建了烤煙葉片臨界氮稀釋曲線 （簡稱兩步法曲線，見圖3）。其具體的擬合方程為Nc =3.52×LDM?0.34，方程決定系數(shù)為0.35。兩步法曲線與貝葉斯曲線存在較大的差異，兩步法曲線 （藍色曲線） 在貝葉斯曲線 （紅色曲線） 之上。在葉片干物重小于1 t/hm2 時，兩步法曲線計算的臨界氮濃度較貝葉斯曲線計算的臨界氮濃度高28%～59%；在葉片干物重為1～2 t/hm2 時，兩步法曲線計算的臨界氮濃度較貝葉斯曲線計算的臨界氮濃度高20%～27%；在葉片干物重為2～4 t/hm2 時，兩步法曲線計算的臨界氮濃度較貝葉斯曲線計算的臨界氮濃度高12%～19%。

通過與表2 中生長受氮限制和不受氮限制處理組數(shù)的結(jié)果進行比對，進一步分析了兩步法曲線和貝葉斯曲線能否區(qū)分氮限制組和不受氮限制組。受氮限制組數(shù)據(jù)大部分在曲線以下，而不受氮限制組數(shù)據(jù)大部分在曲線以上。由圖3 可知，對于兩步法曲線，3 個試驗中，受氮限制組數(shù)量為32 個，與表2中的數(shù)量 （33個） 相當，而不受氮限制組數(shù)量僅為12 個，占表2 中觀測數(shù)量 （35 個） 的34%，兩步法曲線總區(qū)分度為65%。基于貝葉斯曲線計算的3 個試驗中受氮限制組有26 個，占觀測數(shù)據(jù) （33 個） 的79%，生長不受氮限制組數(shù)為22 個，占觀測結(jié)果（35 個） 的63%，貝葉斯曲線總的區(qū)分度為71%。貝葉斯曲線較兩步法曲線，可以更好地區(qū)分烤煙生長受氮限制組與不受氮限制組。

2.4 模型擬合的臨界氮濃度與氮營養(yǎng)指數(shù)的驗證

在兩步法構(gòu)建葉片臨界氮稀釋曲線時，首先將觀測數(shù)據(jù)分為氮受限制組和不受氮限制組，根據(jù)兩組數(shù)據(jù)計算獲得每次取樣的最大葉片干物重和葉片臨界氮濃度，該臨界氮濃度為實際臨界氮濃度。每次取樣的最大葉片干物重，分別帶入兩步法曲線和貝葉斯曲線中，計算得到的臨界氮濃度為擬合臨界氮濃度。兩步法曲線和貝葉斯曲線擬合的臨界氮濃度與實際臨界氮濃度相關(guān)性如圖4 所示，其中兩步法曲線擬合的臨界氮濃度與實際的臨界氮濃度間線性方程的決定系數(shù)為0.35，均方根誤差 RMSE 為0.57，標準化均方根誤差 n-RMSE 為 20%；貝葉斯曲線擬合的臨界氮濃度與實際臨界氮濃度間線性方程的決定系數(shù)為0.36，均方根誤差 RMSE 為 0.77，標準化均方根誤差 n-RMSE 為 27%。貝葉斯曲線在擬合臨界氮濃度的效果上，與兩步法曲線差異不大。

分析比較了根據(jù)兩步法曲線和貝葉斯曲線計算得到的氮營養(yǎng)指數(shù)之間的關(guān)系，如圖5 所示。氮營養(yǎng)指數(shù)根據(jù)方程 （2） 進行計算，其中Na 為每次取樣時，實際測定的烤煙葉片氮濃度；Nc 為根據(jù)每次取樣時的葉片干物重，代入兩步法曲線和貝葉斯曲線計算得到的臨界氮濃度。貝葉斯曲線計算得到的氮營養(yǎng)指數(shù)略大于兩步法曲線計算得到的氮營養(yǎng)指數(shù)，增大的幅度為4% 左右，兩者呈線性相關(guān)，決定系數(shù)R2 為0.96，均方根誤差 RMSE 為 0.17，標準化均方根誤差 n-RMSE 為6%。表明本研究中貝葉斯曲線計算得到的氮營養(yǎng)指數(shù)進行葉片氮營養(yǎng)診斷時，可完全代替兩步法曲線。

3 討論

3.1 烤煙葉片臨界氮稀釋曲線的特征與不確定性

臨界氮稀釋曲線是評價作物氮營養(yǎng)水平的重要方法，可用于作物氮營養(yǎng)快速診斷[15， 17， 31]，為及時調(diào)整施肥方案及實現(xiàn)節(jié)肥增效提供理論依據(jù)[9?10]。本研究中構(gòu)建的貝葉斯曲線擬合的臨界氮濃度與實際臨界氮濃度存在一定差異 （圖4），原因可能是烤煙的氮稀釋特征與其他作物存在明顯的不同，綜合本研究和以往烤煙葉片臨界氮稀釋曲線[30?31]結(jié)果發(fā)現(xiàn)，烤煙葉片干物質(zhì)積累量每增加1 t/hm2，臨界氮濃度下降0.4%～0.5%，而小麥[11]、水稻[32]和玉米[33?34]等糧食作物中，葉片干物質(zhì)積累量每增加1 t/hm2，臨界氮濃度下降0.3%～0.4%，這可能與不同作物的生長特性有關(guān)，在煙草整個生育期內(nèi) （一般在120 天左右），葉片干物重可增加至4 t/hm2 左右[30]，而整株干物重僅能增加至8 t/hm2 左右[35?36]，糧食作物在整個生育期內(nèi)可實現(xiàn)葉片干物重增加至5～6 t/hm2 左右[21]，地上部干物重可增加至13～22 t/hm2 左右[21]，相同時間內(nèi)干物質(zhì)增加速率的偏低可能導(dǎo)致氮稀釋特征不明顯，同時在圖3 中也可以看到，氮限制組和非限制組數(shù)據(jù)為交錯分布，這也證明了烤煙臨界氮稀釋曲線構(gòu)建效果會相對較差，限制了其實際應(yīng)用。基于貝葉斯統(tǒng)計構(gòu)建的烤煙葉片臨界氮稀釋曲線可部分解決該問題，貝葉斯曲線擬合的臨界氮濃度雖然與實際臨界氮濃度仍有一定誤差，但在圖3 中可以觀察到，其對氮限制組和非氮限制組的區(qū)分效果明顯優(yōu)于兩步法曲線，同時其計算的氮營養(yǎng)指數(shù)較兩步法曲線計算的氮營養(yǎng)指數(shù)略大，并線性相關(guān)，因此證明貝葉斯曲線具有更高的可應(yīng)用價值。

貝葉斯統(tǒng)計是利用現(xiàn)有數(shù)據(jù)更新模型參數(shù)先驗分布的不確定性[37]，推導(dǎo)更加確定的參數(shù)后驗分布范圍。Justes 等[6]開發(fā)的兩步法構(gòu)建臨界氮稀釋曲線，需獲得一系列臨界氮濃度點再擬合曲線，曲線擬合過程中未再考慮臨界氮濃度的回歸誤差[25]。基于貝葉斯統(tǒng)計新方法不需要區(qū)分作物生長是否受氮限制，也不需要計算臨界氮濃度，可以直接估計所有不確定性并將其表示為概率分布[24]。在本研究中，貝葉斯曲線的不確定性水平取決于烤煙葉片干物重水平的高低，這與Makowski 等[24]研究結(jié)果一致，整個貝葉斯曲線的不確定性水平在0.16%～0.70%，這與構(gòu)建的其它作物臨界氮稀釋曲線的不確定性水平基本一致[22， 25]，也進一步證明了貝葉斯曲線具有更高的可應(yīng)用價值。

3.2 不同烤煙葉片臨界氮濃度稀釋曲線的差異

由于作物的組織結(jié)構(gòu)、形態(tài)、生理特征不同，不同作物所構(gòu)建的臨界氮濃度稀釋曲線具有差異，模型參數(shù)值具有差異[8]。不同的參數(shù)反映出品種、環(huán)境等對烤煙生長過程中生物量和氮濃度相對變化的綜合效應(yīng)。因此，需要評估臨界氮濃度稀釋曲線的擬合參數(shù)，否則會增加作物氮素營養(yǎng)診斷誤差和施肥管理的決策風險[38?39]。在本研究結(jié)果中，擬合曲線參數(shù) A1，A2 分別為2.74 和0.22，與東南煙稻區(qū)[30]構(gòu)建的葉片臨界氮稀釋曲線參數(shù)相比 （A1=3.23，A2=0.47） 相對較低。不同區(qū)域之間的氣候特征和土壤供氮能力可能是導(dǎo)致A1 值存在差異的原因[35]，本試驗中土壤有機質(zhì)含量較低，土壤供氮能力較弱，所以該區(qū)域烤煙形成相同的生物量僅需要較低的臨界氮濃度。參數(shù) A2 指植株氮濃度隨生物量增加的稀釋過程，是氮濃度稀釋的具體模式。有研究表明，擬合曲線參數(shù)A2 的變異程度高于參數(shù) A1，基因型×環(huán)境×管理對其綜合的影響也大于參數(shù) A1[25]，同時貝葉斯方法與傳統(tǒng)方法構(gòu)建的臨界氮稀釋曲線參數(shù)也會由于構(gòu)建方法不同導(dǎo)致差異，因此，不同煙區(qū)間 A2值的差異原因需進一步深入探究。

3.3 應(yīng)用貝葉斯理論構(gòu)建臨界氮稀釋曲線的可行性

Lemaire 等[4]提出臨界氮濃度稀釋曲線的方法，通過擬合所有的臨界氮濃度，來進行氮營養(yǎng)診斷，雖已在小麥[9?11]、水稻[12?13]、馬鈴薯[14]等作物上得到大量應(yīng)用并驗證了其較高的精度。然而，臨界氮濃度的確認及后續(xù)數(shù)據(jù)處理需要大量計算，不利于實際生產(chǎn)應(yīng)用。此外，用Lemaire 等[4]方法構(gòu)建臨界氮稀釋曲線時，可能會由于氮稀釋效應(yīng)較低出現(xiàn)較多的離散點。這種數(shù)據(jù)分散性使得臨界氮稀釋曲線的擬合具有較大的偏差[24]。因此，在保證臨界氮稀釋曲線診斷能力的同時簡化復(fù)雜的步驟具有重要意義。Makowski 等[24]提出的臨界氮稀釋曲線即應(yīng)用貝葉斯統(tǒng)計構(gòu)建的臨界氮稀釋曲線，可以跳過區(qū)分生長受氮限制組和不受氮限制組的步驟。通過MCMC 方法[24]，可以很容易地從基礎(chǔ)數(shù)據(jù)和先驗知識中獲得參數(shù)的后驗分布，這種方法可以大大提高建模效率。在本研究結(jié)果中，雖然貝葉斯曲線計算獲得的臨界氮濃度與實際臨界氮濃度間的RMSE 和n-RMSE 較兩步法曲線略高，但其對受氮限制和非受氮限制的區(qū)分效果顯著好于兩步法曲線；同時，從圖4 中可以看出，兩步法曲線擬合臨界氮濃度與實際臨界氮濃度相差較大的部分為實際臨界氮濃度lt;2.6% 的區(qū)間，在該區(qū)間內(nèi)，貝葉斯曲線擬合臨界氮濃度效果較好，葉片干物重一般gt;2.1 t/hm2，對應(yīng)烤煙移栽70 天以后 （表2），此時已進入烤煙葉片成熟期[40]，該時期對烤煙品質(zhì)形成至關(guān)重要[41]。由臨界氮稀釋曲線計算得到的氮營養(yǎng)指數(shù)（NNI） 被廣泛用于作物氮素營養(yǎng)診斷，其可以直觀地量化作物氮狀況，以優(yōu)化氮肥管理[12]。本研究結(jié)果表明，貝葉斯曲線計算的NNI 比兩步法計算的NNI 略大，主要原因是兩步法曲線在貝葉斯曲線之上，根據(jù)方程 （2） 進行計算，貝葉斯曲線計算的NNI 應(yīng)大于兩步法，但兩者間差異僅為4%，且呈線性關(guān)系，方程決定系數(shù)R2 為 0.96，穩(wěn)定度較高，在NNI 驗證中取得了很好的結(jié)果，因此應(yīng)用貝葉斯統(tǒng)計構(gòu)建的烤煙葉片氮濃度稀釋曲線計算的NNI，可在一定程度上代替兩步法來進行氮營養(yǎng)狀態(tài)診斷。

4 結(jié)論

應(yīng)用貝葉斯統(tǒng)計法構(gòu)建的烤煙葉片臨界氮稀釋曲線方程為Nc=2.74×LDM?0.22，其不確定水平為0.16%～0.70%，不確定水平隨著干物重的增加先降低后增加。貝葉斯曲線在擬合臨界氮濃度的效果上，與兩步法曲線差異不大，但較兩步法曲線可以更好的區(qū)分生長受氮限制組與不受氮限制組，且計算的氮營養(yǎng)指數(shù)與兩步法曲線的線性變異決定系數(shù)R2 高達0.96，因此，貝葉斯曲線可更簡單準確地用于評價烤煙氮營養(yǎng)狀況。

參 考 文 獻：

[ 1 ]Zhang F S， Cui Z L， Chen X P， et al. Integrated nutrient management for food security and environmental quality in China[J]. Advances in Agronomy， 2012， 116： 1?32.

[ 2 ]榮凡番， 張忠鋒， 張教俠， 等. 不同施氮水平對山東褐土煙田供氮狀況及烤煙質(zhì)量的影響[J]. 中國煙草科學， 2009， 30（2）： 41?46.

Rong F F， Zhang Z F， Zhang J X， et al. Effects of nitrogen rates on brown soil nitrogen status of tobacco fields and tobacco quality in Shandong[J]. Chinese Tobacco Science， 2009， 30（2）： 41?46.

[ 3 ]Chen Y J， Ren K， He X， et al. The response of flue-cured tobacco cultivar K326 to nitrogen fertilizer rate in China[J]. The Journal of Agricultural Science.， 2020， 158（5）： 371?382.

[ 4 ]Lemaire G， Gastal F， Salette J. Analysis of the effect of N nutrition on dry matter yield of a sward by reference to potential yield and optimum N content[A]. Proceedings of the 16th international grassland congress[C]. Nice， France： The French Grassland Society，1989.

[ 5 ]Lemaire G， Sinclair T， Sadras V， Bélanger G. Allometric approach to crop nutrition and implications for crop diagnosis and phenotyping. A review[J]. Agronomy for Sustainable Development， 2019， 39（2）： 27.

[ 6 ]Justes E， Mary B， Meynard J M， et al. Determination of a critical nitrogen dilution curve for winter wheat crops[J]. Annals of Botany，1994， 74（4）： 397?407.

[ 7 ]Greenwood D J， Gosse L G， Cruz P， et al. Decline in percentage N of C3 and C4 crops with increasing plant mass[J]. Annals of Botany，1990， 66（4）： 425?436.

[ 8 ]Chen P， Zhu Y. A new method for winter wheat critical nitrogen curve determination[J]. Agronomy Journal， 2013， 105（6）： 1839?1846.

[ 9 ]趙犇， 姚霞， 田永超， 等. 基于臨界氮濃度的小麥地上部氮虧缺模型[J]. 應(yīng)用生態(tài)學報， 2012， 23（11）： 3141?3148.

Zhao B， Yao X， Tian Y C， et al. Accumulative nitrogen deficit models of wheat above ground part based on critical nitrogen concentration[J]. Chinese Journal of Applied Ecology， 2012， 23（11）：3141?3148.

[10]田興帥， 李松陽， 張羽， 等. 基于臨界氮濃度稀釋曲線的小麥氮肥需求量估測研究[J]. 麥類作物學報， 2019， 39（9）： 1112?1120.

Tian X S， Li S Y， Zhang Y， et al. Estimation of nitrogen fertilizer requirement of winter wheat based on critical nitrogen dilution curve[J]. Journal of Triticeae Crops， 2019， 39（9）： 1112?1120.

[11]馬曉晶， 張小濤， 黃玉芳， 葉優(yōu)良. 小麥葉片臨界氮濃度稀釋曲線的建立與應(yīng)用[J]. 植物生理學報， 2017， 53（7）： 1313?1321.

Ma X J， Zhang X T， Huang Y F， Ye Y L. Establishment and application of critical nitrogen dilution curve for wheat based on leaf dry matter[J]. Plant Physiology Journal， 2017， 53（7）： 1313?1321.

[12]Ataulkarim S T， Yao X， Liu X J， et al. Development of critical nitrogen dilution curve of japonica rice in Yangtze River Reaches[J].Field Crops Research， 2013， 149： 149?158.

[13]Zhang K， Liu X J， Ma Y， et al. A new canopy chlorophyll index based paddy rice critical nitrogen dilution curve in Eastern China[J].Field Crops Research， 2021， 266： 108139.

[14]Belanger G， Walsh J R， Richards J E， et al. Critical nitrogen curve and nitrogen nutrition index for potato in Eastern Canada[J].American Journal of Potato Research， 2001， 78（5）： 355?364.

[15]楊慧， 曹紅霞， 柳美玉， 劉世和. 水氮耦合條件下番茄臨界氮濃度模型的建立及氮素營養(yǎng)診斷[J]. 植物營養(yǎng)與肥料學報， 2015， 21（5）：1234?1242.

Yang H， Cao H X， Liu M Y， Liu S H. Simulation of critical nitrogenconcentration and nitrogen nutrition index of tomato under different water and nitrogen conditions[J]. Journal of Plant Nutrition and Fertilizers， 2015， 21（5）： 1234?1242.

[16]王新， 馬富裕， 刁明， 等. 滴灌番茄臨界氮濃度、氮素吸收和氮營養(yǎng)指數(shù)模擬[J]. 農(nóng)業(yè)工程學報， 2013， 29（18）： 99?108.

Wang X， Ma F Y， Diao M， et al. Simulation of critical nitrogen concentration， nitrogen uptake and nitrogen nutrition index of processing tomato with drip irrigation[J]. Transactions of the Chinese Society of Agricultural Engineering， 2013， 29（18）： 99?108.

[17]付江鵬， 賀正， 賈彪， 等. 滴灌玉米臨界氮稀釋曲線與氮素營養(yǎng)診斷研究[J]. 作物學報， 2020， 46（2）： 290?299.

Fu J P， He Z， Jia B， et al. Critical nitrogen dilution curve and nitrogen nutrition diagnosis of maize with drip irrigation[J]. Acta Agronomica Sinica， 2020， 46（2）： 290?299.

[18]李正鵬， 宋明丹， 馮浩. 關(guān)中地區(qū)玉米臨界氮濃度稀釋曲線的建立和驗證[J]. 農(nóng)業(yè)工程學報， 2015， 31（13）： 135?141.

Li Z P， Song M D， Feng H. Development and validation of critical nitrogen content curve for maize in Guanzhong area[J]. Transactions of the Chinese Society of Agricultural Engineering， 2015， 31（13）：135?141.

[19]Chen P F， Wang J H， Huang W J， et al. Critical nitrogen curve and remote detection of nitrogen nutrition index for corn in the Northwestern Plain of Shandong Province， China[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2013， 6（2）： 682?689.

[20]Colnenne C， Meynard J M， Reau R， et al. Determination of a criticalnitrogen dilution curve for winter oilseed rape[J]. Annals of Botany，1998， 81（2）： 311?317.

[21]Sieling K， Kage H. Organ-specific critical N dilution curves and derived NNI relationships for winter wheat， winter oilseed rape and maize[J]. European Journal of Agronomy， 2021， 130： 126365.

[22]Cheng M H， He J Y， Wang H D， et al. Establishing critical nitrogen dilution curves based on leaf area index and aboveground biomass for greenhouse cherry tomato： A Bayesian analysis[J]. European Journal of Agronomy， 2022， 141： 126615.

[23]Bohman B J， Culshaw-Maurer M J， Ben Abdallah F B， et al. Quantifying critical N dilution curves across G × E × M effects for potato using a partially-pooled Bayesian hierarchical method[J].European Journal of Agronomy， 2023， 144： 126744.

[24]Makowski D， Zhao B， Ata-Ul-Karim S T， et al. Analyzing uncertainty in critical nitrogen dilution curves[J]. European Journal of Agronomy，2020， 118： 126076.

[25]Yao B， Wang X L， Lemaire G， et al. Uncertainty analysis of critical nitrogen dilution curves for wheat[J]. European Journal of Agronomy，2021， 128： 126315.

[26]Ciampitti I A， Fernandez J， Tamagno S， et al. Does the critical N dilution curve for maize crop vary across genotype x environment x management scenarios？ A Bayesian analysis[J]. European Journal of Agronomy， 2021， 123： 126202.

[27]Fernández J A， Lemaire G， Bélanger G， et al. Revisiting the critical nitrogen dilution curve for tall fescue： A quantitative synthesis[J].European Journal of Agronomy， 2021， 131： 126380.

[28]Fernandez J A， Versendaal E， Lacasa J， et al. Dataset characteristics for the determination of critical nitrogen dilution curves： From past to new guidelines[J]. European Journal of Agronomy， 2022， 139：126568.

[29]解夢凡， 賈海江， 曲遠凱， 等. 種植密度和氮肥用量對百色煙區(qū)烤煙葉片發(fā)育及煙葉產(chǎn)量的影響[J]. 作物雜志， 2024， （2）： 189?197.

Xie M F， Jia H J， Qu Y K， et al. Effects of planting density and nitrogen fertilizer application rate on leaf development and yield of flue-cured tobacco in Baise tobacco region.[J]. Crops， 2024， （2）：189?197.

[30]何仲秋， 王曉琳， 張啟明， 等. 東南煙稻輪作區(qū)烤煙臨界氮濃度稀釋曲線的建立與驗證[J]. 植物營養(yǎng)與肥料學報， 2021， 27（11）：2001?2009.

He Z Q， Wang X L， Zhang Q M， et al. Construction and verification of a critical nitrogen dilution curve for flue-cured tobacco under tobacco-rice rotation systems in Southeast China[J]. Journal of PlantNutrition and Fertilizers， 2021， 27（11）： 2001?2009.

[31]周方， 席奇亮， 張思琦， 等. 豫中煙區(qū)水肥一體化條件下烤煙臨界氮濃度稀釋曲線與氮素營養(yǎng)診斷研究[J]. 中國煙草學報， 2019， 25（1）：58?66.

Zhou F， Xi Q L， Zhang S Q， et al. Critical nitrogen concentration dilution curve and nitrogen nutrition diagnosis of flue-cured tobacco under conditions of integration of water and fertilizer in tobaccogrowing areas in Central Henan[J]. Acta Tabacaria Sinica， 2019，25（1）： 58?66.

[32]Yao X， Ata-Ul-Karim S T， Zhu Y， et al. Development of critical nitrogen dilution curve in rice based on leaf dry matter[J]. European Journal of Agronomy， 2014， 55： 20?28.

[33]付江鵬， 賈彪， 楊文偉， 等. 基于葉片干物質(zhì)的滴灌玉米臨界氮稀釋曲線構(gòu)建[J]. 應(yīng)用生態(tài)學報， 2020， 31（3）： 945?952.

Fu J P， Jia B， Yang W W， et al. Development of critical nitrogen dilution curve based on leaf dry matter for maize under drip irrigation[J]. Chinese Journal of Applied Ecology， 2020， 31（3）：945?952.

[34]Zhao B， Ata-Ul-Karim S T， Liu Z D， et al. Development of a critical nitrogen dilution curve based on leaf dry matter for summer maize[J].Field Crops Research， 2017， 208： 60?68.

[35]馬興華， 梁曉芳， 劉光亮， 等. 氮肥用量及其基追施比例對烤煙氮素利用的影響[J]. 植物營養(yǎng)與肥料學報， 2016， 22（6）： 1655?1664.

Ma X H， Liang X F， Liu G L， et al. Effect of nitrogen application rate and base and topdressing ratio on nitrogen utilization of flue-cured tobacco[J]. Journal of Plant Nutrition and Fertilizers， 2016， 22（6）：1655?1664.

[36]劉力瑋， 劉洋， 王藝， 等. 山東煙區(qū)烤煙地上部臨界氮稀釋曲線的構(gòu)建與營養(yǎng)診斷[J]. 中國煙草科學， 2024， 45（3）： 19?25.

Liu L W， Liu Y， Wang Y， et al. Construction and nutritional diagnosis of shoot critical nitrogen dilution curve of flue-cured tobacco in Shandong tobacco growing region[J]. Chinese Tobacco Science， 2024， 45（3）： 19?25.

[37]王婧， 鮑貴. 貝葉斯統(tǒng)計與傳統(tǒng)統(tǒng)計方法的比較[J]. 統(tǒng)計與決策，2021， （1）： 24?29.

Wang J， Bao G. Comparison between Bayesian statistics and traditional statistical methods[J]. Statistics amp; Decision， 2021， （1）：24?29.

[38]Yin M， Li Y， Xu L， et. al. Nutrition diagnosis for nitrogen in winter wheat based on critical nitrogen dilution curves[J]. Crop Science. ，2018. 58： 416-425.

[39]Zhao Y， Chen P F， Li Z H， et al. A modified critical nitrogen dilution curve for winter wheat to diagnose nitrogen status under different nitrogen and irrigation rates[J]. Frontiers in Plant Science， 2020， 11：549636.

[40]孫延國， 馬興華， 黃擇祥， 等. 煙草溫光特性研究與利用： Ⅰ. 氣象因素對山東煙區(qū)主栽品種生育期的影響[J]. 中國煙草科學， 2020，41（1）： 30?37.

Sun Y G， Ma X H， Huang Z X， et al. Investigation and utilization of temperature and light characteristics of tobacco： I. effects of meteorological factors on growth period of main varieties in Shandong tobacco-growing areas[J]. Chinese Tobacco Science， 2020，41（1）： 30?37.

[41]李春儉， 張福鎖， 李文卿， 等. 我國烤煙生產(chǎn)中的氮素管理及其與煙葉品質(zhì)的關(guān)系[J]. 植物營養(yǎng)與肥料學報， 2007， 13（2）： 331?337.

Li C J， Zhang F S， Li W Q， et al. Nitrogen management and its relation to leaf quality in production of flue-cured tobacco in China[J]. Journal of Plant Nutrition and Fertilizers， 2007， 13（2）：331?337.