肖 武 陳佳樂(lè) 笪宏志 任 河 張建勇 張 雷

(1.中國(guó)礦業(yè)大學(xué)(北京)土地復(fù)墾與生態(tài)重建研究所， 北京 100083； 2.浙江大學(xué)公共管理學(xué)院， 杭州 310058)

0 引言

煤炭資源的開(kāi)發(fā)利用為現(xiàn)代化建設(shè)做出了巨大貢獻(xiàn)，但長(zhǎng)時(shí)間、大規(guī)模、高強(qiáng)度的開(kāi)采不可避免地破壞原生的礦床地質(zhì)條件和占用、破壞大量的土地[1]，造成原有生態(tài)系統(tǒng)的嚴(yán)重破壞，引發(fā)植被退化、水土流失、物種減少等一系列生態(tài)效應(yīng)[2]。如何對(duì)受損土地進(jìn)行快速監(jiān)測(cè)、確定受損范圍與程度一直以來(lái)都是研究的重點(diǎn)。

國(guó)內(nèi)外對(duì)采礦區(qū)的生態(tài)損毀監(jiān)測(cè)已有不少研究，但傳統(tǒng)方法大多集中于利用衛(wèi)星遙感技術(shù)對(duì)礦區(qū)土地利用類(lèi)型或積水面積變化進(jìn)行監(jiān)測(cè)，如VENKATARAMAN等利用衛(wèi)星遙感技術(shù)對(duì)礦區(qū)進(jìn)行了監(jiān)測(cè)，結(jié)果表明采礦活動(dòng)的增加與研究區(qū)植被和土地退化程度有關(guān)；CHRISTAIN等為了估算動(dòng)態(tài)系統(tǒng)“沉降-地下水-植被”的變化，進(jìn)行了大面積的時(shí)序性衛(wèi)星觀(guān)測(cè)，得到土地利用和土地覆蓋的光譜特征可能與水文狀況變化引起的特定植被現(xiàn)象有關(guān)的結(jié)論；吳雪茜等應(yīng)用開(kāi)采沉陷預(yù)測(cè)技術(shù)、地理信息及衛(wèi)星遙感技術(shù)，對(duì)淮南礦區(qū)土地、水域演變趨勢(shì)進(jìn)行研究并提出治理對(duì)策；肖武等利用遙感影像解譯了淮南張集礦區(qū)15年內(nèi)積水面積變化，結(jié)合高潛水位礦區(qū)的特點(diǎn)，論證了采煤沉陷地構(gòu)建平原水庫(kù)的可行性。在監(jiān)測(cè)礦區(qū)生態(tài)擾動(dòng)影響中，遙感衛(wèi)星尺度大、分辨率低、時(shí)間周期長(zhǎng)且遙感指標(biāo)宏觀(guān)化、單一化，無(wú)法對(duì)處于開(kāi)采過(guò)程中的礦區(qū)進(jìn)行實(shí)時(shí)、動(dòng)態(tài)、有效的監(jiān)測(cè)。近年來(lái)，隨著無(wú)人機(jī)遙感技術(shù)的不斷成熟，加之無(wú)人機(jī)響應(yīng)快、周期短、精度高、易操作、成本低的特點(diǎn)，無(wú)人機(jī)遙感技術(shù)在作物生態(tài)參數(shù)反演估算、農(nóng)田生態(tài)環(huán)境信息監(jiān)測(cè)等方面得到廣泛應(yīng)用[7-9]，但在礦區(qū)的應(yīng)用尚處于起步階段，大多數(shù)研究都集中于對(duì)礦區(qū)各類(lèi)采礦設(shè)施與土地利用情況分類(lèi)與監(jiān)測(cè)、礦區(qū)地?zé)豳Y源分布調(diào)查、非法與越界開(kāi)采識(shí)別、露天礦工程量(采剝量、堆放量等)計(jì)算等幾個(gè)方面，而對(duì)于礦區(qū)生態(tài)環(huán)境擾動(dòng)識(shí)別、植被長(zhǎng)勢(shì)監(jiān)測(cè)的應(yīng)用還有待于進(jìn)一步挖掘。

采煤擾動(dòng)對(duì)礦區(qū)生態(tài)的影響很大程度上表征在作物生化參數(shù)的改變上，其中，作物生物量指單位面積內(nèi)作物累積有機(jī)物質(zhì)的總量，其與長(zhǎng)勢(shì)和單產(chǎn)密切相關(guān)，是農(nóng)業(yè)監(jiān)測(cè)中應(yīng)用最為廣泛的指標(biāo)之一。玉米生物量直接影響玉米整體的長(zhǎng)勢(shì)，進(jìn)而決定玉米的產(chǎn)量，而在高潛水位采煤沉陷影響下，區(qū)內(nèi)地形受到沉陷的影響導(dǎo)致潛水位相對(duì)升高，產(chǎn)生漬害進(jìn)而影響農(nóng)作物的正常生長(zhǎng)，其直觀(guān)反映為農(nóng)作物生物量降低，作物枯死和減產(chǎn)、絕產(chǎn)。因此實(shí)時(shí)、有效地對(duì)采煤沉陷區(qū)作物的長(zhǎng)勢(shì)監(jiān)測(cè)是礦區(qū)土地復(fù)墾、生態(tài)修復(fù)的前提。同時(shí)作物生化參數(shù)與植被指數(shù)之間也存在著明顯的相關(guān)性，利用植被指數(shù)等遙感參量反演作物生物量是監(jiān)測(cè)作物長(zhǎng)勢(shì)的重要方法[11-12]，以往的研究大多用于精準(zhǔn)農(nóng)業(yè)相關(guān)研究，裴浩杰等基于多種生化參數(shù)指標(biāo)與光譜指數(shù)結(jié)合偏最小二乘回歸構(gòu)建模型判斷出研究區(qū)小麥整體的長(zhǎng)勢(shì)差異；楊粉團(tuán)等針對(duì)玉米粘蟲(chóng)災(zāi)情構(gòu)建了基于重歸一化植被指數(shù)多時(shí)相的葉片生物量定量模型，實(shí)現(xiàn)了對(duì)玉米粘蟲(chóng)災(zāi)情程度的有效監(jiān)測(cè)。需要指出的是，現(xiàn)階段研究主要針對(duì)正常農(nóng)田或受某單一脅迫影響植被的生化參數(shù)反演而展開(kāi)的，并且所用傳感器多為普通數(shù)碼相機(jī)或近紅外多光譜相機(jī)。KROSS等利用Rapid Eye數(shù)據(jù)分析了NDVI、 red-edge NDVI (NDVIre)和Modified Triangular Vegetation Index 2(MTVI2)等7種寬波段指數(shù)與大豆和玉米生物量的相關(guān)性，得出NDVIre與生物量的相關(guān)性最高(R2=0.78)。JIN等則通過(guò)分析發(fā)現(xiàn)，MTVI2相較NDVI和EVI等與生物量的相關(guān)性更強(qiáng)，決定性系數(shù)分別為0.72和0.7，均方根誤差為198.65、227.41 g/m2，可以證明紅邊波段對(duì)作物生長(zhǎng)參數(shù)的響應(yīng)更為顯著。因此，本文嘗試?yán)脽o(wú)人機(jī)搭載含有紅邊波段的多光譜相機(jī)，結(jié)合地面實(shí)測(cè)樣本數(shù)據(jù)，研究以植被指數(shù)反演作物生物量表征耕地?fù)p毀程度的可行性。根據(jù)已有研究，篩選相關(guān)性較好的植被指數(shù)并借助紅邊波段對(duì)作物長(zhǎng)勢(shì)響應(yīng)敏感的優(yōu)勢(shì)加以擴(kuò)展，構(gòu)建適用于高潛水位采煤沉陷地的玉米生物量遙感監(jiān)測(cè)模型，以期為礦區(qū)土地?fù)p毀監(jiān)測(cè)與評(píng)價(jià)、土地復(fù)墾與生態(tài)修復(fù)等提供基礎(chǔ)數(shù)據(jù)與理論支撐。

1 材料與方法

1.1 研究區(qū)概況

東灘煤礦位于山東省濟(jì)寧市境內(nèi)，跨兗州、鄒城、曲阜三市(縣)，地理位置116°50′49″～116°56′56″E，35°24′11″～35°31′25″N。地處魯中低山丘陵到平原洼地的過(guò)渡地帶，整體地勢(shì)由東北向西南逐漸降低，潛水埋深為2 m左右，屬于高潛水位礦區(qū)。礦區(qū)內(nèi)土壤類(lèi)型多為褐土，土質(zhì)較好，土壤肥沃，耕性良好，主要實(shí)行冬小麥與夏玉米的輪作模式，其中，夏玉米一般在當(dāng)年6月種植，10月收獲，是重要的糧食產(chǎn)區(qū)。礦區(qū)內(nèi)主要河流有白馬河與泥河，向南流入南陽(yáng)湖，均為季節(jié)性河流。研究區(qū)為東灘煤礦A1工作面開(kāi)采影響的范圍(圖1)。其中A1大采高綜放工作面位于三采區(qū)南部3號(hào)煤層，地面平均標(biāo)高+50.72 m，開(kāi)采平均標(biāo)高-497.8 m。所采煤層煤厚平均8.30 m，煤層傾角平均為4°，自2014年8月開(kāi)采結(jié)束至實(shí)驗(yàn)日期2017年8月，研究區(qū)已基本達(dá)到穩(wěn)沉狀態(tài)，地表形成了大面積的沉陷區(qū)，加之潛水位較高，地面沉陷深度最深已達(dá)5 m，大量的優(yōu)質(zhì)耕地沉入水中，形成典型的礦區(qū)積水盆地。研究區(qū)內(nèi)沉陷程度自積水中心向四周逐漸減輕，參考地表實(shí)測(cè)數(shù)據(jù)沉降范圍基本分布在0.15～5.25 m。地表主要地物分布類(lèi)型由內(nèi)及外依次為蘆葦、沼澤、灘涂、水淹的農(nóng)作物。結(jié)合該工作面的開(kāi)采特征，依據(jù)采煤沉陷預(yù)計(jì)軟件獲得沉陷區(qū)的下沉等值線(xiàn)，按照傳統(tǒng)定義設(shè)置下沉10 mm等值線(xiàn)為影響邊界。采用走向、傾向、角平分線(xiàn)方向布設(shè)3條實(shí)地采樣線(xiàn)實(shí)地采樣線(xiàn)L1、L2、L3，長(zhǎng)度分別為550、540、620 m，其中近水區(qū)域植被長(zhǎng)勢(shì)受采煤擾動(dòng)影響明顯，為了采集更細(xì)致的樣本信息來(lái)表征植被的生長(zhǎng)變化趨勢(shì)，3條樣線(xiàn)上樣點(diǎn)布設(shè)間距均隨下沉值的減小而增大，依次為5、10、20、30、60 m，記錄所有樣點(diǎn)的空間位置和相應(yīng)環(huán)境信息，共計(jì)54個(gè)采樣點(diǎn)，3條樣線(xiàn)L1、L2、L3最遠(yuǎn)端點(diǎn)布設(shè)在擾動(dòng)邊界之外作為參照點(diǎn)(圖1)。本次無(wú)人機(jī)觀(guān)測(cè)試驗(yàn)的時(shí)間為2017年8月15日(夏玉米灌漿期)，灌漿期是玉米關(guān)鍵生育期中生物量應(yīng)達(dá)到最大和最適宜的時(shí)期，與后期產(chǎn)量有著顯著相關(guān)性，是利用生物量監(jiān)測(cè)玉米長(zhǎng)勢(shì)的最佳時(shí)期。

1.2 數(shù)據(jù)獲取與處理方法

1.2.1多光譜影像數(shù)據(jù)獲取

使用的遙感傳感器為瑞士parrotsequoia無(wú)人機(jī)多光譜相機(jī)(下文簡(jiǎn)稱(chēng)sequoia)，是專(zhuān)為無(wú)人機(jī)在農(nóng)業(yè)科研、調(diào)查而研發(fā)的，能適用多種飛行器。它不僅可以獲取1 600萬(wàn)像素RGB三原色照片，還能獲取120萬(wàn)像素的綠光(波長(zhǎng)550 nm、帶寬40 nm)、紅光(波長(zhǎng)660 nm、帶寬40 nm)、紅邊(波長(zhǎng)735 nm、帶寬10 nm)、近紅外(波長(zhǎng)790 nm、帶寬40 nm)等波段影像(圖2)。依靠自帶的光照傳感器可記錄光照條件并自動(dòng)校準(zhǔn)4個(gè)多光譜傳感器的獨(dú)立亮度，同時(shí)內(nèi)置GPS和IMU。搭載sequoia的遙感平臺(tái)為四旋翼無(wú)人機(jī)大疆M100。試驗(yàn)現(xiàn)場(chǎng)如圖3所示，此時(shí)天氣晴朗無(wú)風(fēng)，視野良好。無(wú)人機(jī)飛行高度110 m，設(shè)定航速9 m/s。傳感器鏡頭視場(chǎng)角15°，鏡頭垂直向下，地面分辨率13 cm，航拍面積1.1 km2，共獲取4 980幅影像，覆蓋整個(gè)研究區(qū)。

圖1 研究區(qū)概況和采樣點(diǎn)布設(shè) Fig.1 Map of study area and layout of sampling points

在地面鋪設(shè)好傳感器自帶的光譜反射率校正板，每個(gè)架次起飛前先手持飛機(jī)在校正板正上方1.5 m處拍照，獲得當(dāng)時(shí)條件下的標(biāo)準(zhǔn)反射率，選用軟件pix4dmapper對(duì)影像進(jìn)行處理，在處理過(guò)程中利用標(biāo)準(zhǔn)反射率校準(zhǔn)所有的航拍影像以得到理想的處理成果，并使用ENVI5.1軟件，以研究區(qū)數(shù)碼正射影像為參考影像在圖像不同位置均勻選取30個(gè)參考點(diǎn)對(duì)多光譜影像進(jìn)行幾何精校正，經(jīng)檢驗(yàn)圖像幾何糾正誤差小于0.5個(gè)像元，根據(jù)生物量地面測(cè)量對(duì)應(yīng)的樣點(diǎn)位置構(gòu)建感興趣區(qū)域(Region of interest，ROI)，以ROI范圍內(nèi)地物的平均反射率光譜值作為該樣點(diǎn)玉米葉片反射率光譜，得到各樣點(diǎn)的反射率光譜數(shù)據(jù)。

圖2 多光譜相機(jī) Fig.2 Multispectral camera

圖3 航測(cè)現(xiàn)場(chǎng) Fig.3 Scene of aerial survey site

1.2.2樣點(diǎn)地表高程獲取

為了獲取沉陷區(qū)下沉情況，在樣點(diǎn)信息采集過(guò)程中，同時(shí)采用經(jīng)礦區(qū)已知點(diǎn)校正過(guò)的南方銀河一號(hào)RTK沿3條樣線(xiàn)L1、L2、L3測(cè)得樣點(diǎn)地表高程。

1.2.3生物量采集

在進(jìn)行航測(cè)的同時(shí)，在地面樣線(xiàn)上同步取樣。如圖1所示分別沿走向、傾向、角平分線(xiàn)方向樣線(xiàn)由內(nèi)及外進(jìn)行采樣。因地面沉陷，部分樣點(diǎn)處于常年積水區(qū)域和濕地中，植被樣本多為蘆葦、雜草。本文研究對(duì)象為耕地，故最后數(shù)據(jù)采集為耕地上39個(gè)玉米樣方的地上干生物量(Above ground biomass, AGB)。為不影響品種正常選育，選取1 m×1 m樣方范圍內(nèi)邊緣行取2株長(zhǎng)勢(shì)均勻的玉米作為試驗(yàn)樣本，并測(cè)量樣方內(nèi)的總株數(shù)、行間距、株高、基徑參考信息。用實(shí)驗(yàn)室高精度天平測(cè)量樣本玉米鮮生物量，平均后得到樣方內(nèi)玉米鮮生物量，在實(shí)驗(yàn)室內(nèi)將其先在105℃高溫條件下殺青1 h，接著在80℃恒溫條件下干燥24 h至恒質(zhì)量，稱(chēng)其質(zhì)量，并根據(jù)實(shí)測(cè)總株數(shù)計(jì)算出樣方內(nèi)的總干生物量，結(jié)合樣方的總面積，計(jì)算出每平方米的地上干生物量(g/m2)[18]。

1.3 植被指數(shù)選取

圖4 健康與受脅迫植被光譜反射曲線(xiàn) Fig.4 Spectral reflectance curves of healthy and stressed vegetations

基于植被光合色素與光譜反射率特征波段具有很強(qiáng)的相關(guān)性，特征波段的選取是需要參考色素的光譜特征信息[19]。如圖4所示,在可見(jiàn)光范圍內(nèi)，葉綠素出現(xiàn)2個(gè)強(qiáng)烈的吸收峰、分別為以450 nm和640～680 nm為中心的藍(lán)波段和紅波段，吸收的峰值出現(xiàn)在670 nm；一個(gè)強(qiáng)烈的反射峰出現(xiàn)在550 nm左右的綠光波段，此時(shí)葉綠素吸收系數(shù)最小。在700 nm附近，也為對(duì)應(yīng)的葉綠素吸收谷值波段，因此550 nm與700 nm附近常常被選作抗干擾的特征波段，用來(lái)削弱非光合作用物質(zhì)引起的光合有效輻射；紅邊區(qū)間700～780 nm是葉綠素在紅邊波段的強(qiáng)吸收到近紅波段多次散射形成的高反射平臺(tái)的爬升嵴，這個(gè)區(qū)間能非常敏感地響應(yīng)植被營(yíng)養(yǎng)、長(zhǎng)勢(shì)、水分。

生化參數(shù)等指標(biāo)的變化已被國(guó)內(nèi)外的大量研究結(jié)果證實(shí)。其中，當(dāng)植被長(zhǎng)勢(shì)良好時(shí)，紅邊會(huì)向長(zhǎng)波方向移動(dòng)(紅移)，當(dāng)植被遭受蟲(chóng)害、污染、水分等因素脅迫時(shí)，紅邊則會(huì)向短波方向移動(dòng)(藍(lán)移)[20]；在780 nm以后，近紅外波段葉綠素對(duì)電磁波的吸收特征微弱，也常常將780～800 nm附近波段選作特征波段。

結(jié)合以上對(duì)光譜曲線(xiàn)的特征波段的分析與本次實(shí)驗(yàn)多光譜傳感器的多通道優(yōu)勢(shì)，初步選取了現(xiàn)有研究中常見(jiàn)的且符合本次實(shí)驗(yàn)特征波段組合的植被指數(shù)，并在此基礎(chǔ)上，利用紅邊波段對(duì)植物葉綠素響應(yīng)敏感的優(yōu)勢(shì)[21]，用紅邊波段替換紅或綠波段來(lái)擴(kuò)展現(xiàn)有植被指數(shù)，共選取植被指數(shù)22個(gè)。表1列出了部分指數(shù)的計(jì)算公式及出處。

表1 部分選用的反演玉米生物量的植被指數(shù)公式 Tab.1 Vegetation index formula inverted from biomass

注：ρNIR、ρR、ρG、ρRed-edge分別表示近紅外、紅光、綠光、紅邊波段的反射率。

與現(xiàn)有研究中針對(duì)某單一地類(lèi)或受某單一脅迫影響的植被生物量反演不同的是，該研究區(qū)為礦區(qū)采煤沉陷耕地，地下采煤擾動(dòng)大，作物生長(zhǎng)情況受脅迫因子較復(fù)雜，選用的22種植被指數(shù)與玉米生物量的相關(guān)性表現(xiàn)可能與前人研究結(jié)果不一致。因此為保證獲得最優(yōu)的反演精度，提前計(jì)算22種植被指數(shù)與實(shí)測(cè)生物量的相關(guān)性并根據(jù)其高低進(jìn)行篩選，選取其中相關(guān)性較好的部分植被指數(shù)作為參照。在多種新的植被指數(shù)都與生物量的相關(guān)性良好時(shí)，若忽視植被指數(shù)間的近似線(xiàn)性關(guān)系則可能會(huì)使得回歸方程不穩(wěn)定，有些植被指數(shù)對(duì)生物量影響的顯著性被隱蔽起來(lái)，某些回歸系數(shù)的符號(hào)與實(shí)際意義不相符合等[22]。改進(jìn)前后植被指數(shù)兩兩之間的方差膨脹因子(Variance inflation factor，VIF)，以分析和控制變量間的多重共線(xiàn)性；再按照相關(guān)性高且顯著的原則，并設(shè)定變量最大VIF小于10，確定最優(yōu)植被指數(shù)及其組合；最后分別將最優(yōu)植被指數(shù)及其組合作為入選變量輸入多元逐步回歸模型、BP神經(jīng)網(wǎng)絡(luò)模型進(jìn)行擬合。

VIF=1/(1-R2)

(1)

式中VIF——植被指數(shù)間的方差膨脹因子

R——植被指數(shù)間的相關(guān)系數(shù)

1.4 數(shù)據(jù)分析方法

在野外實(shí)驗(yàn)中，由于實(shí)驗(yàn)環(huán)境和儀器本身等客觀(guān)因素或偶然因素的影響，可能會(huì)導(dǎo)致獲取的數(shù)據(jù)中存在異常數(shù)據(jù)。為了剔除可能存在的異常數(shù)據(jù)，提高數(shù)據(jù)的可靠性及分析結(jié)果的精度，在建立回歸方程前，對(duì)樣本數(shù)據(jù)進(jìn)行回歸分析與殘差分析[23]。根據(jù)求出的殘差和標(biāo)準(zhǔn)殘差值選擇樣本，最終選取了殘差最小的36個(gè)樣本。表2為剔除粗差后的實(shí)測(cè)樣本的統(tǒng)計(jì)特征。

表2 優(yōu)化后玉米地上干生物量統(tǒng)計(jì)特征 Tab.2 Statistical characteristics of AGB in optimized maize g/m2

在確定好樣本集后，以預(yù)處理后的無(wú)人機(jī)多光譜圖像為數(shù)據(jù)源，基于0.05 m空間分辨率數(shù)碼正射影像精確劃分的36個(gè)實(shí)地樣方范圍，提取各個(gè)樣方的植被指數(shù)，在建立經(jīng)驗(yàn)?zāi)Ｐ椭?，先分別討論初步選取的植被指數(shù)與玉米生物量的相關(guān)性，根據(jù)結(jié)果進(jìn)行初步排序，篩選出相關(guān)性較好的植被指數(shù)并加以改進(jìn)。利用SPSS 22.0軟件對(duì)灌漿期的36個(gè)玉米生物量實(shí)測(cè)值進(jìn)行隨機(jī)抽樣，選出26個(gè)作為建模樣本,10個(gè)作為檢驗(yàn)樣本。針對(duì)改進(jìn)前后的植被指數(shù)利用一元回歸、基于最小二乘法的多元逐步線(xiàn)性回歸、BP神經(jīng)網(wǎng)絡(luò)法進(jìn)行擬合。從中篩選出與玉米生物量相關(guān)性最好的模型。

1.5 模型評(píng)價(jià)方法

選取上述構(gòu)建的各種模型決定系數(shù)、均方根誤差和估測(cè)精度進(jìn)行模型分析檢驗(yàn),其計(jì)算公式為

(2)

(3)

(4)

圖5 采樣線(xiàn)上玉米生物量與地表高程隨距離變化趨勢(shì) Fig.5 Changing trends of maize biomass and land surface elevation with distance in each sampling line of corn

2 結(jié)果與分析

2.1 采煤擾動(dòng)與遙感參數(shù)關(guān)系

采煤沉陷造成采空區(qū)上覆巖層的拉伸、斷裂與彎曲等移動(dòng)變形，致使高潛水位的地區(qū)地下水位相對(duì)上升，沉陷中心區(qū)域耕地變?yōu)檎訚?，玉米在?nèi)的農(nóng)作物無(wú)法繼續(xù)種植，周邊影響區(qū)域地表變形減弱改變了土壤持水能力和通氣狀況，影響有機(jī)物和礦物質(zhì)的分解、淋溶和沉積，以及由此而引起的土壤侵蝕，使土壤保水能力變差，養(yǎng)分流失嚴(yán)重，土質(zhì)惡化，造成耕地上玉米不同程度的減產(chǎn)[25]。本文針對(duì)此類(lèi)現(xiàn)象，選取3條樣線(xiàn)上玉米生物量、地表高程和植被指數(shù)值隨樣點(diǎn)距積水中心距離的變化趨勢(shì)對(duì)比分析，可知受地下采煤擾動(dòng)的影響，樣點(diǎn)上地表高程不同程度的沉陷，導(dǎo)致沉陷耕地土壤質(zhì)量不同程度的受損，傳播到地表，玉米生物量呈現(xiàn)隨地表下沉值減少、與積水中心距離增加而增加至最大值后，略微下降后便穩(wěn)定在某一值的現(xiàn)象(圖5)，表征到樣線(xiàn)上植被指數(shù)值上，植被指數(shù)變化趨勢(shì)與玉米生物量一致，隨地表高程增加而增加至最大值后，略微變化后穩(wěn)定在某一值(圖6)。這說(shuō)明在采煤沉陷耕地的復(fù)雜背景中，地下的采煤擾動(dòng)導(dǎo)致地表作物長(zhǎng)勢(shì)差異進(jìn)而反映到遙感參數(shù)上是一個(gè)連續(xù)的過(guò)程，通過(guò)反演作物生物量來(lái)監(jiān)測(cè)耕地?fù)p毀程度是可行的。

圖6 采樣線(xiàn)上植被指數(shù)值與地表高程隨距離變化趨勢(shì) Fig.6 Changing trends of vegetation index and land surface elevation with distance in each sampling line of corn

2.2 優(yōu)選植被指數(shù)

2.2.1不同光譜植被指數(shù)對(duì)灌漿期玉米地上干物質(zhì)量的敏感性分析

通過(guò)對(duì)不同光譜植被指數(shù)與玉米地上干物質(zhì)量 進(jìn)行相關(guān)分析，結(jié)果發(fā)現(xiàn)22種植被指數(shù)與灌漿期玉米生物量均極顯著相關(guān)(表3)，但不同的植被指數(shù)與地上干物質(zhì)量之間的相關(guān)程度不同，即不同的植被指數(shù)對(duì)玉米地上干物質(zhì)量變化表現(xiàn)出不同的敏感性。其中，結(jié)合表4中各指數(shù)的計(jì)算方式，可以得到初步的規(guī)律：在現(xiàn)有的指數(shù)中近紅波段和紅邊的組合與玉米地上干物質(zhì)量的相關(guān)系數(shù)最高，可達(dá)到0.9以上，而近紅波段和綠、紅波段組合的指數(shù)與玉米地上干物質(zhì)量的相關(guān)性則相對(duì)差一些，相關(guān)系數(shù)在0.6～0.8之間。表明紅邊指數(shù)對(duì)灌漿期玉米地上干物質(zhì)量的敏感性較強(qiáng)，能很好地反映玉米AGB的變化情況。

2.2.2最優(yōu)估算模型變量選擇

表4列出基于22種植被指數(shù)與灌漿期玉米AGB擬合的一元回歸模型。從表4可以看出，利用各種植被指數(shù)擬合的一元曲線(xiàn)回歸模型(SCRM)均達(dá)到較好的擬合效果。在所選植被指數(shù)的最優(yōu)模型中，紅邊植被指數(shù)與玉米AGB擬合的指數(shù)模型效果較好。其中，在原有植被指數(shù)WDRVI、MVI的基礎(chǔ)上改進(jìn)的紅邊指數(shù)WDRVIre、MVIre的擬合效果最優(yōu)，其決定系數(shù)R2達(dá)0.80、0.81；紅邊指數(shù)CIredge、SRredge、MSRredge次之，R2為0.79，其他指數(shù)中擬合效果較好的是NDVI、GNDVI，R2也達(dá)0.68。

表3 不同植被指數(shù)與灌漿期玉米地上生物量的 相關(guān)系數(shù)(n=26) Tab.3 Correlation coefficient between different vegetation indices and corn aboveground biomass at filling stage(n=26)

注：** 表示在0.01水平上顯著相關(guān)。

綜合以上一元最優(yōu)模型的分析結(jié)果，為進(jìn)一步提高模型精度，多元回歸模型的入選參量的篩選應(yīng)在一元最優(yōu)模型精度較高的指數(shù)中進(jìn)行，故對(duì)22種一元模型決定系數(shù)在0.65以上的植被指數(shù)進(jìn)行雙變量相關(guān)性分析，根據(jù)相關(guān)系數(shù)計(jì)算方差膨脹因子，對(duì)選取的植被指數(shù)的多重共線(xiàn)性分析(表5)，由表分析可知，單變量中改進(jìn)后的紅邊指數(shù)與生物量的相關(guān)性最高，但是紅邊波段指數(shù)之間都存在著嚴(yán)重的共線(xiàn)性分析，因此，多變量構(gòu)建模型應(yīng)在除紅邊和近紅波段組合外，考慮加入綠波段、紅波段與近紅波段之間的組合指數(shù)以進(jìn)一步提高模型的精度。

根據(jù)表5的結(jié)果，為排除后期模型可能出現(xiàn)的過(guò)擬合現(xiàn)象且考慮模型計(jì)算的簡(jiǎn)潔性，選取在單元回歸模型精度最高的紅邊指數(shù)MVIre和兩兩之間共線(xiàn)性低的紅、綠波段指數(shù)NDVI、GNDVI作為多元變量入選多元回歸和BP神經(jīng)網(wǎng)絡(luò)反演模型。

2.3 生物量反演建模與驗(yàn)證

基于模型入選變量，分別采用多元回歸、BP神經(jīng)網(wǎng)絡(luò)方法建模并選取其中最佳的反演模型，其中單元最優(yōu)模型已在上文提及過(guò)(表4)；運(yùn)用SPSS 22.0軟件基于入選變量的所有組合可能進(jìn)行MLR分析與建模；采用Matlab R2012a軟件編程構(gòu)建多種輸入特征組合得BPNN模型，分別比較所得到的MLR、BPNN模型的決定系數(shù)和誤差，篩選出各自的最優(yōu)模型組合并基于交叉驗(yàn)證法[26]進(jìn)行模型的檢驗(yàn)和比較(表6)。

表4 各植被指數(shù)和玉米生物量的相關(guān)系數(shù)(n=26) Tab.4 Correlation coefficient between 22 vegetation indices and biomass(n=26)

表5 相關(guān)性較好的植被指數(shù)間的方差膨脹因子 Tab.5 Variance expansion factor among well-correlated vegetation indices

表6 生物量的估算與驗(yàn)證 Tab.6 Biomass estimation and verification

表6和圖7顯示了改進(jìn)后的植被指數(shù)及組合在不同回歸模型下對(duì)玉米生物量的最佳預(yù)測(cè)能力。在所有回歸模型中，由改進(jìn)后的紅邊MVI指數(shù)與綠波指數(shù)GNDVI構(gòu)建的BP神經(jīng)網(wǎng)絡(luò)模型建模精度最高，R2=0.83、RMSE為168.09 g/m2。為了驗(yàn)證模型的精度，利用剩下的樣本個(gè)數(shù)進(jìn)行驗(yàn)證，結(jié)果顯示驗(yàn)證模型R2和RMSE分別為0.64和178.72 g/m2，估測(cè)精度EA為79.4%，比其他模型R2增加了0.1～0.17、RMSE減少了29.65～60.23 g/m2、估測(cè)精度提高了3.3%～7.1%。其建模精度較為理想，由于上述選擇建模樣本和驗(yàn)證樣本具有不確定性，所以本文再次用交叉驗(yàn)證法來(lái)進(jìn)行生物量的估算，結(jié)果顯示R2和RMSE分別為0.78和177.23 g/m2，預(yù)測(cè)值和實(shí)測(cè)值的散點(diǎn)圖具有較好的一致性，該模型與生物量的相關(guān)性達(dá)到極顯著水平，表明利用上述模型顯著提高了采煤沉陷耕地玉米地上干生物量的反演精度。

圖7 各模型生物量預(yù)測(cè)值與實(shí)測(cè)值之間的關(guān)系 Fig.7 Relationship between predicted and measured values of biomass

2.4 反演模型的應(yīng)用

利用上述模型結(jié)合研究區(qū)玉米種植范圍，進(jìn)行采煤沉陷盆地玉米灌漿期地表干生物量實(shí)地反演和制圖(圖8)，獲得研究區(qū)玉米生物量介于52.25～1 575.72 g/m2，平均值為873.08 g/m2，標(biāo)準(zhǔn)差為294.89 g/m2，與研究樣本的描述性統(tǒng)計(jì)結(jié)果(表2)較為一致。利用自然間斷分類(lèi)法基于數(shù)據(jù)中固有的自然分組，對(duì)分類(lèi)間隔加以識(shí)別，在數(shù)據(jù)值差異相對(duì)較大的位置處設(shè)置邊界進(jìn)行分組，將研究區(qū)玉米地表干生物量分為5個(gè)等級(jí)，統(tǒng)計(jì)每個(gè)等級(jí)的像元數(shù)量及所占比例，見(jiàn)表7。

圖8 基于BP神經(jīng)網(wǎng)絡(luò)模型采煤沉陷耕地灌漿期 玉米地上生物量反演圖 Fig.8 Inversion map of aboveground biomass of corn at grain filling stage based on BP neural network model in mining subsidence

遙感識(shí)別結(jié)果表明，研究區(qū)內(nèi)玉米生物量主要集中于592～1 050 g/m2，屬于Ⅲ級(jí)(592～860 g/m2)和Ⅳ級(jí)(860～1 050 g/m2)的區(qū)域占74.4%，地表生物量低于592 g/m2的玉米面積占21.4%，研究區(qū)內(nèi)玉米生物量整體偏低，反映耕地整體受采煤沉陷擾動(dòng)較嚴(yán)重；玉米地表生物量空間分布基本表現(xiàn)為沿3條樣線(xiàn)自擾動(dòng)邊界向盆地中心逐漸增加，至近水區(qū)域達(dá)到最低值。其中擾動(dòng)邊界內(nèi)L1樣線(xiàn)遠(yuǎn)端、L3覆蓋區(qū)域玉米地表生物量分布在1 050 g/m2以上，表明此處玉米長(zhǎng)勢(shì)大多未受影響，中部區(qū)域耕地多屬于中度和輕度損毀，反演結(jié)果與實(shí)地調(diào)查樣本基本一致。擾動(dòng)邊界外部分耕地也呈現(xiàn)輕度擾動(dòng)的現(xiàn)象，究其原因?yàn)檠芯繀^(qū)內(nèi)玉米均為農(nóng)民自然播種，玉米品種、施肥等因素不一致導(dǎo)致了玉米長(zhǎng)勢(shì)差異。

表7 基于BP模型反演的研究區(qū)灌漿期玉米地表干 生物量分等 Tab.7 Statistical analysis of surface dry biomass of maize at grain filling stage based on BP model inversion

3 討論

研究區(qū)受采煤擾動(dòng)地表已基本達(dá)到穩(wěn)沉且受損耕地上的玉米屬于灌漿期，玉米發(fā)育成熟，玉米生物量可以很好地反映作物的產(chǎn)量情況，進(jìn)一步凸現(xiàn)耕地?fù)p毀的狀況。應(yīng)用灌漿期的遙感影像和地面實(shí)測(cè)數(shù)據(jù)能對(duì)作物最終的整體長(zhǎng)勢(shì)進(jìn)行反演。在以后的研究中，應(yīng)在玉米生長(zhǎng)的各個(gè)時(shí)期同步或準(zhǔn)同步進(jìn)行無(wú)人機(jī)多光譜影像和地面數(shù)據(jù)的采集進(jìn)一步提高災(zāi)害過(guò)程實(shí)時(shí)監(jiān)測(cè)的準(zhǔn)確性和時(shí)效性。同時(shí)還可以設(shè)置地面模擬試驗(yàn)進(jìn)行多光譜的采集和分析，加強(qiáng)基于葉片葉面積指數(shù)、葉綠素、作物產(chǎn)量模型的反演機(jī)理研究和精度提高。但是，通過(guò)研究也發(fā)現(xiàn)，由于播種條件、地面微地形及田間管理的差異，玉米生物量在未開(kāi)采損毀區(qū)域也呈現(xiàn)一定的波段，如何去確定基于作物生物量的開(kāi)采沉陷損毀閾值也是下一步研究的關(guān)鍵。

在野外實(shí)地采集樣本的過(guò)程中，發(fā)現(xiàn)研究區(qū)內(nèi)植被的分布呈多樣化，由沉陷中心的沉水植物到浮葉植物再到近水岸邊的挺水植物，最后是長(zhǎng)勢(shì)各異的玉米，此次試驗(yàn)只研究了現(xiàn)有耕地中各樣觀(guān)測(cè)線(xiàn)上玉米生物量的變化趨勢(shì)，導(dǎo)致可用的樣本較少，限制了反演精度的提高。因此如何將觀(guān)測(cè)線(xiàn)上所有植被的生長(zhǎng)參數(shù)如生物量、作物產(chǎn)量等指標(biāo)標(biāo)準(zhǔn)化、統(tǒng)一化，構(gòu)建多指標(biāo)、多植被、時(shí)序性的反演模型是下一步研究工作的重點(diǎn)。

4 結(jié)論

(1)在礦區(qū)采煤沉陷地復(fù)雜的生態(tài)環(huán)境中，紅邊波段能夠在一定程度上提高生物量估算精度。通過(guò)對(duì)傳統(tǒng)指數(shù)計(jì)算中引入紅邊波段進(jìn)行改造，植被指數(shù)與玉米生物量的相關(guān)系數(shù)提高0.06～0.27，篩選得到相關(guān)性最優(yōu)的紅邊指數(shù)MVIre，構(gòu)建的冪指數(shù)模型決定系數(shù)R2達(dá)到0.81，比其他指數(shù)提高0.11～0.15，并在此基礎(chǔ)上，排除多重共線(xiàn)性，構(gòu)建多元線(xiàn)性回歸模型、BP神經(jīng)網(wǎng)絡(luò)模型。其中以指數(shù)MVIre、GNDVI為基礎(chǔ)的BP神經(jīng)網(wǎng)絡(luò)模型決定系數(shù)R2為0.83、預(yù)測(cè)均方根誤差為178.72 g/m2、估測(cè)精度為79.4%，3個(gè)指標(biāo)均為模型中最優(yōu)。表明該方法可有效提高玉米生物量遙感反演模型的精度。

(2)基于最佳模型對(duì)東灘煤礦A1工作面地表采煤沉陷耕地進(jìn)行空間識(shí)別和分析，結(jié)果表明：遙感指標(biāo)識(shí)別獲得的研究區(qū)采煤沉陷耕地玉米地表生物量空間分布情況與地面樣點(diǎn)描述性統(tǒng)計(jì)結(jié)果一致，研究區(qū)內(nèi)生長(zhǎng)期74.4%的玉米生物量分布在592～1 050 g/m2之間，21.4%的玉米生物量集中分布在592 g/m2以下，4.2%的玉米生物量分布在1 050 g/m2以上。空間分布上玉米生物量呈現(xiàn)從沉陷盆地內(nèi)部向邊緣逐漸增加的趨勢(shì)，這說(shuō)明了沉陷深度與玉米生物量的強(qiáng)相關(guān)性。由于農(nóng)民田間管理的不同，導(dǎo)致各測(cè)線(xiàn)之間玉米AGB也存在顯著差異，因此無(wú)法通過(guò)某一確定閾值與區(qū)間來(lái)劃分損毀程度。該研究可為同類(lèi)型其他高潛水位礦區(qū)土地?fù)p毀監(jiān)測(cè)與評(píng)價(jià)、土地復(fù)墾與生態(tài)修復(fù)等提供基礎(chǔ)數(shù)據(jù)與理論支撐。