運用Sentinel-1A時序數據和動態時間規整算法提取稻蝦田空間分布信息

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

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

關鍵詞：Sentinel-1A;動態時間規整算法;稻蝦田空間分布信息;TWDTW算法;VH極化;VV極化

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

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

收稿日期：2019-10-25

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

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

通信作者：熊勤學，碩士，副教授，主要研究方向為農業遙感。E-mail：nxqx@tom.com。

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

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

1?研究區概況和數據預處理

1.1?研究區概況

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

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

1.2?數據的預處理

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

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

1.3?分類結果驗證方法

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

Kappa=Po-Pc1-Pc。

其中：Po=s/n。

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

式中：n為柵格總像元;驗證區柵格數據中為1的象元數為a1，為0的象元數為a0;結果模擬數據中為1的象元數為b1，為0的象元數為b0;2個柵格對應象元值相等的象元數為s。Kappa計算結果中通常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?稻蝦田提取原理

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

2.2?稻蝦田提取方法

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

2.2.1?DTW計算方法

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

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

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

Ri=（ri-u）/σ。

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

（2）用歐氏距離方法計算出數據R中的每個點與數據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;

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

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

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

3?結果與分析

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

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

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

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

3.3?監利縣稻蝦田空間分布特征

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

4?結論與討論

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

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

隨著目前光學衛星數據空間分辨率的提高，面向對象的分類法開始運用[16]。面向對象的分類法既能提取對象的光譜信息，也能對對象的紋理、幾何特征、位置信息進行量化分類。因此，將光學衛星數據和SAR衛星數據進行融合，并結合TWDTW法進行面向對象的分類，是將來稻蝦田分類信息研究的發展方向。

參考文獻：

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

[2]農業農村部漁業漁政管理局全國水產技術推廣總站中國水產學會.中國小龍蝦產業發展報告2019[J]. 中國水產，2019，7（9）：12-19.

[3]曹湊貴，江?洋，汪金平，等. 稻蝦共作模式的“雙刃性”及可持續發展策略[J]. 中國生態農業學報，2017，25（9）：1245-1253.

[4]魏妍冰，陸?苗，吳文斌. 基于水體季相差異的稻蝦共作提取方法研究[J]. 中國農業資源與區劃，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.