基于無人機光譜遙感的田塊尺度蒸散發(fā)空間分布估算

蘇寶峰 王 琮 張茹飛 陳 山

(1.西北農(nóng)林科技大學機械與電子工程學院，陜西楊凌 712100；2.農(nóng)業(yè)農(nóng)村部農(nóng)業(yè)物聯(lián)網(wǎng)重點實驗室，陜西楊凌 712100)

0 引言

蒸散是農(nóng)業(yè)水資源管理中灌溉管理和產(chǎn)量預報的重要參考指標[1-3]。準確的區(qū)域蒸散信息可以為區(qū)域水資源的管理和農(nóng)田作物培養(yǎng)提供信息以及決策支持[4-6]。然而，目前缺少對農(nóng)田蒸散量空間模式的有效量化方法，獲取精準的農(nóng)田蒸散量分布圖對于開展精準農(nóng)業(yè)研究、灌溉作業(yè)非常必要[4]。

目前，獲取農(nóng)田蒸散量的方法主要是渦動協(xié)方差法和衛(wèi)星遙感估算法[7-8]。渦動協(xié)方差法結果較為準確，但數(shù)據(jù)不具備空間尺度差異性[8]；使用衛(wèi)星數(shù)據(jù)，估算結果精度低、分辨率低與實時性差是制約衛(wèi)星技術在農(nóng)業(yè)中應用的重要因素[9-11]。無人機技術監(jiān)測土壤和植被信息具有快速、宏觀、方便等特點，已被廣泛應用于農(nóng)業(yè)領域[11-12]。因此，基于無人機技術可以建立適用于農(nóng)田的蒸散模型，然而合適的模型和數(shù)據(jù)的獲取是建立農(nóng)田蒸散量分布圖必須解決的問題。現(xiàn)存遙感蒸散模型大多基于衛(wèi)星數(shù)據(jù)建立，根據(jù)無人機數(shù)據(jù)建立的蒸散模型較少，需要更深入的研究[13-15]。根據(jù)分析通量交換方式的不同，遙感蒸散量模型可分為兩類[16-19]：一是將土壤和植被當做均勻單層進行研究的單層模型，二是考慮土壤和植被水熱傳輸特性差異及相互作用的多層模型。近年來，適用于稀疏植被覆蓋條件下土壤和植被能量交換的多層模型，尤其是雙層模型被廣泛運用[20]。SAMANI等[21]已經(jīng)證明，RSEB模型可用于植被分布不均勻地區(qū)潛熱通量或蒸散量的估算，這是一種適用于植被分布不均勻區(qū)域的蒸散模型。ORTEGA-FARAS等[14]使用RSEB模型估算了巴西橄欖園中的蒸散量，由于對多光譜數(shù)據(jù)利用不足，使用經(jīng)驗公式獲取土壤熱通量，導致適用性非常差，對于農(nóng)田并不適用。

在數(shù)據(jù)獲取方面，反演蒸散量需要溫度數(shù)據(jù)和多光譜數(shù)據(jù)[22-24]，無人機上搭載小型熱像儀和多光譜相機可以獲取相關數(shù)據(jù)。在使用無人機的小型熱像儀時，由于飛行中距離目標遠等原因極易出現(xiàn)較大的測量誤差，從而影響最終的估算結果，通過校正可以有效減小測量誤差[25]。

本文在分析以上方法與模型優(yōu)缺點的基礎上，利用無人機搭載熱像儀獲取的地表溫度(LST)以及多光譜相機獲取的歸一化植被指數(shù)(NDVI)進行蒸散量的估算，提出將RSEB模型的土壤熱通量計算方式進行基于多光譜數(shù)據(jù)的改進，并對熱像儀數(shù)據(jù)進行溫度校正，建立農(nóng)田蒸散量分布圖，最后，利用渦度系統(tǒng)的通量觀測數(shù)據(jù)驗證結果的準確度，為開展大范圍的農(nóng)田水分監(jiān)測和精準農(nóng)業(yè)灌溉提供支持。

1 材料和方法

1.1 研究區(qū)概況

研究區(qū)位于中國旱區(qū)節(jié)水農(nóng)業(yè)研究院的小麥種植試驗場(34°17′58.47″N，108°4′2.93″E，海拔525 m)，年均蒸發(fā)量1 500 mm，土壤質地為中壤土，1 m土層田間持水率為23%～24%(質量含水率)，凋萎含水率為11%～12%(質量含水率)，平均干容重為1.4 g/cm3。試驗期田間種植冬小麥，在2018年10月播種，試驗區(qū)施肥和灌溉水平與該地區(qū)大田一致。試驗區(qū)地形開闊，附近多為農(nóng)田，氣象監(jiān)測站位于麥田中央。

1.2 硬件系統(tǒng)

基于蒸散模型估算蒸散量的數(shù)據(jù)需求進行硬件系統(tǒng)搭建，硬件系統(tǒng)平臺為M100型大疆多旋翼無人機，無人機質量2 431 g(含電池)，最大起飛質量3 600 g，最大可承受風速10 m/s，最大航行速度22 m/s,無人機平臺如圖1所示。使用Pix4Dcapture進行無人機任務規(guī)劃，飛行高度為15 m。在無人機底部通過云臺搭載熱像儀和多光譜相機，頂部為全球定位系統(tǒng)(GPS)模塊和光照傳感器。

圖1 無人機平臺

試驗所需的溫度數(shù)據(jù)由熱像儀獲得，熱像儀型號為FLIR VueProR，紅外圖像分辨率為640像素×512像素，測量范圍-20～50℃，保存的數(shù)據(jù)為TIFF格式圖像[26]。飛行前使用地面站設備通過藍牙連接熱像儀進行輻射校準。

試驗所需的多光譜數(shù)據(jù)由多光譜傳感器獲得，多光譜傳感器為美國Micasense RedEdge多光譜成像儀[27],鏡頭焦距5.5 mm，視場角47.2°。波段數(shù)5個,分別為藍(中心波長475 nm)、綠(中心波長560 nm)、紅(中心波長668 nm)、紅邊(中心波長717 nm)和近紅外(中心波長840 nm)；單幅影像分辨率為1 280像素×960像素。另配有日光照度計和GPS模塊,總質量僅170 g，適合小型無人機搭載，飛行前通過WiFi連接多光譜傳感器，使用光譜標準面板的圖像作為參考進行輻射校準。

1.3 數(shù)據(jù)獲取和處理

1.3.1氣象數(shù)據(jù)和冠層數(shù)據(jù)

在具有代表性的冬小麥夏玉米輪作旱田下墊面安裝渦度相關系統(tǒng)(Open path eddy covariance，OPEC)自動觀測儀器,如圖2所示，該系統(tǒng)主要組成部分為：CAST3A型三維超聲風速儀、LI75600A型開路CO2/H2O分析儀、CR1000型數(shù)據(jù)采集器、HMP-60型空氣溫濕度探頭、LI200SZ型輻射量表和HFP01SC型熱通量板等。

圖2 渦度相關系統(tǒng)

SunScan冠層分析儀是一款簡便測量和分析冠層中入射和透射光合有效輻射(PAR)的系統(tǒng)，提供關于影響田間作物生長限制因素的有價信息。研究期間使用SunScan在麥田中隨機取4個1 m×1 m的方塊，每個方塊在不同位置測量9次取平均值，每7 d測量葉面積指數(shù)(LAI)。

1.3.2低空遙感數(shù)據(jù)獲取

歸一化差異植被指數(shù)(NDVI)的計算公式為[28]

(1)

式中NIR——近紅外波段反射值

R——紅光波段反射值

相機和無人機之間的時間需同步，以便將記錄的GPS和旋轉位置與每個圖像鏈接，為后續(xù)圖像拼接提供支持。使用Pix4Dcapture在移動設備上進行飛行任務設置，由于熱紅外相機不具備定位設備，熱紅外相機與多光譜相機設置相同的拍照頻率和開啟時間，使用同一個地理坐標數(shù)據(jù)，本研究中拍照間隔為1 s。

為了保證數(shù)據(jù)清晰并去除邊緣影響，設置圖像數(shù)據(jù)區(qū)域比研究區(qū)域邊界多50 m，并在數(shù)據(jù)處理時將多余部分去除。對原始遙感數(shù)據(jù)坐標進行校正，通過GIS技術對數(shù)據(jù)進行配準處理[29]，選取同一時間拍攝的LST和NDVI圖像，首先對兩幅圖像進行特征提取得到特征點；通過進行相似性度量找到匹配的特征點對；然后通過匹配的特征點對得到圖像空間坐標變換參數(shù)；最后由坐標變換參數(shù)進行圖像配準。

2019年夏季，每2 d取一組數(shù)據(jù)，其中在陰雨天氣停止飛行任務。對重疊率進行對照試驗，航向、行間重疊率在60%～95%時每隔5%取值，分析拼圖后的準確度，取拼圖效果最佳的值為試驗標準。結果表明，在大田環(huán)境中，多光譜傳感器在航向重疊率85%、行間重疊率75%時效果最優(yōu)，熱像儀航向重疊率85%、行間重疊率95%時效果最優(yōu)。

每組數(shù)據(jù)包括不同重疊率的兩次飛行，在溫度變化比較明顯的日期適當增加飛行任務。每次飛行都要采集溫度與多光譜數(shù)據(jù)，在5月和6月共飛行80次。飛行時間如表1所示，表中時間采用中國標準時間。

表1 UAV飛行任務

1.3.3無人機數(shù)據(jù)校正

考慮到溫度作為模型中重要的原始數(shù)據(jù)，必須保證其準確性以減小誤差，參考METRIC模型中冷熱點的校正方法，在有代表性的區(qū)域取點將實際值與熱像儀返回的測量值進行對比，通過建立的模型對地表溫度進行校正，減少由于設備和輻射等外界因素產(chǎn)生的系統(tǒng)誤差。

METRIC模型中通過選擇圖像上的“極冷點”(植被覆蓋度高,地表溫度最低)和“極熱點”(植被覆蓋度低,地表溫度最高)求解線性回歸系數(shù)。通過假設地表溫度與溫差的線性關系進行參數(shù)標定。

本研究在溫度分布線性關系的基礎上進行溫度校正。根據(jù)溫度分布的線性關系，選取能代表溫度分布情況的地點進行標定[30-32]，選取的點應當具有差異性、代表性，為溫度校正提供支持。測量極干裸土地表的溫度(裸土)、極干完全植被覆蓋地表的溫度(雜草)、土壤水分充足條件下的裸土(土壤)和完全植被覆蓋地表的溫度(冠層)這4種溫度代表地區(qū)的溫度值。使用UT302C型工業(yè)手持測溫儀進行溫度測量，并利用地表輻射溫度和組分溫度線性組合關系得出溫度的變化規(guī)律。對實際溫度與測量溫度進行建模，進而校正通過熱像儀獲得的地表溫度。

1.4 蒸散量模型

METRIC模型和RSEB模型估算地表瞬時蒸散量是通過計算地表能量平衡相關組分獲得潛熱通量LE，進而估算蒸散量，計算公式為[13,21]

(2)

式中Rn——凈輻射通量，W/m2

G——土壤熱通量，W/m2

H——顯熱通量，W/m2

LE——潛熱通量，W/m2

ET——蒸散量，mm/h

ρw——水密度，kg/m3

δ——蒸發(fā)潛熱，J/kg

1.4.1METRIC模型

METRIC模型優(yōu)點在于利用地表輻射溫度代替空氣動力學溫度,將顯熱通量H與地表輻射溫度和空氣溫度之間的差值通過空氣動力學阻抗建立聯(lián)系，對溫度進行內(nèi)部調(diào)整，但是缺點在于只在植被覆蓋均勻處才能有一個較好的結果，農(nóng)田中植被和裸土交錯分布對結果產(chǎn)生極大的影響。

METRIC模型中Rn計算公式為[14]

Rn=(1-α)Rs↓+RL↑-RL↓-(1-ε)RL↓

(3)

式中α——地表反射率

Rs↓——下行太陽短波輻射，W/m2

RL↑——上行太陽長波輻射，W/m2

RL↓——下行太陽長波輻射，W/m2

ε——地表比輻射率

土壤熱通量在熱量平衡中是一個相對較小的量，直接計算較為困難，METRIC模型估算土壤熱通量和顯熱通量公式為[14]

(1-0.98NDVI4)Rn

(4)

(5)

(6)

(7)

式中Ts——地表溫度，K

ρ——空氣密度，kg/m3

CPair——空氣定壓比熱容，J/(g·K)

z1——裸地地表粗糙度，m

z2——氣象站的高度，m

ΔT——z1與z2處的空氣溫度差，K

rah——熱量傳輸空氣動力學阻抗，s/m

u*——像元的摩擦風速，m/s

k——馮卡門常數(shù)

u200——混合層風速，m/s

zom——像元地表粗糙度，m

1.4.2RSEB模型

RSEB模型[18]可用于估算地表為部分植被覆蓋地區(qū)的潛熱通量。RSEB模型為雙層模型，和單層模型相比可以更準確描述植被分布不均勻地區(qū)的地表通量。RSEB模型的Rn計算公式為[14,22]

Rn=frRnc+(1-fr)Rns

(8)

其中

Rnc=(1-αc)Rse+Lin-Loutc-

(1-εc)Lin

(9)

Rns=(1-αs)Rse+Lin-Louts-

(1-εs)Lin

(10)

式中Rnc——冠層凈輻射分量，W/m2

Rns——土壤凈輻射分量，W/m2

fr——最低覆蓋率

Rse——瞬時入射短波輻射，W/m2

αs——土壤反照率

αc——土壤和冠層反照率

Lin——瞬時傳入長波輻射，W/m2

Loutc——冠層瞬時出射長波輻射，W/m2

Louts——土壤瞬時出射長波輻射，W/m2

εc——冠層表面熱發(fā)射率

εs——土壤表面熱發(fā)射率

G=0.323 6Rn-51.52

(11)

顯熱通量的計算公式為

H=frHc+(1-fr)Hs

(12)

(13)

(14)

式中Hc——冠層顯熱通量分量，W/m2

Hs——土壤顯熱通量分量，W/m2

rac——冠層上方空氣動力學阻抗，s/m

raa——冠層與參考水平之間的空氣動力學阻抗，s/m

ras——土壤與冠層之間的空氣動力學阻抗，s/m

Tc——冠層溫度，K

Ta——空氣溫度，K

1.4.3土壤熱通量的RSEB模型改進

無人機多光譜遙感獲得的數(shù)據(jù)量大而且光譜分辨率高，可以有效表現(xiàn)出地物特性。ORTEGA-FARAS等[14]使用RSEB模型估算蒸散量時沒有引入多光譜數(shù)據(jù)，G通過經(jīng)驗公式得到，由于使用環(huán)境只種植單一的植被，所以不具有普遍適用性，只適用于特定環(huán)境。

為了使估算方法更好地適用于大田試驗環(huán)境，分別考慮METRIC和RSEB模型的優(yōu)點，因為METRIC模型對土壤熱通量的計算是基于多光譜數(shù)據(jù)進行的，可以反映不同環(huán)境下的下墊面差別，所以選擇式(4)進行土壤熱通量的計算[24]，使土壤熱通量的估算能夠適用于不同的環(huán)境。

引入NDVI數(shù)據(jù)計算土壤熱通量可以更好地體現(xiàn)空間的差異性，更好地適用于農(nóng)田尺度?；诙喙庾V數(shù)據(jù)改進的RSEB模型比原始模型適用于更復雜的環(huán)境。

1.5 精度驗證

通過平均絕對誤差(MAE)、均方根誤差(RMSE)和決定系數(shù)(R2)3種指標驗證數(shù)據(jù)的可靠性，分別分析不同蒸散量模型和溫度校正的結果。MAE和RMSE越小，R2越接近1，說明結果越準確。

2 結果與分析

剔除飛行任務中由于環(huán)境和設備影響而產(chǎn)生的異常值，選取其中26組飛行測量數(shù)據(jù)進行對比試驗，選取數(shù)據(jù)的時間序列是在試驗期間均勻分布的。將得到的遙感數(shù)據(jù)進行拼接裁剪等預處理，之后導入GIS中進行計算，得到如下數(shù)據(jù)：NDVI分辨率約為0.01 m，LST分辨率約為0.1 m，最終凈輻射、土壤熱通量、顯熱通量和潛熱通量的空間分布圖以0.1 m的分辨率獲得。

2.1 無人機估算與OPEC測量值對比

2.1.1METRIC模型與RSEB模型比較

試驗區(qū)為大量冬小麥覆蓋的規(guī)則田地，在最初的METRIC模型中對地表溫度的計算是采用熱紅外輻亮度通過對比輻射率的校正所得[14]，本研究中使用熱圖像直接得到LST信息。

圖3 G的測量值與估算值比較結果

由于G在兩種模型中采用的計算方式相同[24]，通過與OPEC系統(tǒng)測量值進行對比，驗證結合多光譜數(shù)據(jù)計算農(nóng)田中G的準確度，對比結果如圖3所示。對數(shù)據(jù)分析后得出，RMSE為5.192 W/m2，MAE為4.511 W/m2，R2為0.94，G的估算結果接近測量值，這種對于G的估算方式適用于無人機數(shù)據(jù)和大田環(huán)境。

然后對比H和LE，LE可以直觀地反映蒸散量的大小[13]，將單層模型的計算結果與雙層模型計算結果分別與渦度系統(tǒng)測量數(shù)據(jù)對比，結果如圖4所示。

圖4 LE和H的測量值與估算值比較結果

分別對METRIC模型和RSEB模型進行分析，RMSE和MAE如表2所示，可以發(fā)現(xiàn)，METRIC模型相對于RSEB模型在估算LE時估算值平均超過實際值20 W/m2，這是因為單層模型在估算時沒有區(qū)分土壤和植被的差異。

表2 METRIC模型和RSEB模型計算結果誤差

在估算LE時，METRIC模型的R2為0.779，RSEB模型的R2為0.788，說明模型估算結果與真實值一致且RSEB模型的結果更為接近，誤差分析如表2所示。

雙層模型的結果一定程度上優(yōu)于單層模型，這是因為雙層模型考慮了冠層與土壤的差別，但總體結果還是有一定的誤差，除了環(huán)境因素產(chǎn)生的誤差外主要誤差來自熱像儀的測溫誤差，飛行的高度和速度都會對其造成影響[26]。通過校正對原始數(shù)據(jù)進行優(yōu)化，以減小誤差，使估算結果能夠更好地接近測量結果。

2.1.2溫度校正對結果的影響

由于飛行任務耗時較長，考慮到溫度隨時間變化，所以飛行前后對不同位置選取點分別取3次數(shù)據(jù)的平均值作為實測值，對熱像儀與測溫儀的值進行比較，平均誤差為1.73℃，部分溫度數(shù)據(jù)示例如表3所示。

表3 實際溫度與熱像儀測量值比較結果

在試驗區(qū)域內(nèi)對經(jīng)過校正后的LST隨機取點與實測值進行比較，平均誤差為0.65℃，說明校正方法對于熱像儀數(shù)據(jù)是可靠的。溫度分布圖如圖5所示，中心部分設備為渦度系統(tǒng)。

圖5 熱像儀測量的地表溫度和校正后地表溫度

根據(jù)溫度校正改進模型,先估算顯熱通量H，進而推算潛熱通量LE，將經(jīng)過校正的數(shù)據(jù)與未經(jīng)過校正的數(shù)據(jù)結果分別與渦度系統(tǒng)測量數(shù)據(jù)對比，結果如圖6所示。

圖6 瞬時LE和H的測量值與估計值比較結果

經(jīng)過溫度校正后，估算LE時METRIC模型的R2為0.848，RSEB模型的R2為0.940，說明經(jīng)過溫度校正后RSEB模型的結果與真實值相關性是最好的，誤差分析如表4所示。

對比數(shù)據(jù)明顯看出，校正溫度數(shù)據(jù)后的H與LE的RMSE比校正溫度數(shù)據(jù)前平均降低了9 W/m2，MAE平均降低了6 W/m2，R2平均提升了0.12，校正后的估算結果更為合理準確，說明校正方法對于無人機的熱紅外數(shù)據(jù)非常有效，經(jīng)過校正后的溫度數(shù)據(jù)更貼合實際溫度。本研究校正方法可以減小無人機熱像儀數(shù)據(jù)獲取過程中的誤差，使結果更為準確。研究結果表明，溫度校正后的RSEB模型估算得到的通量值和監(jiān)測得到的通量值有很高的一致性，R2為0.940，RMSE為40.202 W/m2，MAE為26.017 W/m2, 可以通過本研究方法獲取農(nóng)田中的LE。

表4 校正后的METRIC和RSEB計算結果誤差

2.2 蒸散量空間分布圖

將獲取的數(shù)據(jù)代入模型對遙感數(shù)據(jù)進行計算，NDVI如圖7a所示，得到凈輻射通量Rn、顯熱通量H、土壤熱通量G數(shù)據(jù)分別如圖7b～7d所示。

圖7 NDVI和通量的空間分布圖

根據(jù)式(1)和求得的通量數(shù)據(jù)，進而求得潛熱通量LE如圖8所示，通過估算得到的潛熱通量與測量的潛熱通量非常一致，生成熱通量的空間分布圖進而得到蒸散量的空間分布圖，如圖8b所示。

圖8 潛熱通量和蒸散量的空間分布圖

通過現(xiàn)有其他類型技術很難得到相同結果，渦動協(xié)方差系統(tǒng)測量不能得到空間尺度差異性的數(shù)據(jù)[8]，使用衛(wèi)星數(shù)據(jù)的遙感估算模型無法得到空間分辨率較高數(shù)據(jù)，并且受到時間條件約束[9]。通過無人機數(shù)據(jù)最終得到ET如圖8b所示，圖中可以清楚地看到不同位置蒸散量的差異，本研究結果在空間尺度上表示的農(nóng)田蒸散量可以達到分米級。

3 結論

(1)搭建了獲取蒸散模型基礎數(shù)據(jù)的無人機平臺，并將無人機數(shù)據(jù)代入蒸散模型進行研究，比較METRIC模型與RSEB模型在農(nóng)田環(huán)境下的結果，對RSEB模型的土壤熱通量計算方式進行基于多光譜數(shù)據(jù)的改進，校正無人機溫度數(shù)據(jù)，最終得到農(nóng)田蒸散量分布圖。

(2)對比典型單層模型和典型雙層模型在農(nóng)田中的適用性發(fā)現(xiàn)，雙層模型的結果更優(yōu)。在通量計算方面，對RSEB模型的土壤熱通量計算方式進行基于多光譜數(shù)據(jù)的改進，使模型的適用性更強。在LST獲取方面，對無人機測溫方式進行優(yōu)化，基于FLIR VueProR熱像儀多次試驗，取最佳飛行重疊率與高度，并進行溫度校正，即選取有代表性的溫度點，分別位于裸地、冠層、雜草、土壤，隨機采集5個點的真實溫度，與無人機全局測溫圖進行建模，獲得校正后溫度。溫度校正有效減小了測量中產(chǎn)生的誤差，使熱像儀提供的數(shù)據(jù)更為準確。

(3)將無人機遙感、渦度相關技術、地理信息技術相結合，通過無人機獲取溫度和多光譜數(shù)據(jù)，GIS軟件進行數(shù)據(jù)處理，建立有效的蒸散模型，獲取蒸散量結果。將結果與渦度系統(tǒng)測量真實數(shù)據(jù)進行比較，結果表明，基于無人機與氣象數(shù)據(jù)建立的蒸散模型適用于大田作物，空間分辨率較高(分米級)，能夠實現(xiàn)蒸散量的反演。