徐藝，邵光成*，丁鳴鳴，章二子，何菁

(1. 河海大學農(nóng)業(yè)工程與科學學院，江蘇 南京 210098； 2. 南京市水務局，江蘇 南京 210036； 3. 江寧淳化街道水務站,江蘇 南京 211122； 4. 南京市水利建筑工程檢測中心有限公司,江蘇 南京 210036)

表層土壤全氮含量(total nitrogen，TN)為氮含量的總和，包括有機氮、硝態(tài)氮、亞硝態(tài)氮以及銨態(tài)氮，還包括部分聯(lián)氮、偶氮和疊氮等含氮化合物.傳統(tǒng)測定全氮含量的凱氏定氮法雖然精度高，然而步驟復雜、成本高，不適用于田間實際應用[1]，而遙感監(jiān)測具有以下優(yōu)點：一方面，遙感不需要直接接觸觀測對象，從而避免對土壤表面作物的傷害.另一方面，可以開展長時間連續(xù)監(jiān)測，而且它的數(shù)據(jù)可以很容易地集成到地理信息系統(tǒng)中，以便后續(xù)進行分析[2].目前，常用的遙感技術包括衛(wèi)星遙感、無人機遙感、載人機遙感和地基遙感[3-4].相對于其他遙感平臺，無人機遙感具有成本低[5]、能獲取高空間分辨率影像、適合田塊尺度觀測等優(yōu)勢[6]，在農(nóng)情監(jiān)測中占據(jù)了重要位置[7].目前常采用的多光譜反演方法有偏最小二乘法、LASSO縮減系數(shù)法、主成分回歸等.其中，偏最小二乘回歸將關聯(lián)分析、主成分分析與多元線性回歸分析相結合[8]，使用來自響應變量的信息以及預測因子進行特征轉換[9]，曾用于熱帶森林的光譜和化學分析[10]，并作為處理大型高光譜數(shù)據(jù)集的常用方法[11]，但在采用最小二乘法求解回歸系數(shù)時發(fā)現(xiàn)光譜數(shù)據(jù)的自相關程度高，存在嚴重的多重共線性問題，反演結果不可靠[12].楊福芹等[13]基于多重共線性，篩選出對冬小麥氮營養(yǎng)指數(shù)相關性較高的圖像指數(shù)，再利用偏最小二乘法構建反演模型，建模集均方根誤差RMSE可達0.085 8，驗證集RMSE達0.187 1，預測精度較高；LAURIN等[14]采用偏最小二乘回歸(PLSR)，利用高光譜數(shù)據(jù)與植被指數(shù)，對非洲熱帶雨林生物量進行建模，發(fā)現(xiàn)改進后模型精度(決定系數(shù)R2=0.70)優(yōu)于不考慮多重共線性的回歸模型(R2=0.64).

高光譜數(shù)據(jù)通常由100多個帶寬為10 nm或更小的連續(xù)波段組成，波段之間相關性高，存在多重共線性和“維數(shù)災難”的問題，對反演模型的可靠性造成嚴重影響[3].為彌補高光譜數(shù)據(jù)中的多重共線性問題，改善預測模型精度，文中基于無人機多光譜高精度影像，獲取其光譜反射率數(shù)據(jù)，分析農(nóng)田土壤表層全氮含量實測值與光譜反射率的多元線性回歸中的多重共線性問題，構建基于嶺回歸的無人機遙感影像反演土壤全氮含量預測模型，旨在探索一種兼顧反演精度與光譜數(shù)據(jù)多重共線性問題的方法，以便為無人機遙感土壤氮素營養(yǎng)診斷提供理論依據(jù).

1 材料與方法

1.1 研究區(qū)概況

試驗于2021年4月29日在位于南京市江寧區(qū)淳化街道某小型灌區(qū)插秧前農(nóng)田進行.該農(nóng)田位于119°4′0″～119°4′6″E，31°54′11″～34°54′17″N，農(nóng)田總面積約為6 000 m2，無植被覆蓋，有利于觀測表層土壤光譜反射率.

1.2 數(shù)據(jù)采集與處理

1.2.1 數(shù)據(jù)采集

試驗采用棋盤式布點法，選定10塊條田，每塊再選取6個取樣點，每個點位在表土層深度0～30 cm進行取樣，并使用紫外分光度計法測定土壤樣本的全氮含量，測點分布如圖1所示.

圖1 采樣點分布圖Fig.1 Distribution of measurement sampling points

農(nóng)田遙感影像采用大疆P4 Mulitispectral無人機獲取.該型號無人機共搭載6個1/2.9 in CMOS影像傳感器，其中1個彩色傳感器用于常規(guī)可見光(RGB)成像，5個單色傳感器用于多光譜成像.單色傳感器前濾光片可通過波段：藍光波段(B)，450±16 nm；綠光波段(G)，560±16 nm；紅光波段(R)，650±16 nm；紅邊波段(RE)，730±16 nm；近紅外波段(NIR)，840±26 nm.航拍在4月29日正午11：00—12：00進行，飛行高度50 m，飛行速度5 m/s，規(guī)劃航點299個，航向重疊率90%，旁向重疊率75%，主航線角度262°，主航線9條，云臺俯仰角-90°，共拍攝圖片1 794張，采用ArcGIS 10.7對圖片進行拼接.該航高下遙感影像的地面分辨率為2.65 cm/pixel，滿足精度要求.

1.2.2 遙感影像拼接

精靈4無人機懸停拍攝的圖片為帶有坐標的TIF格式，將其導入ArcGIS后，指定坐標系為WGS1984，利用ArcToolbox中的“鑲嵌”工具進行拼接，最終得到5張不同波段的農(nóng)田影像.

1.2.3 光譜反射率計算

由無人機拍攝的原始影像數(shù)據(jù)點為像元值(digital number).為了得到對應波段的地表反射率，后續(xù)處理參照《P4 MultiSpectral 圖像處理指南》[15]進行.

以求解藍光波段反射率RB為例(按文獻[15]格式)為

(1)

式中：Bluecamera為圖像信號值；pCamBlue為相機參數(shù)，查遙感圖像EXIF中XMP-drone-dji項可得1.355 955；BlueLS為藍光光強傳感器信號值，pLSBlue為藍光傳感器校準參數(shù)，兩者乘積可直接查XMP-drone-dji下的Irradiance得到，BlueLS×pLSBlue=12 629.98.

Bluecamera可計算為

(2)

式中：IBlue和IBlackLevel分別為歸一化到值域[0,1]上的像素值和黑電平的值，其中IBlue計算方式參考式(3)，IBlackLevel可以在遙感圖像信息中的EXIF-IFD0-BlackLevel得到，文中IBlackLevel=4 096/65 635=0.062 4；Bluegain為相機曝光時的增益參數(shù)；Blueetime為曝光時間.

查圖像信息XMP-drone-dji與XMP-Camera項，可得Bluegain=1，Blueetime=0.364.

(3)

式中：DNBlue為藍光波段像元值.

經(jīng)過上述計算，可得到土壤樣本5個波段的光譜反射率R曲線，如圖2所示，圖中λ為波長.

圖2 光譜反射率曲線Fig.2 Spectral reflectance curve

2 結果與分析

2.1 多元線性回歸建模

假定以式(4)形式的回歸模型，反演表土全氮含量，即

TNtopsoil=β1RB+β2RG+β3RR+

β4RRE+β5RNIR+ε，

(4)

式中：TNtopsoil為表土全氮含量，g/kg；RB，RG，RR，RRE，RNIR分別為藍光波段、綠光波段、紅光波段、紅邊波段、近紅外波段的地表光譜反射率；β1—β5為回歸系數(shù)；ε為殘差.

以五波段反射率為自變量，土壤全氮含量為因變量，50組有效數(shù)據(jù)作為樣本，利用SPSS 21.0軟件，采用最小二乘回歸求解系數(shù)β1—β5，結果為

TNtopsoil=-8.437RB+18.310 4RG-28.724RR-
18.569 3RRE+33.995RNIR.

(5)

該模型回歸系數(shù)R2達到了0.504，然而在顯著性t檢驗中，回歸系數(shù)均大于0.05，說明多元線性回歸系數(shù)不具有統(tǒng)計學意義.

2.2 回歸系數(shù)不合理原因

一般多元線性回歸模型可概化為

y=Xβ+ε，

(6)

式中：X={x1,x2,…,xn}為解釋變量矩陣；y為被解釋變量矩陣；β為待估計的回歸系數(shù)；ε為殘差，滿足數(shù)學期望E(ε)=0.

依據(jù)最小二乘法，回歸系數(shù)的參數(shù)估計值可表示為

(7)

多元線性回歸理論假定，參與回歸的各自變量間線性無關.若該假設不滿足，則會導致模型對誤差ε極敏感，回歸系數(shù)不可靠.然而，當自變量之間存在嚴重多重共線性問題時，XTX的行列式值接近于0，導致回歸系數(shù)估計值的解非常不穩(wěn)定.

由實測反射率數(shù)據(jù)可得

(8)

使用最小二乘法計算回歸系數(shù)時，需要使用矩陣XTX的逆，而|XTX|=1.164 6×10-7接近于0，將導致計算出的回歸系數(shù)過大，且當解釋變量發(fā)生微小擾動時，回歸系數(shù)波動劇烈，甚至改變符號.為衡量自變量之間的多重共線性程度，MARQUARDT[16]于1970年提出方差膨脹因子(variance inflation factor, VIF)，計算公式為

(9)

式中：VIFi為對應自變量Xi的方差膨脹系數(shù)，i=1,2,…,6.

方差膨脹因子根據(jù)式(9)計算可得

(10)

由式(10)可知，與波段反射率RB，RG，RR，RRE，RNIR對應的方差膨脹因子分別為17.804 2，254.710 1，1 047.716 9，986.509 9，103.774 1.根據(jù)HOERL等[17]的研究，當VIF>10時，說明自變量之間存在多重共線性問題；若VIF>100，說明多重共線性現(xiàn)象嚴重.本次所測量的光譜數(shù)據(jù)中，有4個波段的VIF在100以上，說明光譜反射率數(shù)據(jù)存在較嚴重的多重共線性問題.

2.3 嶺回歸建模

嶺回歸是一種在自變量高度相關的情況下估計多元回歸模型系數(shù)的方法.該理論由HOERL等[18]首次提出.當線性回歸模型具有一些高度相關的獨立變量時，嶺回歸為最小二乘估計不精確的一種可靠的解決方案，它通過人為引入懲罰項kIp，回避了XTX行列式接近0的問題，即

(11)

k取值在0到1.當k=0時，嶺回歸退化為最小二乘估計.由于引入人為誤差項，該估計為有偏估計.HOERL等[18]證明，存在k值使得嶺回歸參數(shù)估計值的均方誤差小于最小二乘估計，并提出嶺跡法以確定合適的嶺回歸系數(shù)；嶺跡為所有標準化回歸系數(shù)與k的曲線圖，k∈[0,1].此外，HOERL等[18]提出了4個選擇最佳k值應滿足的條件：① 嶺跡線基本穩(wěn)定；② 回歸系數(shù)沒有不合理的數(shù)值大小；③ 回歸系數(shù)不再正負波動，符號變得合理；④ 殘差平方和相較多元線性回歸沒有顯著增加.

圖3為嶺跡圖及各方差膨脹因子隨k的變化圖，圖中VIF1，VIF2，VIF3，VIF4，VIF5分別表示回歸系數(shù)β1—β5對應的方差膨脹因子.其中圖3a為嶺跡圖，βs1—βs5為標準化回歸系數(shù)，是將自變量矩陣與因變量進行z-score標準化[19]后回歸得到的系數(shù)，分別與回歸系數(shù)β1—β5對應.由嶺跡圖可知：① 當k=0時，回歸系數(shù)βs2，βs3，βs5的絕對值較穩(wěn)定時偏大，說明βs2，βs3，βs5被嚴重高估；隨著k值增大，βs1—βs5的絕對值都逐漸減小，且趨于穩(wěn)定.② 隨著k值增大，βs2由正轉負，最后達到穩(wěn)定.③ 在k≥0.025之后，標準化回歸系數(shù)沒有不合理的值.

由圖3b—3f可知，5個方差膨脹因子均隨k增加而迅速減小，當k=0.04時，VIF均小于10，此時可認為多重共線性對預測模型的精度影響較小.

圖4為R2、均方根誤差與P值隨k值的變化情況.嶺回歸與其他加入人為懲罰項的回歸方法類似，模型回歸擬合精度會因懲罰因子的增大而迅速降低，因此k不宜設置過大.圖4c中P1—P5分別對應回歸系數(shù)β1—β5的P值，用于評估不同嶺回歸系數(shù)下的顯著性水平[20].當P≤0.05時，可認為通過回歸系數(shù)的顯著性檢驗；當k≥0.12時，P1—P5均小于0.05，此時所有的回歸系數(shù)都通過顯著性檢驗.

綜上所述，當k=0.12時，回歸決定系數(shù)R2從0.504降至0.408，且此時回歸系數(shù)趨于穩(wěn)定，數(shù)值合理，各波段方差膨脹因子均小于10，且P值小于0.05，表明所有回歸系數(shù)均通過顯著性檢驗.因此，取k=0.12作為最佳嶺回歸參數(shù)，此時表土全氮含量估算值可由五波段反射率表示為

TNtopsoil=1.931 6RB-5.308 4RG-10.538 9RR-
0.324 4RRE+13.746 1RNIR.

(12)

圖3 嶺跡圖和方差膨脹因子隨k值的變化圖Fig.3 Ridge trace and variation of variance inflation factors with k

圖4 R2、均方根誤差與P值隨k值的變化Fig.4 Variation of R2, root mean square error and variation of P with k

2.4 預測模型評價

圖5為擬合結果，圖中TNp，TNt分別為土全氮的預測值、實測值.圖5a—5b是多元線性回歸與嶺回歸的反演結果.

圖5 擬合結果Fig.5 Fitting results

多元線性回歸反演全氮含量的R2為0.504，均方根誤差RMSE為0.472，均優(yōu)于k=0.12時嶺回歸預測值；然而由于多重共線性的影響，線性回歸在模型預測時表現(xiàn)較差.圖5c—5d是利用反演模型對驗證集數(shù)據(jù)進行預測的結果，其中驗證集數(shù)據(jù)由10個條田中各自隨機選取1個樣本點構成.可以看出，多元線性回歸預測模型的R2為0.645，RMSE為0.820，對部分點位全氮含量嚴重高估，建模效果雖好，但驗證效果差.而嶺回歸預測模型RMSE下降不明顯，有較好的預測效果.

通過式(12)并結合農(nóng)田光譜反射率數(shù)據(jù)圖，可以得到如圖6所示的完整表土全氮含量反演結果.利用表土全氮含量分布圖，可以快速獲取土壤全氮含量信息，對無人機遙感土壤氮素營養(yǎng)診斷和精準施肥具有重要意義.

圖6 農(nóng)田表土全氮含量反演結果Fig.6 Inversion results of total nitrogen content in farmland topsoil

3 討 論

反演模型中，近紅外波段反射率與表土全氮含量的相關度最高.這是因為近紅外波段光譜的信息來源于分子振動的倍頻和合頻，常用于含C-H，N-H，O-H等基團的有機物分析[21]，而土壤中氮素絕大部分為有機結合態(tài)，與有機質(zhì)關聯(lián)密切[22]，因此高分辨率的近紅外波段包含了土壤全氮含量的敏感譜區(qū).利用近紅外全譜波段對土壤全氮含量反演可以達到較高的精度.如李頡[23]利用12 500～3 600 cm-1的近紅外光譜數(shù)據(jù)建立的土壤全氮含量偏最小二乘回歸模型，決定系數(shù)可達到89.63%.利用多光譜數(shù)據(jù)進行全氮含量反演時，由于土壤全氮含量的取值區(qū)間較狹窄，方差較小，不利于穩(wěn)定預測模型的建立，但通過相關分析，篩選與全氮含量相關性強、顯著性高的敏感波段進行建模，仍能達到較高精度[24].

文中研究嘗試了利用可見光波段(RGB)與紅邊(RE)、近紅外(NIR)波段光譜反射率數(shù)據(jù)，在不同嶺回歸系數(shù)下對土壤全氮含量進行預測，其中多元線性回歸模型是k=0時的特例.由于上述五波段數(shù)據(jù)獲取便利，雖然模型精度有所下降，但有助于大面積田塊全氮含量的快速、定性診斷.綜合反演精度與回歸系數(shù)顯著性兩方面考慮，選取k=0.12時的嶺回歸模型作為反演結果.

4 結 論

基于低空無人機搭載多光譜傳感器，通過獲取表土全氮含量和光譜反射率，研究分析了多元線性回歸在表土全氮含量光譜反演問題上的不足，揭示了光譜反射率數(shù)據(jù)特有的多重共線性問題及其應用限制.通過嶺回歸方法，以損失一定回歸精度為代價，得到一組對多重共線性不敏感的且穩(wěn)定的回歸系數(shù).

研究結果表明，通過表土全氮含量與光譜反射率多元線性回歸，發(fā)現(xiàn)波段RG，RR，RRE，RNIR對應的方差膨脹因子均大于100，自變量之間存在嚴重的多重共線性問題.而基于嶺回歸建模，當k=0.12時，回歸R2雖然有所降低，從0.504降至0.408，然而回歸系數(shù)趨于穩(wěn)定，數(shù)值合理，方差膨脹因子均小于10，且回歸系數(shù)之間差異具有統(tǒng)計學意義，說明基于光譜反射率信息反演土壤全氮含量的效果較好.但在實踐中獲取無人機影像進行應用時，還需要綜合考慮無人機不同飛行高度而引起的不同分辨率對反演結果的影響.