基于嶺回歸的土壤全氮含量反演模型

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

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

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

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

1 材料與方法

1.1 研究區(qū)概況

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

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

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

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

圖1 采樣點(diǎn)分布圖Fig.1 Distribution of measurement sampling points

農(nóng)田遙感影像采用大疆P4 Mulitispectral無(wú)人機(jī)獲取.該型號(hào)無(wú)人機(jī)共搭載6個(gè)1/2.9 in CMOS影像傳感器，其中1個(gè)彩色傳感器用于常規(guī)可見(jiàn)光(RGB)成像，5個(gè)單色傳感器用于多光譜成像.單色傳感器前濾光片可通過(guò)波段：藍(lán)光波段(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進(jìn)行，飛行高度50 m，飛行速度5 m/s，規(guī)劃航點(diǎn)299個(gè)，航向重疊率90%，旁向重疊率75%，主航線(xiàn)角度262°，主航線(xiàn)9條，云臺(tái)俯仰角-90°，共拍攝圖片1 794張，采用ArcGIS 10.7對(duì)圖片進(jìn)行拼接.該航高下遙感影像的地面分辨率為2.65 cm/pixel，滿(mǎn)足精度要求.

1.2.2 遙感影像拼接

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

1.2.3 光譜反射率計(jì)算

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

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

(1)

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

Bluecamera可計(jì)算為

(2)

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

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

(3)

式中：DNBlue為藍(lán)光波段像元值.

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

圖2 光譜反射率曲線(xiàn)Fig.2 Spectral reflectance curve

2 結(jié)果與分析

2.1 多元線(xiàn)性回歸建模

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

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

β4RRE+β5RNIR+ε，

(4)

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

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

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

(5)

該模型回歸系數(shù)R2達(dá)到了0.504，然而在顯著性t檢驗(yàn)中，回歸系數(shù)均大于0.05，說(shuō)明多元線(xiàn)性回歸系數(shù)不具有統(tǒng)計(jì)學(xué)意義.

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

一般多元線(xiàn)性回歸模型可概化為

y=Xβ+ε，

(6)

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

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

(7)

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

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

(8)

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

(9)

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

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

(10)

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

2.3 嶺回歸建模

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

(11)

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

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

由圖3b—3f可知，5個(gè)方差膨脹因子均隨k增加而迅速減小，當(dāng)k=0.04時(shí)，VIF均小于10，此時(shí)可認(rèn)為多重共線(xiàn)性對(duì)預(yù)測(cè)模型的精度影響較小.

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

綜上所述，當(dāng)k=0.12時(shí)，回歸決定系數(shù)R2從0.504降至0.408，且此時(shí)回歸系數(shù)趨于穩(wěn)定，數(shù)值合理，各波段方差膨脹因子均小于10，且P值小于0.05，表明所有回歸系數(shù)均通過(guò)顯著性檢驗(yàn).因此，取k=0.12作為最佳嶺回歸參數(shù)，此時(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 預(yù)測(cè)模型評(píng)價(jià)

圖5為擬合結(jié)果，圖中TNp，TNt分別為土全氮的預(yù)測(cè)值、實(shí)測(cè)值.圖5a—5b是多元線(xiàn)性回歸與嶺回歸的反演結(jié)果.

圖5 擬合結(jié)果Fig.5 Fitting results

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

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

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

3 討 論

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

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

4 結(jié) 論

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

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