張 漫 苗艷龍 仇瑞承 季宇寒 李 寒 李民贊

(1.中國農業(yè)大學現代精細農業(yè)系統集成研究教育部重點實驗室， 北京 100083;2.中國農業(yè)大學農業(yè)農村部農業(yè)信息獲取技術重點實驗室， 北京 100083)

0 引言

葉面積指數(Leaf area index, LAI)是重要的作物表型參數，它與光合作用、水循環(huán)、碳循環(huán)密切相關，是描述作物生長狀況和產量的重要指標[1-3]，因此，葉面積指數的快速測量具有重要意義。LAI測量方法分為直接測量和間接測量[4]，直接測量是對植株葉片面積進行測量，具有破壞性，消耗大量人力，且效率低下。間接測量是通過測量冠層輻射的透過率計算LAI，使用的傳感器有彩色相機、多光譜傳感器和激光雷達。

彩色圖像分析技術為快速測量LAI提供了一種快速、高效的方法[5-8]，但其易受外界環(huán)境影響，尤其是光照影響，測量精度還可以進一步提高。使用多光譜或高光譜數據可以實現作物LAI的準確、快速測量[9-14]，不過同樣存在受外界光照的影響。此外還有一些學者使用紅外圖像或衛(wèi)星遙感數據進行作物葉面積指數測量研究[15-20]。

激光雷達受外界光照影響較小，具有分辨率高、掃描速度快等優(yōu)點[21]，可以為LAI快速測量提供一種方法。文獻[22-23]使用機載激光雷達數據建立回波點云密度、激光穿透力指數、孔隙度和葉傾角與LAI的回歸模型，使用地面激光雷達數據，提出運用體素內葉片及其投影數學關系的真實葉面積獲取方法。文獻[24]使用地面激光雷達數據，運用球極平面投影和Lambert方位角等面積投影方法，計算點云數據的孔隙度，進而計算單木真實LAI。文獻[25]使用地面雷達數據，運用點云切片算法，計算森林的角度孔隙度和消光系數，進而估算森林的LAI。文獻[26]使用多個不同方向的車載二維激光雷達，使用迭代最近點算法配準點云數據，對室內玉米植株進行三維重建。地面雷達存在數據量大、算法復雜、測量范圍較小等問題；機載激光雷達分辨率較小，LAI計算精度有待提高；車載三維激光雷達可以有效彌補以上缺點。

本文以大喇叭口期的玉米為研究對象，使用車載16層三維激光雷達采集點云數據，提出一種以玉米分層點云數量或分層點云數量與地面點云數量比值為自變量計算玉米LAI的方法。

1 材料與方法

1.1 試驗材料與設備

選用Velodyne LiDAR公司的VLP-16型三維激光雷達采集玉米生長期的三維點云數據。VLP-16型三維激光雷達水平視角范圍0°～360°，分辨率0.1°～0.4°；垂直視角范圍0°～30°，分辨率2°；測量范圍為0.5～100 m；采集頻率為10 Hz。在Windows 7操作系統下，以Visual Studio 2013為平臺，安裝點云庫PCL
1.8.0(Point cloud library)和Cmake 3.8.0，使用C++語言進行軟件編程實現對三維點云數據的處理。葉面積指數真值獲取的方法有網格法、復印稱量法、回歸方程法等破壞性測量方法和LI葉面積儀、LAI冠層分析儀等非破壞測量方法。結合試驗需求，使用LI-COR公司的LAI-2200型植物冠層分析儀測量葉面積指數作為真實值。LAI-2200型植物冠層分析儀通過測定透過冠層后輻射能量的削弱程度確定葉面積指數；利用“魚眼”光學傳感器，垂直視野范圍0°～148°，水平視野范圍0°～360°，測量冠層上下5個角度光線強度的變化。

圖1 玉米種植及LAI-2200型植物冠層分析儀測量位置圖Fig.1 Location map of maize planting and LAI-2200 plant canopy analyzer

試驗地點為北京市中國農業(yè)大學上莊實驗站。試驗時間為2017年7月13日09:00—18:00，玉米品種為農大84和京農科728，其中：農大84為舒展型品種，京農科728為緊湊型品種。玉米種植區(qū)長64.8 m，寬11 m，共分為21個種植小區(qū)，每個小區(qū)種植3行，1個小區(qū)種植1個玉米品種，2種玉米交替種植，如圖1a所示。小區(qū)種植行間距為90 cm，株間距為30 cm，每行種植36株玉米，如圖1b所示，生長期為大喇叭口期。使用冠層分析儀測量玉米的葉面積指數。每行以0.9 m等間距確定6個測量樣本區(qū)域，6個測量區(qū)域中心點把該行均勻分為7個區(qū)域，2個測量中心點的距離為1.5 m，每個測量樣本區(qū)域長度為0.6 m，測量區(qū)域中心位置為玉米種植點，前后邊緣為相鄰玉米種植點，樣本測量區(qū)域包含測量中心點的1株完整玉米，和前后各1株玉米的部分。每個樣本區(qū)域測量3次，取平均值作為該區(qū)域的葉面積指數真值，共計獲得12個區(qū)域的葉面積指數，測量位置如圖1b所示。VLP-16型三維激光雷達安裝在車載測量平臺頂端，如圖2所示，采集平臺長為2.0 m，寬為1.3 m，垂直采集玉米三維點云數據，其與玉米冠層頂端的距離約為60 cm。采集玉米點云數據之前，為了降低地面起伏和松軟造成采集平臺俯仰傾斜和通行困難，對地面進行平整處理，保證激光雷達始終垂直采集。車載測量平臺以1 m/s左右的速度向前行進，VeloView軟件同時進行三維點云數據的采集，采集頻率為10 Hz，存儲為.pcap格式，用于后續(xù)處理。本次試驗采集了2行玉米的點云數據，每個玉米品種1行；其中：京農科728有32株，農大84有28株，總共獲得3 000幀點云數據。測量LAI真值區(qū)域共有400幀點云數據。

圖2 玉米三維點云數據采集Fig.2 3D point cloud data acquisition map of maize1.VLP-16型激光雷達 2.數據采集終端

1.2 三維點云數據處理

首先通過預處理去除背景點云，主要有其他區(qū)域玉米點云和測量平臺車體點云。預處理之后的點云數據包括玉米植株點云和地面點云；使用隨機采樣一致性平面分割算法分離玉米植株點云和地面點云；對玉米植株點云使用直通濾波算法進行分層處理，并統計每層點云個數。對得到的每層點云數量或其與地面點云數量比值和測量得到的葉面積指數，進行數據分析和模型建立。

1.2.1數據預處理

使用PCL軟件把三維點云數據，從.pcap格式轉換為.pcd格式。在每幀三維點云數據中，非測量區(qū)域玉米植株和測量平臺為背景點云，對后續(xù)處理無用，且影響處理時間。因此，本文采用簡單快速的直通濾波算法去除非測量區(qū)域玉米植株和測量平臺點云數據，獲得葉面積指數測量區(qū)域的三維點云數據。三維點云數據X軸與作物行垂直，Y軸與地面垂直，Z軸與作物行平行。田間采集試驗時，測量VLP-16型激光雷達的安裝位置，通過位置數據確定測量玉米區(qū)域的二維空間范圍，進而分別設置點云坐標X、Z的閾值范圍。

1.2.2點云分割

預處理后的點云由玉米植株點云和地面點云組成。因為地面起伏較小，地面點云組成一個近似平面的形狀，符合平面模型。選擇基于隨機采樣一致性(Random sample consensus, RANSAC)平面分割算法對點云數據進行分割。點云變化如圖3所示，平面分割步驟如下：

(1)讀入預處理之后的點云，對點云進行一次直通濾波，把點云分為兩部分：一部分為地面點云與少量植株點，另一部分為大量植株點。因為在植株稠密處點云數據，地面點數量占點云數量比例較低，難以分割出地面點云，分離出的數據多是幾個位于同一平面葉子的點云，所以要先進行一次直通濾波把點云數據分為上下兩層。

(2)對地面點與少量植株點數據進行基于隨機采樣一致性平面分割。一部分為分割平面點，即地面點云，另一部分為非平面點云，即植株點和凹陷地面點。

(3)對非平面點云Y軸數據與地面點平均值比較，大于平均值，即為地面點插入平面點云，得到最終的地面點云；并從非平面點云中刪除。

(4)把非平面點云插入到步驟(1)得到的植株點云中，獲得最終的植株點云。

圖3 平面分割點云變化Fig.3 Change of plane segmentation point cloud

1.2.3點云分層

玉米植株垂直結構分布一般為3層，植株上層、中層和下層，其中植株中層也稱為高效葉層[27]。玉米灌漿期植株達到最高，平均高度在300 cm，上層為240 cm以上，中層為120～240 cm，下層為120 cm以下。本次試驗采集玉米生長期為大喇叭口期，玉米平均高度沒有達到300 cm，需要對每層邊界高度進行等比例變化。

根據分層邊界高度使用直通濾波對植株點云進行分層，首先分離下層點云數據，對中層和上層點云進行兩種處理：①濾波處理，去除葉片邊緣噪聲點，對點云數據進行中層和上層分離。②不做處理，對點云數據進行中層和上層分離，記錄每層點云的數量，作為建立葉面積指數模型的自變量。通過對濾波處理和不做濾波處理建立的模型進行對比分析，判斷濾波處理的必要性。

1.2.4葉面積指數模型

模型的自變量分為兩類：一類是以玉米植株上、中、下每層點云數量H、M、L，或它們的組合作為自變量；另一類是以玉米上、中、下每層的點云數量與地面點云數量的比值Hr、Mr、Lr或它們的組合作為自變量。因變量為植物冠層分析儀測得的玉米葉面積指數LAI。對樣本進行數據分析，包括共線性檢測、線性回歸分析和驗證集檢驗。

對自變量進行共線性檢測，使用自變量之間的相關系數和方差膨脹因子(VIF)。如果VIF值小于10，則認為自變量之間不存在共線性。方差膨脹因子計算公式為

式中r——相關系數

DVIF——方差膨脹因子

對樣本數據進行一元線性回歸和多元線性回歸分析，建立玉米大喇叭口期的葉面積指數測量模型。使用決定系數R2、均方根誤差(RMSE)和相對均方根誤差(rRMSE)作為評價模型精度的指標。

一般R2值較大、RMSE和rRMSE值較小表示回歸模型測量的值準確度高。使用F檢驗和t檢驗對自變量和因變量之間的線性顯著性進行檢驗。

2 結果分析

2.1 數據預處理

測量VLP-16型激光雷達安裝位置與車架邊緣的距離，LAI-2200型植物冠層分析儀測量位置與車架邊緣的距離。確定點云坐標X和Z的閾值分別是X∈[-0.55 m, 0.45 m]、Z∈[-0.50 m, 0.50 m]，坐標Y不作處理。對玉米三維點云數據進行直通濾波處理，其中京農科728有208幀數據，農大84有216幀數據。直通濾波前有非測量區(qū)域玉米植株和測量平臺的點云數據，濾波后去除了這些無關數據，獲得測量區(qū)域的玉米植株點云數據和地面點云內數據，平均每幀點云數據有2 288個點，其中植株點2 008個，地面點280個。如圖4所示。

圖4 玉米點云直通濾波處理Fig.4 Maize point cloud pass filter processing

2.2 玉米植株點云和地面點云分割

讀入預處理后的點云數據，進行平面分割算法步驟(1)操作，設置坐標軸Y直通濾波的閾值為[2.7 m, 5.0 m]，獲得Y坐標大于2.7 m的地面點云數據和少量植株數據，Y坐標小于2.7 m的大量玉米點云數據。對地面點數據和少量植株數據進行步驟(2)操作，設置分割模型為平面模型，分割方法為隨機采樣一致性算法；根據地面起伏程度設置距離閾值為0.06 m，獲得分割平面點云數據和非平面點云數據。進行步驟(3)操作，京農科728和農大84地面點云數據Y坐標的平均值分別為2.95 m和2.96 m，非平面點云數據中Y坐標大于地面平均值的點，插入平面點云數據，得到所有地面點云，記錄地面點云數量；并在非平面點云中刪除這些點。進行步驟(4)操作，合并步驟(1)獲得的大量植株點云數據和步驟(3)獲得的非平面點云數據，得到完整植株點云數據。

2.3 玉米植株點分層

經過測量京農科728和農大84玉米植株的平均高度分別為250 cm和220 cm。依據玉米植株垂直結構分層方法，確定京農科728和農大84的分層邊界高度：京農科728上層為200 cm以上，中層為100～200 cm，下層為100 cm以下；農大84上層為176 cm以上，中層為88～176 cm，下層為88 cm以下。

讀入玉米植株點云數據，根據分層邊界高度，使用直通濾波算法對玉米植株點云數據進行分層，并記錄每層點云數量。對中層和上層進行基于統計分析濾波，臨近點數量設置為6，標準差倍數設置為1.2。對濾波后的點云進行分層處理，記錄每層點云數量。

2.4 玉米葉面積指數模型建立

分別將京農科728和農大84的200多組數據分為訓練集(Train set)和驗證集(Validation set)。隨機選擇其中的40組數據作為驗證集，其余數據作為訓練集。對每類處理自變量進行共線性檢測，計算自變量之間的相關系數和方差膨脹因子。對訓練集數據進行一元線性回歸和多元線性回歸分析，建立玉米大喇叭口期的葉面積指數測量模型，計算訓練集和驗證集的R2、RMSE和rRMSE，作為模型評價指標。

2.4.1共線性檢測

對玉米植株上、中、下3層點云數量H、M、L和上、中、下每層的點云數量與地面點云數量的比值Hr、Mr、Lr進行共線性檢測。針對不同的處理類型，京農科728中H、M、L或Hr、Mr、Lr間的相關系數和VIF如表1所示，其中H、M、L間的VIF的取值范圍為1.040～1.081，Hr、Mr、Lr間的VIF的取值范圍為1.197～1.414?？梢钥闯觯┺r科728引入地面點后，VIF值有小幅增加。農大84中H、M、L或Hr、Mr、Lr間的相關系數和VIF如表2所示，其中H、M、L

表1 京農科728中H、M、L或Hr、Mr、Lr間的相關系數和VIFTab.1 Correlation coefficients and VIF values among H, M, L or Hr, Mr, Lr of Jingnongke 728

表2 農大84中H、M、L或Hr、Mr、Lr間的相關系數和VIFTab.2 Correlation coefficients and VIF values among H, M, L or Hr, Mr, Lr of Nongda 84

間的VIF的取值范圍為1.093～1.994，Hr、Mr、Lr間的VIF的取值范圍為1.001～3.815?？梢钥闯鲛r大84引入地面點后有一定的增加，增幅大于京農科728。但是自變量之間的VIF都小于10，表明它們之間的共線性可以忽略，都可以用來建立線性回歸模型。

2.4.2模型及評價

分別以H、M、L為自變量，或者幾個變量的組合為自變量。以測量的葉面積指數值為因變量。使用訓練集數據建立一元線性回歸或多元線性回歸模型，計算訓練集和驗證集模型評價指標R2、RMSE和rRMSE。京農科728、農大84中H、M、L葉面積指數模型及評價如表3、4所示。對Hr、Mr、Lr自變量進行相同的處理，京農科728、農大84中Hr、Mr、Lr葉面積指數模型及評價如表5、6所示。

上述葉面積指數測量模型沒有通過F檢驗和t檢驗的有：表3中編號為3、11的模型，表5中編號為1、9的模型。模型通過F檢驗，沒有通過t檢驗的有：表3中編號為6、14的模型，表4中編號為5、13的模型，表5中編號為6、14的模型，表6中編號為

表3 京農科728中H、M、L葉面積指數模型及評價Tab.3 LAI model and evaluation indicators of H, M and L variables in Jingnongke 728

表4 農大84中H、M、L葉面積指數模型及評價Tab.4 LAI model and evaluation indicators of H, M and L variables in Nongda 84

表5 京農科728中Hr、Mr、Lr葉面積指數模型及評價Tab.5 LAI model and evaluation indicators of Hr, Mr and Lr variables in Jingnongke 728

表6 農大84中Hr、Mr、Lr葉面積指數模型及評價Tab.6 LAI model and evaluation indicators of Hr, Mr and Lr variables in Nongda 84

5、7、13、15的模型；其他模型均通過了F檢驗和t檢驗。從表中可以看出，沒有通過檢驗的模型，基本評價指標顯示的回歸效果不好，只有表6中的編號為7和15模型評價Lr變量沒有通過t檢驗，可以把Lr變量從模型中去除。

對通過F檢驗和t檢驗的模型進行對比分析可以得出，在每個表格中最優(yōu)三元線性回歸模型回歸效果好于最優(yōu)二元線性回歸模型，最優(yōu)二元線性回歸模型回歸效果好于最優(yōu)一元線性回歸模型。其中R2最大、RMSE和rRMSE最小的模型是三元線性回歸模型中編號為7或15的回歸模型。對比每個表格中無濾波處理和濾波處理后相應模型，R2、RMSE和rRMSE只有微小變化，濾波效果不顯著，說明在使用分層點云數量建立葉面積指數測量模型的方法中，濾波處理作用不大。原因是離群點數量占點云數量的比例很小，所以對于點云數量建立的模型影響?。粸V波會對中、上層點云數量產生一定比例的降低，可以通過系數的變化進行補償。因此，以下只對無濾波處理進行分析。

對表3～6進行分析，在一元線性回歸模型中，京農科728的M和Mr變量，農大84的H和Hr變量線性擬合度高于其他一元線性回歸模型。產生差異的原因是，京農科728是緊湊型玉米，葉傾角大，葉片緊湊，三維點云主要分布在中層，農大84是舒展型玉米，葉傾角小，葉片舒展，三維點云主要分布在上層。對比京農科728的M和Mr變量模型與農大84的H和Hr變量模型，引入地面點后，京農科728編號2的模型R2從0.660提升至0.759，農大84編號1的模型R2從0.839降至0.795，RMSE和rRMSE變化很小，所以引入地面點對一元線性回歸模型貢獻不大。

在二元線性回歸模型中，京農科728的H、M和Hr、Mr兩個變量組合，農大84的H、M和Hr、Mr兩個變量組合線性擬合度高于其他二元線性回歸模型。兩個品種的二元線性回歸模型中最優(yōu)模型自變量組合相同，說明在不同品種的葉面積指數測量模型中貢獻最大的層雖然不同，但是最優(yōu)變量組合卻是相同的。其原因是中、上層點云的數量占全部點云數量的90%以上，基本包含了點云數量的全部信息。并且H、M和Hr、Mr的二元線性回歸模型與其對應最優(yōu)一元線性回歸模型比較，R2有0.2左右的提升，RMSE和rRMSE降低了一半左右，明顯優(yōu)于一元線性回歸。對比表3和表5中編號4的模型，訓練集R2從0.812升至0.931，RMSE從0.109%降至0.066%，rRMSE從0.067%降至0.040%，驗證集R2從0.884升至0.949，RMSE從0.097%降至0.065%，rRMSE從0.060%降至0.042%；對比表4和表6中編號4的模型，訓練集R2從0.934升至0.979，RMSE從0.085%降至0.048%，rRMSE從0.035%降至0.020%，訓練集R2從0.926升至0.984，RMSE從0.087%降至0.040%，rRMSE從0.037%降至0.016%。可以看出引入地面點對不同玉米品種二元回歸模型的線性擬合度和測量精度都有很大程度的提高。所以選擇二元線性回歸模型中的Hr、Mr的變量組合，放棄H、M變量組合。

在三元線性回歸模型中，全部植株點云數量作為自變量的模型8回歸的效果很差，說明葉面積指數測量模型需要對點云植株進行分層處理。編號為7的模型在表3、4、5中通過了F檢驗和t檢驗，所以只對這3個模型進行分析。H、M、L和Hr、Mr、Lr三元線性回歸模型與其對應最優(yōu)二元線性回歸模型比較，R2只有0.04左右提高，RMSE和rRMSE只降低了20%左右。模型的線性擬合度和測量精度提高不大，且增加了變量，使模型變得復雜，且適應性降低。

綜合分析，選擇最優(yōu)二元線性回歸模型，放棄三元線性回歸模型。

3 討論

使用玉米的垂直結構分層，及每層點云的數量或每層點云數量與地面點云數量的比值，建立了LAI測量模型，為LAI田間快速測量提供了一種新的思路。本研究目前只建立了玉米一個生長期的LAI測量模型，后續(xù)還需要對玉米不同生長期的LAI測量模型進行研究，驗證該方法的通用性。

本文主要針對大喇叭口期的玉米開展研究，后續(xù)研究中，需進一步針對玉米的不同生長階段，采取不同的點云分析處理方法。例如在玉米苗期、拔節(jié)期，玉米植株低矮，葉片較小，就需要先考慮玉米植株是否同樣是存在3層的垂直分布，并根據植株的垂直分布對植株點云進行分層處理；在抽雄期、吐絲期等生育后期，植株完全長成，葉片交錯縱橫，遮擋嚴重，造成地面點云數量減少，可能會導致地面點云與植株點云分離出現問題。在模型建立部分，由于不同生育時期的植株點云和地面點云數量都存在不同程度的變動，在苗期時地面點云數量多，在吐絲期時地面點云數量少，依然要考慮是否要引入地面點云數量的問題。

4 結論

(1)提出了一種基于車載三維激光雷達的葉面積指數測量算法，通過對玉米植株點云和地面點云進行分離，并根據垂直結構對玉米植株進行分層處理，建立了較高精度的葉面積指數計算模型。

(2)通過對點云濾波處理前、后建立的LAI測量模型進行分析比較可知，R2、RMSE和rRMSE沒有顯著優(yōu)化。即有無濾波處理對基于分層點云數量建立LAI測量模型的影響可忽略不計。

(3)通過對建立的葉面積指數測量模型進行綜合分析可知，京農科728和農大84雖然株型不同，但都選擇了以Hr、Mr為變量建立的葉面積指數測量模型，說明此方法在玉米的大喇叭口期具有一定的通用性。