盧澤如 熊勤學(xué) 周雨順 鄧琪云 紀紹威 丁璐

摘要：針對稻蝦田水稻移植期的不確定性給稻蝦田空間分布信息提取帶來的困難，運用SAR VH極化年時序數(shù)據(jù)、VV極化年時序數(shù)據(jù)以及調(diào)查區(qū)稻蝦田各自極化的年時序數(shù)據(jù)通過標準化處理后，進行TWDTW方法（time-weighted dynamic time warping）計算，得到二波段（VH極化、VV極化）與調(diào)查區(qū)稻蝦田年時序數(shù)據(jù)相似度空間分布柵格數(shù)據(jù)，再利用最大似然法進行監(jiān)督分類，得到2016—2018年3年的監(jiān)利縣稻蝦田空間分布柵格數(shù)據(jù)，分類結(jié)果與4塊區(qū)域（龔場鎮(zhèn)、福田寺鎮(zhèn)、紅城鄉(xiāng)、白螺鎮(zhèn)，面積分別是373.33、2 200.00、126.67、240.00 hm2）稻蝦田實際分布圖比較，其Kappa系數(shù)分別是0.79、0.80、0.75、0.90，介于0.75～0.90之間，說明運用Sentinel-1A時序數(shù)據(jù)和動態(tài)時間規(guī)整算法的提取方法是準確的。由于SAR傳感器具有全天侯以及不受云層影響等特點，能定期獲取研究區(qū)上空數(shù)據(jù)，相比光學(xué)傳感器獲得無云或者少云數(shù)據(jù)的不確定性，此方法能在實際操作中運用。

關(guān)鍵詞：Sentinel-1A;動態(tài)時間規(guī)整算法;稻蝦田空間分布信息;TWDTW算法;VH極化;VV極化

中圖分類號：S127?文獻標志碼： A

文章編號：1002-1302（2020）18-0230-07

收稿日期：2019-10-25

基金項目：國家自然科學(xué)基金（編號：31871516）;公益性行業(yè)（農(nóng)業(yè)）科研專項（編號：201203032）;濕地生態(tài)與農(nóng)業(yè)利用教育部工程研究中心開放基金（編號：KF201701、KFT201906）;長江大學(xué)第十一批大學(xué)生創(chuàng)新創(chuàng)業(yè)訓(xùn)練計劃（編號：2018288）。

作者簡介：盧澤如（1999—），女，湖北黃岡人，主要研究方向為農(nóng)業(yè)信息。E-mail：1284563860@qq.com。

通信作者：熊勤學(xué)，碩士，副教授，主要研究方向為農(nóng)業(yè)遙感。E-mail：nxqx@tom.com。

近幾年來，被農(nóng)業(yè)農(nóng)村部譽為“現(xiàn)代農(nóng)業(yè)發(fā)展的成功典范”的稻蝦共作生態(tài)種養(yǎng)高效模式在全國飛速發(fā)展[1]。根據(jù)全國水產(chǎn)技術(shù)推廣總站發(fā)布的《2017中國小龍蝦產(chǎn)業(yè)發(fā)展報告》，全國小龍蝦（克氏原螯蝦，Procambarus clarkii）養(yǎng)殖面積超過80萬hm2[2]。如何研究稻蝦共作模式的生態(tài)效應(yīng)，為稻蝦共作模式的可持續(xù)發(fā)展提供依據(jù)，是當(dāng)前亟待解決的問題[3]。基于遙感影像數(shù)據(jù)提取稻蝦田空間分布信息，被視為最便捷、最高效的研究手段，而目前對稻蝦田空間分布信息的遙感提取研究還不多。魏妍冰等根據(jù)稻蝦田獨特的水域變化特征，利用自動水域提取指數(shù)（automated water extraction index）分別獲取不同季節(jié)的水體，最后根據(jù)水體季相差異特征，用2017年8月和12月2景Landsat 8 OLI影像數(shù)據(jù)提取了2017年潛江市稻蝦田的空間分布;這種冬季與夏季相比會增加的水域，視為稻蝦共作田的方法有明顯的異物同譜的錯誤，如低洼地帶的中稻田和水生植被區(qū)[4]。造成提取稻蝦田空間分布信息困難的原因主要有3個：一是同一地區(qū)稻蝦田的水稻移植期差異很大。小龍蝦市場價格因素和優(yōu)良的水熱氣候條件使得水稻可生長的日期遠大于實際生長期，水稻移植期為5月中旬至7月中旬，過長的水稻移植期導(dǎo)致不同稻蝦田光譜特征差異很大，增加了遙感影像數(shù)據(jù)的提取難度;二是長重訪周期和雨熱同季的氣候特征限制了光學(xué)衛(wèi)星數(shù)據(jù)真實表達物候光譜特征的可能性[5];三是復(fù)雜的種植制度和粗光學(xué)遙感影像空間分辨率造成的混合像素嚴重制約信息提取的精度。

合成孔徑雷達因不受光照和氣候等條件限制，具備實現(xiàn)全天時、全天候?qū)Φ赜^測的特點，特別是歐空局Sentinel-1A衛(wèi)星C波段合成孔徑雷達（SAR），實現(xiàn)了多極化（HV、HH極化）、高重訪周期（12 d）、高空間分辨率[干涉寬幅（IM）模式：250 km，5 m×20 m分辨率]對地觀測，海量數(shù)據(jù)加上免費政策，為稻蝦田空間分布信息的提取提供了足夠的基礎(chǔ)數(shù)據(jù)。本研究意在以全國稻蝦種養(yǎng)面積最大的縣——湖北省監(jiān)利縣為研究對象，在充分分析稻蝦田Sentinel-1A衛(wèi)星C波段SAR數(shù)據(jù)時序特征的基礎(chǔ)上，改進能夠?qū)r序數(shù)據(jù)進行動態(tài)時間歸整（dynamic time warping，DTW）法，來消除長時間水稻移植期對不同稻蝦田光譜特征差異的影響，高精度地提取監(jiān)利縣稻蝦田空間分布信息。

1?研究區(qū)概況和數(shù)據(jù)預(yù)處理

1.1?研究區(qū)概況

監(jiān)利縣（圖1）位于湖北省中南部，長江北岸，面積3 508 km2，四季分明，熱量豐富，光照適宜，雨水充沛，雨熱同季，無霜期長，是全國稻蝦種養(yǎng)面積最大的縣（2017年達到3.33萬hm2）[2]，也是水稻種植面積最大的縣（市）。

2017年8月，為分析稻蝦田SAR數(shù)據(jù)的時序特征和驗證結(jié)果的精度，在監(jiān)利縣的龔場鎮(zhèn)、福田寺鎮(zhèn)、紅城鄉(xiāng)、白螺鎮(zhèn)分別選取4塊驗證區(qū)，并用GPS儀確定了每塊稻蝦田的邊界（圖1），4塊區(qū)域稻蝦田種植的面積分別是373.33、2 200.00、126.67、240.00 hm2。

1.2?數(shù)據(jù)的預(yù)處理

從歐洲航天局的網(wǎng)站（https：//scihub.copernicus.eu/dhus/#/home）上下載2016年1月至2018年12月涵蓋監(jiān)利縣區(qū)域的Sentinel-1A衛(wèi)星C波段SAR GRDH的格式數(shù)據(jù)，2景數(shù)據(jù)能涵蓋監(jiān)利縣所有區(qū)域，共下載79 d共158景（2016年40景、2017年56景、2018年62景）數(shù)據(jù)，運用ESA Sentinel 1 Toolbox（Ver 1.1.1）軟件做數(shù)據(jù)的預(yù)處理，數(shù)據(jù)預(yù)處理步驟如下：

（1）SAR數(shù)據(jù)合并：監(jiān)利縣跨越2景數(shù)據(jù)，因此在ESA Sentinel 1 Toolbox軟件中打開同一天的2景GRDH數(shù)據(jù)，運用SAR MOSAIC功能，將2景數(shù)據(jù)合并;（2）SAR數(shù)據(jù)裁剪：加載SHAPEFILE格式的監(jiān)利縣行政區(qū)劃電子數(shù)據(jù)，在Navigation窗口中縮放到合適位置后，利用spatial subset from view功能將涵蓋監(jiān)利縣的數(shù)據(jù)裁剪下來;（3）后向散射系數(shù)計算：運用Radiometric下的Calibrate功能，選擇output sigma band選項，計算VV、VH 2種極化、裁剪區(qū)域的后向散射系數(shù);（4）圖像降噪和濾波處理：為降低圖像噪聲水平，運用Multi-Looking功能進行多視處理，運用Speckle Filter功能，選取Refine Lee濾波方式進行濾波處理;（5）地形校正：運用Range Doppler Terrain Correction功能，選取SRTM 3Sec數(shù)字高程，坐標系統(tǒng)選取UTM WGS84 49N進行高程校正;（6）數(shù)據(jù)導(dǎo)出：將地形校正后的數(shù)據(jù)以ENVI格式導(dǎo)出，生成VV、VH 2波段數(shù)據(jù);（7）多時序VV、VH極化數(shù)據(jù)生成：將2016—2018年數(shù)據(jù)按上述步驟進行處理，在ENVI軟件中按年打開所有處理數(shù)據(jù)，運用Layer Stacking功能，將1年內(nèi)所有VV極化或VH極化數(shù)據(jù)按日期順序排列，生成一個多波段時序數(shù)據(jù)，波長用DOY表示（一年中的第n天），這樣一年中有2個多波段時序遙感數(shù)據(jù)，一個是VV極化的，另一個是VH極化的。

1.3?分類結(jié)果驗證方法

將結(jié)果數(shù)據(jù)和驗證區(qū)調(diào)查空間數(shù)據(jù)運用GIS中的聚類分析生成只含0和1數(shù)字的柵格數(shù)據(jù)，0代表非稻蝦區(qū)，1代表稻蝦區(qū)，然后進行Kappa系數(shù)計算，公式[6]如下：

Kappa=Po-Pc1-Pc。

其中：Po=s/n。

Pc=（a1×b1+a0×b0）。

式中：n為柵格總像元;驗證區(qū)柵格數(shù)據(jù)中為1的象元數(shù)為a1，為0的象元數(shù)為a0;結(jié)果模擬數(shù)據(jù)中為1的象元數(shù)為b1，為0的象元數(shù)為b0;2個柵格對應(yīng)象元值相等的象元數(shù)為s。Kappa計算結(jié)果中通常Kappa值是在0～1之間，0.0～0.20表示極低一致性，0.21～0.40表示一般一致性，0.41～0.60表示中等一致性，0.61～0.80表示高度一致性，0.81～1.00表示幾乎完全一致。

2?稻蝦田提取原理與方法

2.1?稻蝦田提取原理

影響地物后向散射系數(shù)的要素很多，有地物的粗糙度、土壤含水量、作物結(jié)構(gòu)（生物量、葉面積系數(shù)、葉面、莖桿粗細等）、極化方式、頻率和入射角等[7-8]。相同地物的后向散射系數(shù)變化規(guī)律基本一致，這是利用微波后向散射系數(shù)的時序特征進行地物分類的基本原理。水稻田后向散射系數(shù)的時序特征非常明顯，即營養(yǎng)生長期（5—6月），因水稻淹水，后向散射系數(shù)偏低，與水體的時序特征相似，而進入生殖生長期（7—9月），其后向散射系數(shù)明顯增加，到成熟早期達到最大[9];無論是稻蝦共作還是連作方式種養(yǎng)的稻蝦田，水稻收獲后到下輪水稻移栽期（11月至次年4月），為了方便養(yǎng)殖小龍蝦，稻田里一直會保持適當(dāng)水位的水，因此這段時間的后向散射系數(shù)和水體一樣，是偏低的。圖2為2017年稻蝦田、傳統(tǒng)水稻田、水體、旱地、城市5種地物VH極化后向散射系數(shù)的年變化曲線。水體后向散射系數(shù)一年四季都很低（小于0.01），而城市后向散射系數(shù)一直偏高（大于0.12），旱田后向散射系數(shù)介于0.03～0.07之間，而傳統(tǒng)稻田4—6月后向散射系數(shù)小于0.02，其他時間介于0.03～0.07之間。稻蝦田后向散射系數(shù)的時序變化特征較其他地物有明顯不同，其后向散射系數(shù)只會在6—9月偏高，其他時間都偏低。

2.2?稻蝦田提取方法

盡管稻蝦田后向散射系數(shù)的時序變化特征較其他地物有明顯不同，但因為監(jiān)利縣光熱資源豐富，大于20 ℃的氣溫持續(xù)133 d，屬雙季稻種植區(qū)，水稻可生長期大于實際生長期。而不同種植戶因受龍蝦養(yǎng)殖情況和市場因素影響，插秧期差異很大，從5月底至7月初都有可能插秧。2017年在監(jiān)利縣調(diào)查的4個地點，插秧期為5月23日至7月2日，每個點的HV極化后向散射系數(shù)的時序變化曲線有明顯的平移變化（圖3），而且最大值也有差異，因此常規(guī)的遙感分類方法不適用于稻蝦田的空間分布信息提取。而DTW算法可以計算2個時間序列的相似度，尤其適用于不同長度、不同節(jié)奏的時間序列，能消除插秧期差異帶來的后向散射系數(shù)時序變化曲線的平移影響。因此，本研究擬采用DTW方法提取稻蝦田的空間分布信息。

2.2.1?DTW計算方法

DTW方法最初用來辨別2組語音的相似性，即計算2組時間序列的相似度，尤其適用于不同長度、不同節(jié)奏時間序列相似度的計算，近幾年常用于遙感數(shù)據(jù)時序分析中。Petitjean等從理念的角度介紹了DTW方法，并說明了如何應(yīng)用在遙感數(shù)據(jù)時序分析中[10]。Baumann等運用DTW方法消除了年際間物候差異對NDVI指數(shù)時序的影響，并成功對森林進行分類[11]。Gristina等利用DTW方法識別出Landsat NDVI時間序列中的異常值，并用時空插值方法代替，從而改進森林分類精度[12]。Maus等將時間限制條件引入DTW方法中，提出TWDTW方法（time-weighted DTW）來提高分類精度[13]。由此可知，DTW主要用來計算遙感數(shù)據(jù)歸一化指數(shù)時序數(shù)據(jù)的相似度。為將DTW方法運用于SAR后向散射系數(shù)相似性，提出先將數(shù)據(jù)進行標準化處理，然后再進行TWDTW的方法，具體計算方法如下。

如計算2個點后向散射系數(shù)時序數(shù)據(jù)r=（r1，r2，…，rn）與 t=（t，t2，…，tm）之間的相似度。

（1）對后向散射系數(shù)時序數(shù)據(jù)進行標準化處理：

Ri=（ri-u）/σ。

u、σ分別為數(shù)組的r平均值和方差。

（2）用歐氏距離方法計算出數(shù)據(jù)R中的每個點與數(shù)據(jù)T中各點的距離，生成一個n行m列的距離矩陣D[n，m]。

（3）DTW相似度計算

DTW（i，j）=DTW（i-1，j-1）+min{D（i-1，j），D（i-1，j-1），D（i，j-1）}。

其中DTW（0，0）=0;

當(dāng)i-j>K時，DTW（i，j）=0，K為時間限制，根據(jù)稻蝦田中水稻移栽期不確定的特點，K值確定為50 d。

當(dāng)i或者j等于n或者m時，DTW（i，j）為最終2組時間序列的相似度。

2.2.2?稻蝦田提取計算流程?稻蝦田提取計算具體流程見圖4。

3?結(jié)果與分析

3.1?時序數(shù)據(jù)標準化處理對稻蝦田空間分布信息提取的影響

圖5為調(diào)查點2經(jīng)過標準化處理與未標準化處理的TWDTW計算結(jié)果空間分布圖，其值越低，表明2組數(shù)據(jù)的相似度越高，圖上顯示為黑色。從圖5可以看出，未標準化處理的TWDTW計算結(jié)果（圖5-B）與稻蝦田空間分布差異很大，主要原因是影響SAR后向散射系數(shù)的因子很多，除了地物差異外，還有土壤水分、土壤類型、植被信息等，導(dǎo)致盡管其時序特征曲線非常相似，但不同區(qū)域的時序曲線中最高值、最低值差異很大，直接運用TWDTW計算，其他特征引起的差異高于時序特征引起的差異;而時序特征曲線經(jīng)過標準化處理后，將時序特征曲線置于同一量級，能準確反映物候引起的時序差異（圖5-A）。因此，對SAR時序數(shù)據(jù)作TWDTW處理之前，一定要進行標準化處理。

3.2?稻蝦田空間分布信息提取結(jié)果驗證

圖6為4塊區(qū)域（龔場鎮(zhèn)、福田寺鎮(zhèn)、紅城鄉(xiāng)、白螺鎮(zhèn)，面積分別是373.33、2 200.00、126.67、240.00 hm2）稻蝦田實際分布圖和提取結(jié)果圖。通過計算得出它們的Kappa系數(shù)分別是0.79、0.80、0.75、090，介于0.75～0.90之間，說明運用TWDTW計算后，再用最大似然法提取的結(jié)果還是能真實反映稻蝦田空間分布的。從對比圖中可以看出，其主要誤差是邊界模糊不清和稻蝦田有中間鏤空現(xiàn)象，這主要是由于SAR圖像中存在著顯著的相干斑噪聲[14]，盡管圖像進行過降噪和濾波處理，但它對提取結(jié)果的影響還是很大。

3.3?監(jiān)利縣稻蝦田空間分布特征

監(jiān)利縣稻蝦田種養(yǎng)面積飛速發(fā)展，2016年面積為1.88萬hm2，2017年為3.22萬hm2，至2018年為4.71萬hm2，與監(jiān)利縣農(nóng)業(yè)統(tǒng)計年鑒的統(tǒng)計數(shù)據(jù)（2016年種養(yǎng)面積為2.00萬hm2、2017年為3.33萬hm2、2018年為5.00萬hm2）基本一致。圖7為2016—2018年監(jiān)利縣稻蝦田空間分布，可以看出，監(jiān)利縣稻蝦田主要分布在洪湖西邊的低洼地帶（監(jiān)利縣的東邊），沿長江民垸區(qū)，也就是由以前傳統(tǒng)的水稻田和魚池改建而來。2018年稻蝦田種養(yǎng)面積超過0.27萬hm2的鄉(xiāng)鎮(zhèn)有黃歇口鎮(zhèn)、白螺鎮(zhèn)、棋盤鄉(xiāng)、上車灣鎮(zhèn)、尺八鎮(zhèn)。2016—2018年各鄉(xiāng)鎮(zhèn)稻蝦田種養(yǎng)面積統(tǒng)計見圖8。

4?結(jié)論與討論

本研究運用SAR VH極化年時序數(shù)據(jù)和VV極化年時序數(shù)據(jù)，與驗證區(qū)稻蝦田的相同極化年時序數(shù)據(jù)進行TWDTW計算后，得到二波段DTW數(shù)據(jù)，再用最大似然法成功提取2016—2018年監(jiān)利縣稻蝦田空間分布信息，其Kappa系數(shù)介于 0.75～0.90之間，與實際狀況基本一致。由于SAR傳感器具有全天侯監(jiān)測以及不受云層影響等特點，能定期獲取研究區(qū)上空數(shù)據(jù)，相比光學(xué)傳感器獲得無云或者少云數(shù)據(jù)的不確定性，此方法可在實際操作中運用。由于數(shù)據(jù)來源于VV極化和VH極化數(shù)據(jù)，相比用單極化數(shù)據(jù)提取，其結(jié)果精度有一定的提高。TWDTW法能夠?qū)r間軸進行動態(tài)拉伸或壓縮，對作物種植時間存在提前或延后的耕地時間序列也有較好的適應(yīng)性，特別適合稻蝦田這種水稻移栽期不確定的農(nóng)田分類。本研究提取稻蝦田空間分布信息的方法，克服了前言中提到的提取稻蝦田空間分布信息的難點，不失為一種行之有效的方法。

不可否認，盡管SAR數(shù)據(jù)的空間分辨率和時間分辨率較以前有很大的提升，但運用SAR數(shù)據(jù)進行分類，與同級別光學(xué)數(shù)據(jù)的分類結(jié)果相比，精度偏低[15]。主要誤差來源是SAR圖像相干斑噪聲，這會導(dǎo)致分類結(jié)果的邊界模糊不清;而光學(xué)衛(wèi)星數(shù)據(jù)盡管獲取的概率低，但如果能獲得無云數(shù)據(jù)，可以得到清晰的邊界。因此，如果合理地運用兩者之間的優(yōu)點，則既能有清晰的邊界，又能利用SAR定時獲取數(shù)據(jù)的特點而得到高精度的分類結(jié)果，是提高稻蝦田空間分布信息精度主要的研究方向。

隨著目前光學(xué)衛(wèi)星數(shù)據(jù)空間分辨率的提高，面向?qū)ο蟮姆诸惙ㄩ_始運用[16]。面向?qū)ο蟮姆诸惙饶芴崛ο蟮墓庾V信息，也能對對象的紋理、幾何特征、位置信息進行量化分類。因此，將光學(xué)衛(wèi)星數(shù)據(jù)和SAR衛(wèi)星數(shù)據(jù)進行融合，并結(jié)合TWDTW法進行面向?qū)ο蟮姆诸?，是將來稻蝦田分類信息研究的發(fā)展方向。

參考文獻：

[1]秦尊文. 以“蝦稻共作”模式為抓手推進體制機制創(chuàng)新——潛江市全國中小城市綜合改革的觀察與思考[J]. 中國發(fā)展，2016，160（6）：51-56.

[2]農(nóng)業(yè)農(nóng)村部漁業(yè)漁政管理局全國水產(chǎn)技術(shù)推廣總站中國水產(chǎn)學(xué)會.中國小龍蝦產(chǎn)業(yè)發(fā)展報告2019[J]. 中國水產(chǎn)，2019，7（9）：12-19.

[3]曹湊貴，江?洋，汪金平，等. 稻蝦共作模式的“雙刃性”及可持續(xù)發(fā)展策略[J]. 中國生態(tài)農(nóng)業(yè)學(xué)報，2017，25（9）：1245-1253.

[4]魏妍冰，陸?苗，吳文斌. 基于水體季相差異的稻蝦共作提取方法研究[J]. 中國農(nóng)業(yè)資源與區(qū)劃，2019，40（3）：14-20，34.

[5]Guan X D，Chong H，Liu G H，et al. Mapping rice cropping systems in Vietnam using an NDVI-based timeseries similarity measurement based on DTW distance[J]. Remote Sensing，2016，8（1）：19.

[6]Bradley A P. The use of the area under the ROC curve in the evaluation of machine learning algorithms[J]. Pattern Recognition，1997，30（7）：1145-1159.

[7]Graham A J，Harris R. Extracting biophysical parameters from remotely sensed radar data：a review of the water cloud model[J]. Progrom Physics Geography，2003，27（2）：217-229.

[8]Inoue Y，Kurosu T，Maeno H，et al. Season-long daily measurements of multifrequency （Ka，Ku，X，C，and L） and full-polarization backscatter signatures over paddy rice field and their relationship with biological variables[J]. Remote Sensing of Environment，2002，81（2/3）：194-204.

[9]Toan T L，Laur H，Mougin E，et al. Multitemporal and dual-polarization observations of agricultural vegetation covers by X-band SAR images[J]. IEEE Transactions on Geoscience and Remote Sensing，1989，27（6）：709-718.

[10]Petitjean F，Inglada J，Gancarski P. Satellite image time series analysis under time warping[J]. IEEE Transactions on Geoscience and Remote Sensing，2012，50（8）：3081-3095.

[11]Baumann M，Mutlu O，Richardson A D，et al. Phenology from landsat when data is scarce：using MODIS and dynamic time-warping to combine multi-year landsat imagery to derive annual phenology curves[J]. International Journal of Applied Earth Observation and Geoinformation，2017，54（2）：72-83.

[12]Cristina G，Joanne C W，Michael A W，et al. Integrated object-based spatiotemporal characterization of forest change from an annual time series of landsat image composites[J]. Canadian Journal of Remote Sensing，2015，41（4）：271-292.

[13]Maus V，Cmara G，Cartaxo R，et al. A Time-weighted dynamic time warping method for land-use and land-cover Mapping[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing， 2016，9（8）：3729-3739.

[14]Moreira A. Improved multi-look techniques applied to SAR and scan SAR images[J]. IEEE Trans on Geoscience and Remote Sensing，1991，23（4）：529-534.

[15]McNairn H，Shang J，Jiao X，et al. The contribution of ALOS PALSAR multi-polarization and polarimetric data to crop classification[J]. IEEE Geoscience and Remote Sensing，2009，47（12）：3981-3992.

[16]Ovidiu C，Mariana B，Gregory P A，et al. Object-based time-constrained dynamic time warping classification of crops using sentinel-2[J]. Remote Sensing，2019，11（10）：1257.