孫華生, 張 遠(yuǎn), 史云飛, 趙 敏

1. 臨沂大學(xué)資源環(huán)境學(xué)院(山東省水土保持與環(huán)境保育研究所)， 山東 臨沂 276000

2. 華東師范大學(xué)地理科學(xué)學(xué)院， 崇明生態(tài)研究院， 地理信息科學(xué)教育部重點(diǎn)實(shí)驗(yàn)室， 上海 200241

引 言

無人機(jī)遙感具有操作簡便、 數(shù)據(jù)獲取成本低、 時(shí)空分辨率高、 可在云下實(shí)現(xiàn)數(shù)據(jù)采集等諸多優(yōu)勢， 彌補(bǔ)了衛(wèi)星或大型航空遙感系統(tǒng)的缺點(diǎn)， 所以在小范圍、 高時(shí)間和高空間分辨率遙感監(jiān)測中， 發(fā)揮了極為重要的作用。 目前， 無人機(jī)遙感技術(shù)已被廣泛應(yīng)用于農(nóng)業(yè)、 林業(yè)、 環(huán)境領(lǐng)域的定量監(jiān)測[1]。 反射率是地物較為穩(wěn)定的物理屬性， 也是定量遙感模型的重要輸入?yún)?shù)之一。 在遙感圖像定量化應(yīng)用過程中， 需要將相機(jī)記錄的數(shù)字量化(digital number， DN)值轉(zhuǎn)換為地表反射率。 因此， 得到準(zhǔn)確的地表反射率是實(shí)現(xiàn)遙感數(shù)據(jù)定量化應(yīng)用的基礎(chǔ)。

在環(huán)境光照穩(wěn)定的條件下， 通常只要在無人機(jī)起飛前或航拍完成后， 拍攝兩個(gè)或更多具有不同反射率的參考板， 通過對每個(gè)波段記錄的DN值進(jìn)行線性拉伸， 即可得到每個(gè)波段的地表反射率[2-4]。 然而， 實(shí)際的狀況卻是異常復(fù)雜的[5-7]， 這是因?yàn)椋?(1)受云層遮擋和太陽高度角變化的影響， 數(shù)據(jù)獲取時(shí)的光照條件可能會(huì)隨時(shí)發(fā)生變化； (2)為了得到高質(zhì)量的圖像， 相機(jī)曝光時(shí)間、 光圈大小和ISO值， 需要根據(jù)拍攝時(shí)的光照條件自動(dòng)調(diào)整； (3)在很多情況下， 在地面上鋪設(shè)參考板是很難或無法實(shí)現(xiàn)的(如森林、 水域或沼澤地中)， 而且要想在每幅圖像中都拍攝到一組參考板也是不現(xiàn)實(shí)的。

為了擺脫對參考板的依賴， 一些研究者采用經(jīng)驗(yàn)線法得到近似的地表反射率[2]； 還有的建議在無人機(jī)上安裝一個(gè)光強(qiáng)傳感器， 來同步記錄相機(jī)拍攝時(shí)每個(gè)波段對應(yīng)的光照強(qiáng)度， 并將該測量結(jié)果用于計(jì)算地表反射率[3-5, 7-9]， 所以獲取可靠的地面接收的太陽輻照度是得到準(zhǔn)確地表反射率的前提。 如果光強(qiáng)傳感器不是安裝在一個(gè)精密的增穩(wěn)平臺(tái)上， 受飛行器姿態(tài)變化的影響， 導(dǎo)致傳感器的測量結(jié)果與到達(dá)地面的太陽輻射強(qiáng)度并不一致。 對此， 一些研究者采用區(qū)塊校正法來降低傳感器的傾斜效應(yīng)[8-10]， 但其得到的只是一種近似的結(jié)果。 還有采用余弦函數(shù)對直接測量的結(jié)果進(jìn)行傾斜校正， 但實(shí)驗(yàn)測試表明， 余弦校正的結(jié)果并不理想[3, 6]。 這是因?yàn)椋?太陽輻射不僅包括直射， 還包括散射， 在不同天氣狀況下， 太陽輻射的組成有著很大的差異， 而傳感器傾斜對太陽直射和散射的影響是完全不同的。 為了獲取更加準(zhǔn)確的太陽輻照度測量結(jié)果， 可將光強(qiáng)傳感器安裝在一個(gè)自動(dòng)調(diào)平裝置上。 例如， Markelin等利用一個(gè)光學(xué)平衡系統(tǒng)， 可以對光強(qiáng)傳感器進(jìn)行15°以內(nèi)的傾斜補(bǔ)償[4]， 但對多旋翼無人機(jī)來說， 在飛行狀態(tài)下， 其傾斜角度有時(shí)甚至?xí)^30°。 如果只是將光強(qiáng)傳感器搭載在一個(gè)簡易的調(diào)平裝置上， 并不能完全消除飛行器傾斜造成的影響； 而如果單獨(dú)為光強(qiáng)傳感器安裝一個(gè)精密的調(diào)平裝置(如一個(gè)三軸增穩(wěn)云臺(tái))， 將會(huì)大大增加飛行器的負(fù)擔(dān)和制造成本。 在實(shí)際應(yīng)用中， 在沒有精密調(diào)平裝置的條件下， 為了使光強(qiáng)傳感器的測量結(jié)果更加可靠， 常用的做法包括： (1)要求航拍時(shí)太陽高度角較大， 并且航線與太陽方位角垂直， 以降低傳感器傾斜對太陽輻照度測量結(jié)果的影響， 但這種方法并不能消除傳感器傾斜的影響； (2)采用懸停方式進(jìn)行拍攝， 但這種方式只有在沒有風(fēng)的情況下才能使傳感器保持水平， 否則， 并不能完全消除傳感器傾斜的影響， 而且其作業(yè)效率極低。 因此， 目前利用無人機(jī)光強(qiáng)傳感器對太陽輻照度的測量， 仍然存在一定的缺陷， 或者只能在一定的限定條件下使用。

如何在傾斜狀態(tài)下得到準(zhǔn)確的太陽輻照度， 進(jìn)而得到準(zhǔn)確的地表反射率， 是無人機(jī)遙感領(lǐng)域多年來一直未能很好解決的技術(shù)難題。 本研究對該問題進(jìn)行了深入探討， 并提出了一種太陽輻照度準(zhǔn)確測量的新方法。 該方法可以在任意光照條件下， 直接將多光譜圖像記錄的DN值轉(zhuǎn)換為準(zhǔn)確的地表反射率。

1 實(shí)驗(yàn)部分

1.1 無人機(jī)平臺(tái)和傳感器

利用大疆精靈4 RTK無人機(jī)作為遙感數(shù)據(jù)獲取平臺(tái)。 該無人機(jī)搭載了一臺(tái)一體式的多光譜相機(jī)， 共有6個(gè)1/2.9 英寸的CMOS傳感器(詳見圖1)， 包含1個(gè)可見光波段的RGB傳感器和5個(gè)多光譜傳感器。 5個(gè)多光譜波段參數(shù)的描述， 如表1所示。 此外， 在該無人機(jī)的頂部還集成了一個(gè)光強(qiáng)傳感器， 可同步檢測5個(gè)波段的太陽輻照度， 并將其記錄在圖像的XMP(Extensible Metadata Platform)元數(shù)據(jù)中[11]。

1.2 多光譜圖像轉(zhuǎn)換為地表反射率的方法

地表反射率的計(jì)算方法如式(1)所示

Rλ=Er(λ)/Eg(λ)

(1)

式(1)中，Rλ為波長為λ的地表反射率，Er(λ)為地物對波長為λ的電磁波在單位波長單位面積的反射輻射功率(W·m-2·nm-1)，Eg(λ)為到達(dá)地面的入射輻射輻照度(W·m-2·nm-1)。

對Er(λ)的定量描述如式(2)所示

(2)

式(2)中，L為每個(gè)波段每個(gè)像元對應(yīng)的地物反射輻照度(W·m-2·Sr-1·nm-1)，θ為太陽輻射的天頂角，φ為太陽輻射的方位角。

在計(jì)算反射率時(shí)， 需要首先將DN值轉(zhuǎn)換為對應(yīng)像元的輻射亮度， 然后再將輻射亮度進(jìn)一步轉(zhuǎn)換為地表反射率。 由于無人機(jī)多光譜相機(jī)拍攝的原始圖像記錄的DN值會(huì)受曝光時(shí)間、 ISO值、 暗角效應(yīng)等因素的影響， 所以需要事先對原始圖像記錄的DN值進(jìn)行修正處理， 以得到統(tǒng)一標(biāo)準(zhǔn)下的結(jié)果。 具體轉(zhuǎn)換方法如式(3)所示

DN′=(DN-DNblack)/N/Gsensor/Te×A×V

(3)

式(3)中，DN′為轉(zhuǎn)換后的每個(gè)波段的數(shù)字量化值，DN為原始圖像記錄的每個(gè)波段的數(shù)字量化值，DNblack為黑電平的數(shù)字量化值，N為圖像可表示的最大灰度值，Gsensor為傳感器增益值，Te為曝光時(shí)間，A為傳感器敏感度修正系數(shù)，V為相機(jī)鏡頭的暗角效應(yīng)修正值。 以上參數(shù)都是在實(shí)驗(yàn)室中測定的， 其中DNblack，Gsensor，Te和A可以從XMP元數(shù)據(jù)中查找到， 而暗角效應(yīng)修正值V可按照式(4)計(jì)算得到

V(x,y)=(1+k1r+k2r2+k3r3+k4r4+k5r5+k6r6)

(4)

修正后的DN′(λ)與每個(gè)波段的輻射亮度L(λ)之間滿足簡單的線性轉(zhuǎn)換關(guān)系， 其表達(dá)式如式(5)所示

L(λ)=Gλ×DN′(λ)+Bλ

(5)

式(5)中，Gλ為增益值，Bλ為偏移值，DN′為修正后的數(shù)字量化值。

因此， 根據(jù)式(1)、 式(2)和式(5)， 地表反射率Rλ可表示為式(6)的形式

Rλ=Er(λ)/Eg(λ)=πL(λ)/Eg(λ)=

π[Gλ×DN′(λ)+Bλ]/Eg(λ)

(6)

如果能夠準(zhǔn)確地得到圖像拍攝時(shí)到達(dá)地面的入射輻射輻照度Eg(λ)， 通過在一張圖像中拍攝到兩張或者更多具有不同反射率的參考板(其反射率是已知的)， 即可得到每個(gè)波段的定標(biāo)參數(shù)Gλ和Bλ， 其計(jì)算方法如式(7)所示

DN′(λ)×Gλ+Bλ=RλEg(λ)/π

(7)

由此可見， 利用傾斜的光強(qiáng)傳感器獲取準(zhǔn)確的Eg(λ)， 是需要解決的關(guān)鍵問題， 而這也是長期以來一直未能很好解決的技術(shù)難題。 在本研究中， 提出了一種實(shí)現(xiàn)太陽輻照度準(zhǔn)確測量的新方法， 對其具體描述如下：

首先， 需要對光強(qiáng)傳感器接收的太陽輻射進(jìn)行定量化描述。 其接收的太陽輻射包括： 太陽直射、 天空散射和少量的地表反射。 對太陽直射來說， 其輻照度測量結(jié)果受傳感器傾斜的影響， 如式(8)所示

(8)

(9)

(10)

根據(jù)以上分析可知， 在光強(qiáng)傳感器平面傾斜時(shí)， 直接測量的每個(gè)波段的太陽輻照度Em的表達(dá)式， 如式(11)所示

(11)

根據(jù)式(11)可知， 如果能夠通過某種方法獲取朝向兩個(gè)或者更多方向的太陽輻照度， 就可以直接求出直射Edir和散射Esca。 需要注意的是： 如果z≥90°， 即cos(z)≤0， 就會(huì)導(dǎo)致直射輻射無法到達(dá)光強(qiáng)傳感器的正面， 從而使傳感器正面處于陰影中。 此時(shí)， 需要令cos(z)=0。

在得到直射Edir和散射Esca之后， 地面接收的太陽輻照度Eg可表示為式(12)的形式

Eg=Edircosθ+Esca

(12)

在得出地面接收的太陽輻照度Eg后， 即可根據(jù)式(6)得到各個(gè)波段每個(gè)像素對應(yīng)的地表反射率Rλ。

實(shí)際上， 如果能夠得到直射輻射占總太陽輻射的比例p(p=Edir/Esun， 而Esun=Edir+Esca)， 那么就能將光強(qiáng)傳感器測量值Em， 直接轉(zhuǎn)換為地面接收的太陽輻照度Eg。 把Edir=pEsun，Esca=(1-p)Esun代入式(11)， 即可得到Em和Esun之間的關(guān)系式， 其結(jié)果如式(13)所示

Em=[pcosz+(1-p)cos2(s/2)+

(13)

因此， 地面接收的太陽輻照度Eg可以進(jìn)一步轉(zhuǎn)換為式(14)的形式

Eg=Edircosθ+Esca=(pcosθ+1-p)Esun=

(14)

1.3 實(shí)驗(yàn)測試

1.3.1 實(shí)驗(yàn)區(qū)

實(shí)驗(yàn)區(qū)位于山東省臨沂市， 其地理位置為35°9′53.5″N—35°10′02.2″N， 118°16′00.5″E—118°16′08.5″E。 該實(shí)驗(yàn)區(qū)的面積約250 m×200 m。 航線規(guī)劃是通過DJI GS Pro應(yīng)用軟件實(shí)現(xiàn)的， 設(shè)定的航高為35 m， 航線重疊度為67%， 旁向重疊度為55%， 共生成14條航線， 每條航線拍攝31張圖像。 飛行速度約為3.8 m·s-1， 每次航拍約持續(xù)16 min。 從2020年7月開始， 一直持續(xù)到11月底， 共獲取了4個(gè)時(shí)期的航拍圖像。

1.3.2 相機(jī)定標(biāo)

選取黑、 灰、 白三張參考板(尺寸為18 cm×25 cm)作為靶標(biāo)， 以實(shí)現(xiàn)無人機(jī)多光譜圖像的反射率定標(biāo)和反射率測量結(jié)果的驗(yàn)證。 各個(gè)參考板在不同波段的反射率通過PSR-1100高光譜儀測量得到。 該高光譜儀的波譜覆蓋范圍為320～1 100 nm， 在600 nm處的光譜分辨率為3 nm。 黑、 灰、 白參考板及其高光譜測試結(jié)果， 如圖2所示。 黑、 灰、 白三張參考板在不同波段的平均反射率， 如表2所示。

表2 利用高光譜儀測量的參考板的平均反射率

圖2 黑、 灰、 白參考板及其高光譜測試結(jié)果

在得到無人機(jī)拍攝的參考板的多光譜圖像之后， 根據(jù)式(7)構(gòu)建Rλ和DN′(λ)之間的線性拉伸關(guān)系， 即可得到每個(gè)波段的Gλ和Bλ的結(jié)果。 對一臺(tái)多光譜成像系統(tǒng)來說， 一旦得到其定標(biāo)參數(shù)， 該定標(biāo)參數(shù)通常是較為穩(wěn)定的， 所以在實(shí)際應(yīng)用中就不必再依靠參考板， 即可得到準(zhǔn)確的地表反射率。

1.3.3 驗(yàn)證方案設(shè)計(jì)

如果飛行器自身攜帶兩個(gè)或者更多的朝向不同方向的光強(qiáng)傳感器， 即可根據(jù)式(11)實(shí)現(xiàn)直射Edir和散射Esca的求解。 然而， 由于本研究采用的無人機(jī)只有一個(gè)光強(qiáng)傳感器， 而且是固定在無人機(jī)頂部的， 所以無法直接實(shí)現(xiàn)以上解算。 為了驗(yàn)證本方法的可行性， 設(shè)計(jì)了如下的實(shí)驗(yàn)方案：

首先， 選取一塊開闊的實(shí)驗(yàn)場地， 確定其周圍沒有高大的樹木或建筑物。 然后， 用手舉起飛行器并手動(dòng)調(diào)整其姿態(tài)， 以讓光強(qiáng)傳感器平面朝向不同的方向(在操作時(shí)可以使傳感器平面分別處于水平、 前傾、 后傾、 左傾、 右傾的狀態(tài))， 并拍照以記錄不同方向上各波段接收的太陽輻照度， 從而可實(shí)現(xiàn)對各個(gè)波段的直射Edir和散射Esca的求解。 為了檢驗(yàn)得到的地表反射率的準(zhǔn)確度， 在實(shí)驗(yàn)場中還布置了一組黑、 灰、 白參考板。

分別在不同日期不同光照條件下獲取航拍數(shù)據(jù)， 而每次實(shí)驗(yàn)時(shí)的天氣都較為穩(wěn)定。 在這種情況下， 在一段較短的時(shí)間內(nèi)， 可認(rèn)為太陽輻射的組分是不變的(即使太陽輻照度有小幅度的變化， 但其組分的變化較小)。 因此， 可利用本方法， 得到在航拍前或航拍后的直射比例p， 并根據(jù)式(14)計(jì)算出所有圖像在拍照時(shí)刻地面接收的太陽輻照度Eg， 從而可根據(jù)該結(jié)果得到所有圖像的地表反射率。 最后， 利用參考板來驗(yàn)證計(jì)算結(jié)果的準(zhǔn)確度。

2 結(jié)果與討論

由于每次實(shí)驗(yàn)都是在光照條件較為穩(wěn)定的情況下進(jìn)行的(但不同日期的光照條件是不同的)， 所以航拍時(shí)段的直射比例可以在航拍前或航拍后， 利用本方法來確定。 航拍時(shí)段不同波段直射比例的計(jì)算結(jié)果， 如表3所示。

表3 航拍時(shí)段不同波段的直射比例p

在得到直射比例p之后， 通過式(14)即可將光強(qiáng)傳感器直接測量的結(jié)果Em， 轉(zhuǎn)換為地面接收的太陽輻照度Eg。 不同日期的具體計(jì)算結(jié)果如圖3所示。

圖3 光強(qiáng)傳感器直接測量的太陽輻照度Em及對其進(jìn)行校正后得到的地面接收的太陽輻照度Eg的結(jié)果

由圖3可以看出， 光強(qiáng)傳感器的直接測量結(jié)果受飛行器自身姿態(tài)變化(因航向發(fā)生了改變)的影響很大， 尤其是當(dāng)太陽高度角較小(即入射角較大)時(shí)， 其影響極為顯著。 利用本研究提出的方法可以很好地消除測量結(jié)果的劇烈波動(dòng)。 為了進(jìn)一步檢驗(yàn)該校正結(jié)果的準(zhǔn)確度， 利用事先鋪設(shè)的黑、 灰、 白參考板， 來檢驗(yàn)無人機(jī)多光譜圖像得到的反射率的準(zhǔn)確度， 其結(jié)果如表4所示。

表4 利用無人機(jī)多光譜圖像得到的參考板在不同波段的平均反射率

利用平均絕對誤差(MAE)， 以及標(biāo)準(zhǔn)差(S)兩種指標(biāo)對測量結(jié)果的準(zhǔn)確度進(jìn)行評價(jià)。 計(jì)算結(jié)果顯示， 藍(lán)(B)、 綠(G)、 紅(R)、 紅邊(RE)和近紅外(NIR)波段在所有日期所有參考板的MAE的結(jié)果分別為3.34%， 1.36%， 1.45%， 1.81%和1.63%， 而各自的標(biāo)準(zhǔn)差S分別為2.11%， 0.85%， 1.21%， 1.88%和0.59%。 黑、 灰、 白三張參考板在所有日期的所有波段中的MAE的結(jié)果分別為0.82%， 1.99%和2.94%， 而各自的標(biāo)準(zhǔn)差S分別為0.55%， 1.31%和1.84%。 引起地表反射率測量誤差的因素主要包括以下幾個(gè)方面： (1)參考板不是嚴(yán)格的漫反射體， 得到的反射率會(huì)因太陽入射方向和相機(jī)觀測方向的不同而存在一定的誤差； (2)參考板的高光譜測量結(jié)果， 以及得到的相機(jī)定標(biāo)參數(shù)存在一定的誤差； (3)光照條件不是絕對穩(wěn)定的， 因?yàn)樵诤脚臅r(shí)段內(nèi)太陽高度角和云量發(fā)生了一定的變化， 從而導(dǎo)致真實(shí)的直射比例p不是絕對穩(wěn)定的。

盡管得到的地表反射率存在著一定的誤差， 但從以上統(tǒng)計(jì)結(jié)果可以看出， 在不同光照條件下得到的不同波段的地表反射率， 仍然具有較高的準(zhǔn)確度。 由此可以推斷出， 利用本方法得到的地面接收的太陽輻照度Eg的結(jié)果是可靠的。

3 結(jié) 論

提出了一種全新的可以利用多光譜無人機(jī)圖像實(shí)現(xiàn)地表反射率直接測量的方法。 與傳統(tǒng)的基于光強(qiáng)傳感器實(shí)現(xiàn)太陽輻照度測量方法不同的是： 本方法不是直接消除光強(qiáng)傳感器傾斜的影響， 而是借助傳感器的傾斜效應(yīng)來實(shí)現(xiàn)對光照強(qiáng)度的校正。 即使在光照強(qiáng)度變化的條件下， 利用該方法仍然能夠得到準(zhǔn)確的測量結(jié)果， 從而可據(jù)此進(jìn)一步得到準(zhǔn)確的地表反射率。 本研究解決了多年來一直未能得到很好解決關(guān)鍵技術(shù)問題， 對多光譜無人機(jī)的設(shè)計(jì)， 以及無人機(jī)遙感數(shù)據(jù)的定量化應(yīng)用， 都具有重要意義。 受硬件條件的限制， 本研究并沒有直接利用無人機(jī)搭載的多個(gè)傳感器， 來獲取每個(gè)拍照時(shí)刻的太陽輻照度。 在未來的研究中， 將對飛行器進(jìn)行硬件上的改進(jìn)， 使其攜帶多個(gè)朝向不同方向的光強(qiáng)傳感器， 以實(shí)現(xiàn)每個(gè)拍照時(shí)刻的直射和散射輻照度的準(zhǔn)確測量， 并檢驗(yàn)其實(shí)際測量效果。

附注1: 幾個(gè)角度的計(jì)算方法

(1)太陽直射的天頂角θ可以通過式(A1)來確定

cos(θ)=sin(Ψ)sin(δ)+cos(Ψ)cos(δ)cos(H)

(A1)

式(A1)中，Ψ為某一點(diǎn)的緯度(弧度)；δ為太陽赤緯，δ=-23.44/180πcos(2π/365(N+10))，N為某日期在一年中的天數(shù)；H為時(shí)角(以南為參考)。H可以通過式(A2)計(jì)算出來

H=π/12(Tsolar-12)

(A2)

式(A2)中，Tsolar為地方太陽時(shí)， 即根據(jù)太陽的具體位置所確定的時(shí)刻， 與地方標(biāo)準(zhǔn)時(shí)間TLst的關(guān)系為Tsolar=TLst+12/π(Lonloc-Lonst)+E，Lonloc為某一點(diǎn)的經(jīng)度(弧度)，Lonst為地方標(biāo)準(zhǔn)時(shí)間所采用的標(biāo)準(zhǔn)經(jīng)度(弧度)，E為地球繞太陽公轉(zhuǎn)時(shí)的進(jìn)動(dòng)和轉(zhuǎn)速變化對地方太陽時(shí)的修正，E=[9.87sin(2B)-7.53cos(B)-1.5sin(B)]/60， 而B=2π(N-81)/364。

(2)傳感器傾斜的坡度s可以利用式(A3)來確定

(A3)

式(A3)中，vZ為Z軸的方向向量(vZ=(0, 0, 1)T)；vn為斜面的法線向量，vn=RvZ(R為旋轉(zhuǎn)矩陣， 可以利用俯仰角、 側(cè)滾角和偏航角計(jì)算出來， 它們都記錄在XMP文件中)。

(3)太陽直射在傾斜傳感器坡面上的入射角z可以利用式(A4)來確定

cos(z)=cos(θ)cos(s)+sin(θ)sin(s)cos(φ-a)

(A4)

式(A4)中，θ為太陽直射的天頂角；s為傳感器傾斜的坡度；φ為太陽方位角，a為傳感器傾斜的坡向。 而φ和a可以分別通過式(A5)和式(A6)計(jì)算出來。 注：φ和a的參考方向需要與式(A2)中的時(shí)角H一致， 即以南為參考。

φ=arctan2(sinφ, cosφ)

(A5)

a=arctan2(sina, cosa)+π

(A6)