基于分層抽樣的水稻種植面積回歸估計(jì)

申克建　王飛　裴志遠(yuǎn)

摘要：隨著農(nóng)業(yè)統(tǒng)計(jì)需求和遙感技術(shù)的發(fā)展，中國農(nóng)業(yè)部即將使用中分辨率遙感數(shù)據(jù)開展中國主要農(nóng)作物的全覆蓋遙感調(diào)查，如何修正全覆蓋調(diào)查結(jié)果使其更接近農(nóng)作物面積真值是需要研究的問題。以水稻為例，在區(qū)域尺度上利用基于樣本無人機(jī)數(shù)據(jù)解譯得到的水稻種植面積修正基于中分辨率數(shù)據(jù)解譯得到的水稻種植面積，分別比較分層分別回歸估計(jì)和分層合并回歸估計(jì)，利用基于全覆蓋無人機(jī)數(shù)據(jù)解譯得到的水稻種植面積驗(yàn)證回歸估計(jì)值。在給定5%誤差和95%置信水平下，抽樣比為036%（最小樣本量公式換算）的抽樣結(jié)果表明：（1）2種方法的估計(jì)精度都大于95%；（2）合并比估計(jì)的精度好于分別比估計(jì)。

關(guān)鍵詞：水稻種植面積；中分辨率數(shù)據(jù)；無人機(jī)數(shù)據(jù)；分層抽樣；回歸估計(jì)

中圖分類號： S127文獻(xiàn)標(biāo)志碼：

文章編號：1002-1302（2017）16-0206-05

收稿日期：2016-10-11

基金項(xiàng)目：國家自然科學(xué)基金（編號：41301506）；國家重點(diǎn)研發(fā)計(jì)劃（編號：2016YFB0501505）。

作者簡介：申克建（1982—），男，河北廊坊人，博士，工程師，研究方向?yàn)檗r(nóng)業(yè)遙感。E-mail：ashenkejian@126com。

農(nóng)作物面積是各級管理部門進(jìn)行生產(chǎn)指導(dǎo)、農(nóng)業(yè)補(bǔ)貼、農(nóng)業(yè)保險(xiǎn)核查等決策的重要依據(jù)，也是農(nóng)作物產(chǎn)量預(yù)測和環(huán)境管理等研究的重要依據(jù)。這些管理和研究需要解決的重要問題是提高農(nóng)作物面積估算精度。遙感技術(shù)結(jié)合地理信息系統(tǒng)（GIS）和全球定位系統(tǒng)（GPS）等現(xiàn)代地球信息技術(shù)能提供及時(shí)、客觀的農(nóng)作物面積信息。目前，30 m格網(wǎng)尺度、20～30個(gè)類別的土地覆蓋類型監(jiān)測總體精度大約在80%～85%之間，高分辨數(shù)據(jù)的分類精度超過95%[2-3]。隨著地球資源衛(wèi)星對地觀測能力不斷提高，中分辨率數(shù)據(jù)（10～60 m）正逐漸可以保證大尺度全覆蓋，但是由于混合像元等問題不能滿足農(nóng)業(yè)統(tǒng)計(jì)等行業(yè)應(yīng)用的精度要求[4]。高分辨率數(shù)據(jù)（<5 m）可以得到符合要求的行業(yè)應(yīng)用精度，但它受幅寬和重訪周期的局限，無法滿足調(diào)查區(qū)域的影像全覆蓋，甚至不能提供滿足抽樣要求的影像[5]。在3S技術(shù)（指遙感、地理信息系統(tǒng)、全球定位系統(tǒng)）的支持下，遙感和抽樣調(diào)查相結(jié)合的空間抽樣調(diào)查方法得到了廣泛應(yīng)用，地面調(diào)查是獲取抽樣樣本的重要手段，但地面調(diào)查存在強(qiáng)度大、費(fèi)用高等不足。近些年來，無人機(jī)被廣泛應(yīng)用于民用領(lǐng)域，它具有機(jī)動(dòng)、靈活的特點(diǎn)，有降低地面調(diào)查強(qiáng)度和保障高分辨數(shù)據(jù)獲取的優(yōu)勢。Laliberte 等應(yīng)用無人機(jī)遙感數(shù)據(jù)進(jìn)行牧場監(jiān)測，證明無人機(jī)遙感數(shù)據(jù)可以補(bǔ)充或者替代部分地面調(diào)查數(shù)據(jù)[8]。Breckenridge研究證明，無人機(jī)可以在較短的時(shí)間內(nèi)采集大范圍影像，比地面調(diào)查的效率高很多，他預(yù)見無人機(jī)對牧草地監(jiān)測的影響，將像全球定位系統(tǒng)影響導(dǎo)航和野外數(shù)據(jù)采集一樣[9]。

目前遙感與抽樣相結(jié)合的空間抽樣方案設(shè)計(jì)被廣泛應(yīng)用[10-13]。例如，美國農(nóng)業(yè)部現(xiàn)行農(nóng)作物面積估算方法是將遙感全覆蓋結(jié)果和6月調(diào)查結(jié)果進(jìn)行回歸分析，回歸模型是Y=a+b×X。其中Y是估計(jì)面積；X是全覆蓋作物遙感分類面積（cropland data layer classified acres，簡稱CDL），6月調(diào)查數(shù)據(jù)是全美11 000個(gè)地塊的地面調(diào)查信息（約占國土面積 250%，耕地地塊約為161 km2，城市區(qū)域地塊約為 016 km2，未利用地塊為 644～1288 km2）；a和b是基于6月調(diào)查數(shù)據(jù)和CDL樣本數(shù)據(jù)使用最小二乘法得到的[14]。孫佩軍等使用分層分別回歸估計(jì)估算河南省冬小麥面積，使用無人機(jī)解譯的冬小麥面積樣本進(jìn)行事后分層，輔助信息來自全覆蓋中分辨率數(shù)據(jù)（Landsat TM5和HJ）解譯的冬小麥面積[15]。

隨著農(nóng)業(yè)統(tǒng)計(jì)需求和遙感技術(shù)的發(fā)展，中國農(nóng)業(yè)部即將使用中分辨率遙感數(shù)據(jù)開展中國主要農(nóng)作物的全覆蓋遙感調(diào)查，如何修正全覆蓋調(diào)查結(jié)果使其更接近農(nóng)作物面積的真實(shí)值是需要研究的問題。本研究以水稻作物為例，研究在區(qū)域尺度上利用基于樣本無人機(jī)數(shù)據(jù)解譯得到的水稻種植面積修正基于中分辨率數(shù)據(jù)解譯得到的水稻種植面積，旨在為修正中分辨率全覆蓋遙感調(diào)查結(jié)果提供參考抽樣方法。

1材料與方法

11研究數(shù)據(jù)

研究數(shù)據(jù)位于江蘇省鹽城市，選取理由：（1）該區(qū)域位于水稻主產(chǎn)區(qū)；（2）在水稻生長季節(jié)內(nèi)獲取了無云覆蓋的光學(xué)遙感數(shù)據(jù)SPOT5；（3）該區(qū)域不是無人機(jī)禁飛區(qū)。

研究數(shù)據(jù)主要包括無人機(jī)數(shù)據(jù)和SPOT5多光譜數(shù)據(jù)。無人機(jī)數(shù)據(jù)拍攝于2014年9月24日至2014年9月27日，經(jīng)過處理后覆蓋面積12 16600 hm2，一共5條樣帶，每條樣帶長度約300 km，寬度07～10 km，通過人工目視解譯得到水稻種植面積5 23774 hm2，基于樣本無人機(jī)數(shù)據(jù)解譯得到的水稻解譯結(jié)果簡稱UAV-Rice，SPOT5多光譜數(shù)據(jù)獲取于2014年8月3日，該數(shù)據(jù)經(jīng)過裁剪后覆蓋和無人機(jī)相同的范圍，通過人工目視解譯得到水稻種植面積6 00504 hm2，基于SPOT5數(shù)據(jù)解譯得到的水稻解譯結(jié)果簡稱SPOT5-Rice，數(shù)據(jù)示意見圖1，數(shù)據(jù)用途說明見表1。

12研究方法

121總體思路

為了用樣本UAV-Rice修正SPOT5-Rice，本研究首先基于SPOT5-Rice構(gòu)建抽樣框，并按水稻種植面積規(guī)模分為6層；然后按回歸估計(jì)最小樣本量公式計(jì)算最小樣本量，按等比例分配將樣本分配到各層；隨后用分別比估計(jì)和聯(lián)合比估計(jì)進(jìn)行回歸估計(jì)，最后用UAV-Rice進(jìn)行回歸估計(jì)精度評估（圖2）。

122抽樣框準(zhǔn)備

基于SPOT5-Rice構(gòu)建100 m×100 m格網(wǎng)，將格網(wǎng)和SPOT5-Rice進(jìn)行空間分析，統(tǒng)計(jì)每個(gè)格網(wǎng)內(nèi)的SPOT5-Rice種植面積。

具體是通過ArcGIS Identity工具， 用SPOT5-Rice矢量

數(shù)據(jù)疊合成100 m×100 m矢量格網(wǎng)，這樣對100 m×100 m內(nèi)的每個(gè)地塊都給定了格網(wǎng)編號；最后用Dissolve工具以格網(wǎng)編號為依據(jù)，使每個(gè)100 m×100 m中的水稻合并成一個(gè)整體；最后統(tǒng)計(jì)每個(gè)100 m×100 m中水稻的種植面積，完成抽樣框構(gòu)建，N=8 627，這樣每個(gè)抽樣單元都有來自SPOT5-Rice的種植面積。同理將UAV-Rice矢量數(shù)據(jù)疊合抽樣框后的數(shù)據(jù)再融合，這樣每個(gè)抽樣單元都有來自UAV-Rice的種植面積。

最后，抽樣框中的每個(gè)抽樣單元屬性表包含抽樣單元編碼、來自SPOT5-Rice的種植面積和來自UAV-Rice的種植面積，如圖1所示。

123抽樣框分層標(biāo)志與層數(shù)的確定

本研究選擇以抽樣單元中水稻種植面積為分層標(biāo)志。一般分層數(shù)越多，抽樣統(tǒng)計(jì)越準(zhǔn)確，但抽樣方差的降低是與分層數(shù)的平方成反比的，分層抽樣理論要求每層至少必須抽取2個(gè)樣本單元，層數(shù)不超過總樣本量的一半，但在層數(shù)大于6層時(shí)方差的減少幅度將大為減緩，理論與實(shí)踐研究表明，層數(shù)以不超過6層為宜[16]。本試驗(yàn)將層數(shù)定為6層，分層界限確定采用累計(jì)頻數(shù)直方圖法。

124抽樣樣本量設(shè)計(jì)

在給定研究變量均值允許誤差（本研究設(shè)定為5%）和置信水平1-α（本研究設(shè)定為α=005）下，按最小樣本量[17]計(jì)算公式：

[JZ（]n=[SX（]n01+[SX（]n0N[SX）][SX）]；n0=[SX（]u2α/2·S2y·（1-ρ2）Δ2[SX）]。[JZ）][JY]（1）

式中：n為無放回抽樣的樣本量；n0為有放回抽樣得到的樣本量；N為總體單元數(shù)；uα/2為正態(tài)分布的上側(cè)分位數(shù)；Δ為研究變量均值允許誤差；S2y=[SX（]1N-1[SX）]·∑[DD（]Ni=1[DD）]（Yi-Y[TX-]）2，為研究變量（無人機(jī)水稻種植面積）總體方差；ρ=[SX（]SxySx·Sy[SX）]，為總體相關(guān)系數(shù)，其中輔助變量（中分辨率SPOT5水稻種植面積）總體標(biāo)準(zhǔn)差為Sx=[KF（][SX（]1N-1[SX）]·∑[DD（]Ni=1[DD）]（Xi-X[TX-]）2[KF）]，總體協(xié)方差為Sxy=[SX（]1N-1[SX）]·∑[DD（]Ni=1[DD）]（Xi-X[TX-]）·（Yi-Y[TX-]），其中X[TX-]為輔助變量總體均值，Y[TX-]為研究變量總體均值。

125分層隨機(jī)抽樣下的回歸估計(jì)與精度評估

分層回歸估計(jì)按分別回歸估計(jì)和合并回歸估計(jì)2種方法分別計(jì)算[17]，2種方法各層樣本量都按等比例分配，2種方法都按相同精度評估方法評估。

1251分別回歸估計(jì)

設(shè)第i層的目標(biāo)變量和輔助變量分別為yi和xi，其均值分別為y[TX-5]i和x[TX-5]i，層權(quán)重為Wi=Ni/N，輔助變量層總體為Ni，每一層的回歸系數(shù)為bi，輔助變量層樣本方差為s2xi=[SX（]1n-1[SX）]·∑[DD（]ni=1[DD）]（xi-x[TX-5]i）2，sxiyi=[SX（]1n-1[SX）]·∑[DD（]ni=1[DD）]（xi-x[TX-5]i）·（yi-y[TX-5]i）為層樣本協(xié)方差。

各層均值的回歸估計(jì)量：

式中：r為估計(jì)誤差；Y[DD（-1][HT6]^[DD）]為分別比估計(jì)或合并比估計(jì)的估計(jì)值；Y為真值（UAV-Rice）。

2結(jié)果與分析

在給定均值允許誤差5%和95%置信水平下，根據(jù)公式（1）計(jì)算得到n為31，將31個(gè)樣本按等比例分配到6層中，按2種回歸估計(jì)的公式計(jì)算得到表2所示結(jié)果，詳細(xì)計(jì)算過程見表3。

從表2可以看出，|025%|<|-103%|，|-333%|<|-457%|，說明合并比估計(jì)結(jié)果都好于分別比估計(jì)。同時(shí)2種方法估計(jì)精度都大于95%。

SPOT5-Rice解譯水稻存在漏分，通過無人機(jī)數(shù)據(jù)對比

分析，漏分面積為18709 hm2，漏分面積占實(shí)際面積 5 23774 hm2 的357%。357%與合并比估計(jì)的誤差 -333% 很接近，說明中分辨率數(shù)據(jù)的漏分精度直接影響合并比回歸估計(jì)的精度，那么在實(shí)際中中分辨率數(shù)據(jù)解譯水稻時(shí)要特別減少漏分。

3結(jié)論與討論

本研究使用中分辨率分類數(shù)據(jù)構(gòu)建抽樣框和分層，目的是在實(shí)際應(yīng)用中可指導(dǎo)無人機(jī)樣帶的布設(shè)。本研究在以下幾個(gè)方面值得探討：抽樣框覆蓋誤差，本研究方法采用中分辨率數(shù)據(jù)（SPOT5）解譯的水稻空間分布設(shè)計(jì)抽樣框，為了達(dá)到95%的農(nóng)業(yè)統(tǒng)計(jì)精度要求，根據(jù)上一段合并比回歸估計(jì)誤差和中分辨率數(shù)據(jù)（SPOT5）漏分誤差接近的分析，其漏分誤差保守估計(jì)要≤5%。

樣本量的確定：本研究采用基于無人機(jī)分類數(shù)據(jù)的總體方差S2y，計(jì)算得到最小樣本量是31個(gè)，轉(zhuǎn)為最小抽樣比是036%?；跓o人機(jī)分類數(shù)據(jù)的總體方差S2y在實(shí)際抽樣調(diào)查中常是難以獲取的，由于中分辨率結(jié)果和無人機(jī)結(jié)果具有很強(qiáng)的相關(guān)性，研究用中分辨率（SPOT5）的總體方差S2x替代無人機(jī)的總體方差時(shí)，計(jì)算得到的最小樣本量是34個(gè)，比原來的31個(gè)增加了3個(gè)，增漲了968%。

樣本量的分配：本研究采用等比例法分配樣本，將來試驗(yàn)可以考慮奈曼分配和調(diào)查費(fèi)用的最優(yōu)分配。

分別比估計(jì)和合并比估計(jì)比較：本研究中分別比估計(jì)精度不如合并比估計(jì)精度的結(jié)論，部分驗(yàn)證了倪加勛的結(jié)論[17]：分層分別回歸的效率高于分層合并回歸的估計(jì)效率，但在樣本量較小的情況下，分層分別回歸估計(jì)可能使偏差增大。

合并比估計(jì)回歸系數(shù)的確定：本研究中可以采用基于總體數(shù)據(jù)得到的準(zhǔn)確回歸系數(shù)，但考慮到本研究方法實(shí)際用于推廣，采用了基于樣本數(shù)據(jù)的回歸系數(shù)估計(jì)值。

本研究設(shè)計(jì)了基于中分辨率數(shù)據(jù)和無人機(jī)數(shù)據(jù)水稻種植面積分層回歸估計(jì)方法，該方法在抽樣框設(shè)計(jì)、樣本量確定、樣本量分配、分別比估計(jì)和合并比估計(jì)的選用、合并比回歸估計(jì)系數(shù)確定方面都給出了參考，可以為農(nóng)作物面積遙感抽樣調(diào)查方案設(shè)計(jì)提供參考。

參考文獻(xiàn)：

Chauhan H J，Arora M K，Agarwal A Estimating land cover class area from remote sensing classification[J] Journal of Applied Remote Sensing，2008，2（1）：183-198

張磊，吳炳方 關(guān)于土地覆被遙感監(jiān)測的幾點(diǎn)思考[J] 國土資源遙感，2011，3（1）：15-20

[3]劉慧平，朱啟疆 應(yīng)用高分辨率遙感數(shù)據(jù)進(jìn)行土地利用與覆蓋變化監(jiān)測的方法及其研究進(jìn)展[J] 資源科學(xué)，1999，21（3）：25-29

[4]Gallego F J Remote sensing and land cover area estimation[J] International Journal of Remote Sensing，2004，25（15）：3019-3047

[5]吳炳方，李強(qiáng)子 基于兩個(gè)獨(dú)立抽樣框架的農(nóng)作物種植面積遙感估算方法[J] 遙感學(xué)報(bào)，2004，8（6）：551-569

[6]吳炳方，蒙繼華，李強(qiáng)子 國外農(nóng)情遙感監(jiān)測系統(tǒng)現(xiàn)狀與啟示[J] 地球科學(xué)進(jìn)展，2010，25（10）：1003-1012

[7]范承嘯，韓俊，熊志軍，等 無人機(jī)遙感技術(shù)現(xiàn)狀與應(yīng)用[J] 測繪科學(xué)，2009，34（5）：214-215

[8]Laliberte A S，Herrick J E，Rango A，et al Acquisition，orthorectification，and object-based classification of unmanned aerial vehicle （UAV） imagery for rangeland monitoring[J] Photogrammetric Engineering & Remote Sensing，2010，76（6）：661-672

[9]Breckenridge R P Improving rangeland monitoring and assessment： integrating remote sensing，GIS，and unmanned aerial vehicle systems[D] Idaho Falls： University of Idaho，2007

[10]張煥雪，李強(qiáng)子，文寧，等 農(nóng)作物種植面積遙感抽樣調(diào)查的誤差影響因素分析[J] 農(nóng)業(yè)工程學(xué)報(bào)，2014（13）：176-184

[11]劉國棟，鄔明權(quán)，牛錚，等 基于GF-1衛(wèi)星數(shù)據(jù)的農(nóng)作物種植面積遙感抽樣調(diào)查方法[J] 農(nóng)業(yè)工程學(xué)報(bào)，2015，31（5）：160-166[HJ17mm]

[12]Gallego F J，Kussul N，Skakun S，et al Efficiency assessment of using satellite data for crop area estimation in Ukraine[J] International Journal of Applied Earth Observation and Geoinformation，2014，29：22-30

[13]王迪，周清波，陳仲新，等 玉米種植面積空間抽樣調(diào)查方案優(yōu)化設(shè)計(jì)[J] 農(nóng)業(yè)工程學(xué)報(bào)，2014，30（8）：117-125

[14]Claire B，Yang Z W，Rick M，et al Monitoring US agriculture： the US department of agriculture，national agricultural statistics service，cropland data layer program[J] Geocarto International，2011，26（5）：341-358[HJ]

[15]孫佩軍，張錦水，潘耀忠，等 基于無人機(jī)樣方事后分層的作物面積估算[J] 中國農(nóng)業(yè)資源與區(qū)劃，2016，37（2）：1-10

[16]杜子芳 抽樣技術(shù)及其應(yīng)用[M] 北京：清華大學(xué)出版社，2005：187-191

[17]倪加勛 抽樣調(diào)查[M] 桂林：廣西師范大學(xué)出版社，2003：61-84

[18]劉建紅，朱文泉 耕地變化空間抽樣調(diào)查方案的精度與效率分析[J] 農(nóng)業(yè)工程學(xué)報(bào)，2010，26（10）：331-336