王玥璞， 秦 偉， 楊文濤， 王玉杰?， 王云琦， 蔣 濤

(1.北京林業(yè)大學(xué)水土保持學(xué)院，100083，北京；2.中國水利水電科學(xué)研究院泥沙研究所，100044，北京)

近年來，無人機(jī)遙感技術(shù)憑借其可操作性強(qiáng)、分辨率高的優(yōu)勢(shì)被廣泛應(yīng)用在植被資源調(diào)查和地形監(jiān)測(cè)工作中[1-2]。其中，基于SfM(Structure from Motion)方法的攝影測(cè)量技術(shù)為地形研究提供新的數(shù)據(jù)獲取手段[3-4]。其數(shù)據(jù)采集方式多基于無人機(jī)平臺(tái)開展，且相較于Li DAR 技術(shù)具有操作簡(jiǎn)單靈活，數(shù)據(jù)獲取效率高，技術(shù)成本低等特點(diǎn)。目前無人機(jī)遙感技術(shù)在水土保持領(lǐng)域得到了廣泛的應(yīng)用。例如無人機(jī)遙感技術(shù)可用于水土保持監(jiān)測(cè)項(xiàng)目中土石方量的精確測(cè)算[5]；該技術(shù)還可為水土保持信息化監(jiān)測(cè)提供服務(wù)等[6]。

近期，有研究發(fā)現(xiàn)，無人機(jī)SfM技術(shù)在某些土地類型的地形獲取精度上具有和激光雷達(dá)相媲美的精度。Liao等[7]在東北較平坦地形對(duì)比同一地塊對(duì)大疆精靈4pro所拍攝的可見光生成的DEM與激光雷達(dá)生成的DEM對(duì)比發(fā)現(xiàn)：2類高程間的相關(guān)性在裸地最高(R2=0.998 7)，荒地(R2=0.965 6)僅次于裸地，林地效果最差(R2= 0.340 5)。林鑫等[8]分析無人機(jī)密集匹配點(diǎn)云和機(jī)載激光雷達(dá)點(diǎn)云的異同性，對(duì)密集林分(郁閉度0.85)、稀疏林分(郁閉度0. 55)和未成林地的 2 種點(diǎn)云對(duì)比分析，發(fā)現(xiàn)對(duì)于稀疏林分，無人機(jī)密度匹配點(diǎn)云的統(tǒng)計(jì)特征參數(shù)與機(jī)載激光雷達(dá)點(diǎn)云的統(tǒng)計(jì)特征參數(shù)基本一致，對(duì)于新造未成林地，無人機(jī)密集匹配點(diǎn)云對(duì)幼樹三維結(jié)構(gòu)的刻畫優(yōu)于機(jī)載激光雷達(dá)點(diǎn)云。研究均表明中覆蓋(FVC<0.6)條件下采用無人機(jī)SfM技術(shù)具有較高精度。

無人機(jī)SfM方法已逐步被應(yīng)用于地形研究領(lǐng)域[9]，但其研究更關(guān)注數(shù)據(jù)本身情況，并未對(duì)數(shù)據(jù)獲取時(shí)間以及地表覆蓋情況能否適用于地形研究進(jìn)行探討。在實(shí)際應(yīng)用中，SfM方法受到以下條件的限制：1)在植被發(fā)育好的區(qū)域，由于光線不能穿透植被，SfM測(cè)量方法難以測(cè)量密閉植被下的地形。2)SfM 測(cè)量技術(shù)還容易受不利天氣的影響。例如，大風(fēng)天氣會(huì)較大程度干擾螺旋槳的氣流，影響小型無人機(jī)控制平臺(tái)的穩(wěn)定，易造成飛行器失穩(wěn)從而影響拍攝相片質(zhì)量。降水天氣飛行時(shí)，雨水可能會(huì)隨著螺旋槳?dú)饬鞯男齽?dòng)進(jìn)入無人機(jī)機(jī)身內(nèi)部從而對(duì)電池或電路造成損壞，致使短路或斷路的情況發(fā)生[10]。因此，開展氣象、植被覆蓋等的分析對(duì)預(yù)先規(guī)劃和順利實(shí)施無人機(jī)外業(yè)地形勘測(cè)任務(wù)至關(guān)重要[11]。筆者以東北黑土區(qū)寒溫帶、中溫帶和暖溫帶為例，通過分析不同溫度帶的溫度、風(fēng)速、降水、植被覆蓋度、雪蓋等條件嘗試挑選出常用無人機(jī)的地形勘察適宜期，以期為水土保持領(lǐng)域野外地形勘測(cè)規(guī)劃提供參考。

1 研究區(qū)概況

試驗(yàn)區(qū)東北黑土區(qū)是世界僅存的“三大黑土區(qū)”之一，地理位置介于E 41°01′～53°05′、N 115°03′～135°05′ 之間,是世界3大黑土區(qū)中面積最小，氣候條件最寒冷的區(qū)域。本文東北黑土區(qū)參照水土保持區(qū)劃分區(qū)[12]，其主要分布在松嫩流域腹地，北起嫩江、北安，南至四平(圖1)。研究區(qū)氣候干燥寒冷，降雨集中，年平均降雨量400～700 mm，水資源總量 1 415.99億m3，人均占有量1 213.98 m3。東北黑土區(qū)從南向北具有暖溫帶、中溫帶和寒溫帶的熱量變化，自東向西具有濕潤、半濕潤和半干旱的濕度分異，具有獨(dú)特的植被分布格局，是全球變化研究的敏感區(qū)域之一。東北黑土區(qū)植被覆蓋度較高，主要植被類型為農(nóng)業(yè)植被、森林和草原，還有少部分草甸、灌叢等植被。

圖1 東北黑土區(qū)與代表氣象站地理位置示意圖Fig.1 Geographical location map of the Black Soil Region of Northeast China and meteorological stations

2 數(shù)據(jù)與方法

2.1 不同氣候帶的氣象數(shù)據(jù)

本文的試驗(yàn)區(qū)主要涉及2個(gè)溫度帶：寒溫帶和中溫帶，但是黑土區(qū)最南端毗鄰暖溫帶[13]。為了反映研究區(qū)從北到南的氣候梯度狀況，從北向南依次選擇了大興安嶺地區(qū)的漠河站、齊齊哈爾市的三家子氣象站和最南端的大連市周水子站，并在這3個(gè)氣象站附近各選擇了一個(gè)哨兵二號(hào)數(shù)據(jù)瓦塊分析地表覆被狀況。氣象數(shù)據(jù)來源于中國氣象科學(xué)數(shù)據(jù)服務(wù)網(wǎng)(http:∥data.cma.gov.cn)，選擇含有溫度、風(fēng)速、降雨量等內(nèi)容的數(shù)據(jù)。時(shí)間范圍為1988—2017年的有效數(shù)據(jù)。

2.2 不同氣候帶植被動(dòng)態(tài)數(shù)據(jù)

植被覆蓋是干擾無人機(jī)可見光攝影測(cè)量獲取地形的重要因素。為了分析適宜無人機(jī)獲取真實(shí)地形的有效窗口，在3個(gè)氣候帶選取3個(gè)瓦塊的哨兵二號(hào)數(shù)據(jù)[14]，其土地利用類型均為農(nóng)地。寒溫帶的哨兵二號(hào)數(shù)據(jù)覆蓋了漠河氣象站，數(shù)據(jù)瓦塊編號(hào)是51UVU；中溫帶的哨兵數(shù)據(jù)覆蓋三家子氣象站附近的鶴山農(nóng)場(chǎng)，數(shù)據(jù)瓦塊編號(hào)是51UXQ；暖溫帶的哨兵數(shù)據(jù)覆蓋周水子氣象站附近，數(shù)據(jù)瓦塊編號(hào)是51TVE。在漠河地區(qū)、鶴山農(nóng)場(chǎng)和周水子地區(qū)，分別獲取了2020年的17幅、51幅和12期無云、高質(zhì)量的哨兵二號(hào)衛(wèi)星影像(影像拍攝時(shí)間請(qǐng)參見表1)。

表1 所使用哨兵二號(hào)影像的拍攝日期

對(duì)所有影像進(jìn)行歸一化植被指數(shù)(normalized difference vegetation index, NDVI)計(jì)算，見式(1)[15]；并在此基礎(chǔ)上，采用最小二乘法對(duì)全年的NDVI時(shí)間序列做了平滑濾波，以消除不利天氣的可能影響[16]。然后根據(jù)式(2)進(jìn)行植被覆蓋度(fractional vegetation cover, FVC)計(jì)算。

根據(jù)文獻(xiàn)[17]提到的FVC的計(jì)算公式，式中NDVI的最大值和最小值應(yīng)該來自每景影像中的完全被植被覆蓋和沒有任何植被的像素位置。為了簡(jiǎn)化計(jì)算，本文中的NDVI最大和最小值來自固定像素上NDVI時(shí)間序列的最大和最小值。我們選取的位置主要是農(nóng)地，從全年來看其最大值和最小值也表征了該像素位置完全被植被覆蓋和沒有任何植被的兩種情況。東北研究區(qū)有密閉的常綠林分布，完全被該林種覆蓋(100 %覆蓋)的像素上的NDVI在年內(nèi)數(shù)值較穩(wěn)定且應(yīng)與完全被農(nóng)作物覆蓋的NDVI接近。從這種意義上講，筆者采用的從農(nóng)田全年NDVI時(shí)間序列中選取最大和最小NDVI計(jì)算FVC的方法存在合理性。

(1)

式中：rnir和rred分別為影像的近紅外和紅光波段，μm。

(2)

式中：NDVImax和NDVImin分別為給定像素上最大和最小的NDVI值，范圍為-1～1。

此外地表積雪也是影響無人機(jī)搭載可見光傳感器勘測(cè)地形的一個(gè)重要因素。為了反映研究區(qū)積雪的影響，對(duì)所有影像進(jìn)行了歸一化雪蓋指數(shù)(normalized difference snow index, NDSI)[18]計(jì)算。

(3)

式中：Green和SWIR1分別為影像的綠光和中紅外波段，μm。

2.3 無人機(jī)承受的氣象條件

根據(jù)無人機(jī)結(jié)構(gòu)分類，可簡(jiǎn)單分成固定翼式無人機(jī)、旋翼式及多旋翼式無人機(jī)、撲翼式無人機(jī)、傘翼無人機(jī)和其他特形無人飛行器[19]。而水保行業(yè)無人機(jī)多采用多旋翼無人機(jī)。分析水保行業(yè)常用無人機(jī)進(jìn)行平臺(tái)限制因素調(diào)查，得出現(xiàn)有多旋翼無人機(jī)飛行抗風(fēng)能力多在四級(jí)(8 m/s)以下,如表2所示。

表2 水保行業(yè)常用無人機(jī)平臺(tái)運(yùn)行條件

2.4 數(shù)據(jù)分析法

對(duì)各站點(diǎn)的大量數(shù)據(jù)進(jìn)行篩選和分析，具體根據(jù)無人機(jī)適宜溫度和風(fēng)速條件篩選大量數(shù)據(jù)，由于無人機(jī)飛行最高可達(dá)500 m，因此運(yùn)行環(huán)境溫度按照溫度垂直遞減率進(jìn)行計(jì)算，以此分析。此外雨雪天氣則取其30年間同一天降水次數(shù)占30年的比例作為日降雨概率，以估算遭遇降雨天氣的概率。根據(jù)林草植被覆蓋度劃分標(biāo)準(zhǔn)，裸地的植被覆蓋度<0.1，中覆蓋度<0.6。植被覆蓋度則參照前人研究結(jié)果[8]，以中覆蓋作為篩選條件。

3 結(jié)果與分析

3.1 溫度和風(fēng)速適宜性分析

對(duì)3個(gè)典型氣象站1988—2017年氣象數(shù)據(jù)日值資料進(jìn)行日均值計(jì)算，其溫度和風(fēng)速如圖2所示。其中，寒溫帶的漠河站最高溫和最低溫跨度大，最高氣溫均沒有超過40 ℃，但氣溫0 ℃以下時(shí)間較長，占據(jù)全年近50 %(圖2a)。同時(shí)滿足飛行器工作環(huán)境溫度以及飛行至500 m溫度的時(shí)間范圍是4月26日—10月1日。除此以外風(fēng)速在一年中波動(dòng)變化明顯，且呈現(xiàn)雙峰趨勢(shì)。由圖2a可知，平均風(fēng)速均未超過最大可承受風(fēng)速，但最大持續(xù)風(fēng)速一年中有半數(shù)以上超過了8 m/s的限制，滿足飛行條件的時(shí)間范圍為6月29日—9月5日。1—3月和11—12月風(fēng)速雖均在8 m/s以下，但其溫度不符合無人機(jī)飛行器要求，因此不適宜無人機(jī)室外飛行作業(yè)。

中溫帶的三家子站1988—2017年溫度日均值波動(dòng)呈二次函數(shù)曲線變化，最高最低溫差近40 ℃，0 ℃以上時(shí)間較多，對(duì)無人機(jī)野外試驗(yàn)有利(圖2b)。根據(jù)無人機(jī)工作運(yùn)行溫度篩選，滿足其飛行工作溫度的時(shí)間為4月3日—10月23日。對(duì)風(fēng)速數(shù)據(jù)進(jìn)行處理發(fā)現(xiàn)，平均風(fēng)速均在限制條件以下，均符合要求，而最大持續(xù)風(fēng)速的適宜范圍也與溫度相一致。雖有極個(gè)別極端風(fēng)速出現(xiàn)，但時(shí)間不多，且為最大持續(xù)風(fēng)速，可酌情根據(jù)當(dāng)天風(fēng)速情況進(jìn)行作業(yè)。由此可看出中溫帶氣象站點(diǎn)的氣象條件明顯較寒溫帶更適宜無人機(jī)飛行。

對(duì)位于暖溫帶的周水子站30年氣象數(shù)據(jù)進(jìn)行分析發(fā)現(xiàn)，該站點(diǎn)的溫度在0 ℃以上的時(shí)間大于50%，全年氣溫變化趨勢(shì)同樣呈現(xiàn)凸曲線，在8月份溫度達(dá)到最大值(圖2c)。其中同時(shí)滿足飛行器工作環(huán)境溫度以及飛行至500 m溫度的時(shí)間范圍為3月10日—11月25日，明顯較寒溫帶和中溫帶時(shí)間增多。最大持續(xù)風(fēng)速和平均風(fēng)速的變化趨勢(shì)一致，周水子站平均風(fēng)速均未大于8 m/s，最大持續(xù)風(fēng)速在3、4月有部分超過運(yùn)行風(fēng)速，但最大持續(xù)風(fēng)速僅代表當(dāng)天的瞬時(shí)風(fēng)速，具體還應(yīng)參照平均風(fēng)速進(jìn)行判斷，實(shí)際操作過程中應(yīng)避開當(dāng)天風(fēng)速最大的時(shí)間段。綜合而言不適宜溫度和風(fēng)速的時(shí)間較少。

3.2 降水適宜性分析

由于豐水年和平水年情況不同，且降水量的日均值并不具有顯著代表性，故僅對(duì)30年間的日降水概率進(jìn)行統(tǒng)計(jì)。如圖3所示，漠河站降水頻率波動(dòng)較大，在3月至10月呈現(xiàn)明顯的凸曲線，降水頻率較大的時(shí)段多處在夏季，但降水頻率超過0.5的時(shí)間并不多。降水頻率的變化并不能作為估測(cè)某天降雨可能性的切實(shí)依據(jù)，具體選擇應(yīng)參照氣溫和風(fēng)速因素綜合考量，盡可能選在降雨概率小的時(shí)間進(jìn)行無人機(jī)外業(yè)調(diào)查。

圖3 漠河、三家子和周水子氣象站1988—2017年日降雨概率Fig.3 Probability of daily precipitation in 1988-2017 for the Mohe, Sanjiazi and Zhoushuizi meteorological stations

分析中溫帶三家子站1988—2017年日降水頻率發(fā)現(xiàn)，其日降水概率波動(dòng)較大，降水概率在0.5以上的時(shí)間不超過10%，且大多分布在7、8月。無人機(jī)在野外飛行作業(yè)時(shí)應(yīng)選在日降雨頻率小的時(shí)間，提前關(guān)注天氣預(yù)報(bào)，避開陰雨天氣，防止降雨對(duì)無人機(jī)機(jī)身造成損壞。

位于暖溫帶的周水子站30年間日降水概率超出0.5的時(shí)間出現(xiàn)在夏季7月，且時(shí)間很少。春季、秋季和冬季降水概率明顯少于夏季。從5月上旬開始降水頻率逐漸增大，8月初開始逐漸下降。根據(jù)對(duì)日降水頻率進(jìn)行統(tǒng)計(jì)發(fā)現(xiàn)降水頻率并沒有為0的情況，說明降水頻率只能代表其30年間的降水次數(shù)，并不能作為估測(cè)窗口期是否降雨的標(biāo)準(zhǔn)。

3.3 植被覆蓋度、地表雪動(dòng)態(tài)和適宜性分析

從2020年全年NDSI時(shí)間序列可以看出，NDSI值在冬季變大，說明冬季地表有積雪覆蓋，而進(jìn)入春夏秋則保持相對(duì)穩(wěn)定的趨勢(shì)。3類溫度帶中的寒溫帶(漠河)和中溫帶(三家子)的冬季(10月—翌年3月)持續(xù)被積雪覆蓋，不適宜無人機(jī)地形的勘測(cè)，而暖溫帶冬季受雪蓋影響很小。此外，春初秋末還需要考慮突然的降雪影響。以2020年鶴山地區(qū)為例，2020年4月23日積雪指數(shù)突然升高(圖4b)，經(jīng)檢查哨兵影像發(fā)現(xiàn)該區(qū)域出現(xiàn)了大范圍降雪。

圖4 漠河、三家子(以鶴山實(shí)驗(yàn)農(nóng)場(chǎng)表示)和周水子氣象站附近FVC(黑線)和NDSI(橘色線)時(shí)間分布Fig.4 Time series of the FVC (black lines) and NDSI (orange lines) near the Mohe, Sanjiazi (represented by the Heshan Experimental Farm) and Zhoushuizi meteorological station

從植被覆蓋度時(shí)間序列來看，3類溫度帶中的植被覆蓋年際變化曲線總體呈現(xiàn)“幾”字形。植被覆蓋度總體先增加，后趨于穩(wěn)定，最后呈現(xiàn)降低的趨勢(shì)。植被覆蓋度最大值都在夏季7—8月之間達(dá)到最大。對(duì)3類溫度帶而言，5—10月植被蓋度普遍超過0.6，可能對(duì)無人機(jī)SfM地形勘測(cè)有較大影響。綜合地表雪和植被覆蓋兩種因子來看，在春季和秋季共有2個(gè)時(shí)間窗口適宜無人機(jī)SfM地形勘測(cè)。根據(jù)中覆蓋條件進(jìn)行篩選，并綜合考慮NDSI的變化，在3—4月底前及9月底—11月前均可進(jìn)行無人機(jī)飛行。

4 討論

利用無人機(jī)SfM技術(shù)開展地形勘察是以較低成本獲取高精度地形數(shù)據(jù)的有效手段[20]。但是在實(shí)際開展野外地形勘察時(shí)，無人機(jī)平臺(tái)容易受到溫度、風(fēng)速和降水的影響，而SfM技術(shù)容易受到植被覆蓋度和雪蓋的影響，從而難以獲得真實(shí)的地表數(shù)據(jù)。本文嘗試提出一種基于氣候和覆被要素考量的無人機(jī)地形勘察適宜期挑選方法，并對(duì)東北的不同氣候區(qū)開展實(shí)驗(yàn)。通過提前分析研究區(qū)溫度、風(fēng)速、降水、植被蓋度和雪蓋的時(shí)間規(guī)律，有利于指導(dǎo)避開外業(yè)飛行的不利期，對(duì)規(guī)劃和開展野外測(cè)量有一定的現(xiàn)實(shí)意義。此外，筆者所在的研究組曾于2019年5月初在中溫帶的鶴山農(nóng)場(chǎng)附近進(jìn)行無人機(jī)地形測(cè)繪作業(yè)。在野外測(cè)量工作期間其氣象條件和植被覆蓋度條件均與本文結(jié)論較為一致。因此，筆者提出的這種基于氣象和覆被條件室內(nèi)分析方法可有效的指導(dǎo)無人機(jī)的外業(yè)地形勘測(cè)工作。但在具體實(shí)施前仍需要根據(jù)當(dāng)?shù)貙?shí)際的氣象和覆被條件選擇合適的時(shí)間窗口開展野外作業(yè)，以縮減不利天氣的影響、降低對(duì)無人機(jī)的損耗。

考慮到同溫度帶內(nèi)氣候條件的類似性特點(diǎn)，在對(duì)溫度、風(fēng)速等氣象要素的分析上選取的代表性站點(diǎn)的日均值。對(duì)降雨數(shù)據(jù)的分析上，選取了近30年的日降雨頻率進(jìn)行評(píng)估。受可利用數(shù)據(jù)的限制，僅分析了日平均風(fēng)速，但在實(shí)際飛行作業(yè)中瞬間最大風(fēng)速可能更影響無人機(jī)的外業(yè)工作。同一氣候帶的溫度和降水的空間分布規(guī)律一致性較高，但氣候帶對(duì)風(fēng)速的空間指示性較弱。因此，在對(duì)其他地區(qū)的地形勘察適宜期研究中，應(yīng)該著重注意對(duì)風(fēng)速和覆被條件的分析，例如在未來的研究中，應(yīng)該采用基于面狀的氣象預(yù)測(cè)要素開展適宜期分析。

采用的植被覆蓋度是基于10 m分辨率的哨兵二號(hào)影像和NDVI經(jīng)驗(yàn)公式計(jì)算得到的。哨兵二號(hào)有兩顆極軌衛(wèi)星，在赤道地區(qū)的時(shí)間重訪周期是5 d，在東北高緯度地區(qū)的重放周期小于5 d。在文中，從全年哨兵二號(hào)數(shù)據(jù)中挑選出了無云或少云的高質(zhì)量影像，以去除不利天氣等的影響，以對(duì)地表植被實(shí)現(xiàn)連續(xù)的有效觀測(cè)。需要注意的是，依據(jù)10 m分辨率的哨兵數(shù)據(jù)計(jì)算得到的植被覆蓋度可能仍與實(shí)地測(cè)量的植被覆蓋度存在一定的區(qū)別[21]。為了便于確定時(shí)間窗口，采用的植被覆蓋度為0.6的閾值，而實(shí)際閾值的確定需要結(jié)合野外的實(shí)地工作，這也是下一步計(jì)劃開展的一個(gè)研究?jī)?nèi)容。

在使用經(jīng)驗(yàn)公式從NDVI轉(zhuǎn)化為FVC時(shí)采用的是每個(gè)NDVI時(shí)間序列最大、最小值開展的計(jì)算。之所以采用這種方法，是因?yàn)镹DVI容易飽和[22]，給定地區(qū)一年內(nèi)NDVI的最大值可以當(dāng)作植被蓋度100%的情況，而最小的NDVI當(dāng)作該點(diǎn)植被蓋度為0的情況。因?yàn)?，針?duì)每個(gè)研究區(qū)不知道其真實(shí)的地表覆蓋情況，所以難以針對(duì)每景影像確定其最小和最大植被蓋度的NDVI。這一點(diǎn)也可能會(huì)給結(jié)果帶來誤差。盡管如此，植被蓋度應(yīng)該和NDVI在一年內(nèi)的時(shí)間序列都有類似的變化規(guī)律，因此對(duì)最終的結(jié)果影響應(yīng)該比較有限。

5 結(jié)論

1)寒溫帶開展無人機(jī)外業(yè)的窗口期為4月中旬—5月初；

2)中溫帶的飛行窗口期較寒溫帶長，在4月中旬—5月中旬；

3)暖溫帶可選擇的窗口期范圍較寒溫帶和中溫帶大，時(shí)間范圍分別為3月中旬—4月底，9月底—11月初。