孔 麗，劉美玲，劉湘南，鄒信裕

（中國地質大學（北京）信息工程學院，北京 100083）

0 引 言

隨著經濟的不斷發(fā)展，工業(yè)化不斷推進，土壤重金屬污染成為當今世界面臨的重大生態(tài)環(huán)境問題之一[1]。水稻是中國的主要農作物，產量位居世界第二位。水稻重金屬污染不僅會影響水稻的糧食產量，同時也會威脅人類的身體健康。傳統(tǒng)的實地調查雖然能夠準確檢測水稻中的重金屬濃度[2]，但費時費力，且難以滿足大區(qū)域尺度監(jiān)測的要求。由于遙感技術能夠實時、大范圍、連續(xù)動態(tài)獲取目標地物對環(huán)境脅迫的光譜響應信息，因而利用遙感技術對水稻重金屬污染現(xiàn)狀進行監(jiān)測和評估逐步受到國內外學者的廣泛關注。

目前，許多學者對水稻重金屬脅迫遙感監(jiān)測開展了一系列的研究。趙利婷等[3]利用遙感同化WOFOST（World Food Studies）模型模擬水稻根重進行水稻重金屬污染脅迫分析。Liu等[4]利用多時相的Sentinel-2數據融合光譜信息和時空信息監(jiān)測水稻的重金屬脅迫。近年來，為提高農作物重金屬遙感識別精度，一些學者利用了多尺度小波變換的光譜分析方法對農作物污染脅迫狀況進行識別和評價[5-6]；另一些學者將農作物生理功能特征變化（如蒸騰作用、光合作用、呼吸作用等）在光譜上的響應作為農作物重金屬脅迫評估指標來揭示農作物重金屬污染狀況，如Jin等[7]構建了基于農作物形態(tài)和生理功能變化參數的重金屬脅迫理論模型。Tian等[8]將脅迫效應的譜特征與時間特征結合去識別重金屬脅迫。Tang等[9]提出了基于時空特征指標的水稻重金屬脅迫識別模型，構建了融合水稻脅迫的年內特征、年際特征以及空間特征的時空指標，實現(xiàn)了水稻重金屬脅迫的監(jiān)測。上述農作物重金屬污染監(jiān)測研究主要是在較小的空間尺度且已知重金屬污染條件下，假設光譜變異僅由重金屬脅迫引起（忽略其他脅迫影響），如何剝離其他脅迫對重金屬脅迫的影響是開展重金屬脅迫遙感準確識別亟待解決的問題。

在復雜的農田生態(tài)系統(tǒng)中，農作物會受到復雜脅迫的影響，包括短期脅迫（如病蟲害、水分脅迫等）和長期脅迫（如重金屬脅迫等）[10]。據研究表明，重金屬污染具有隱蔽性、持久性和不可逆性等特點[11-12]，因此，重金屬脅迫存在于農作物整個生長周期，對農作物造成的影響是持續(xù)穩(wěn)定的，且多年都呈現(xiàn)一個穩(wěn)定相似的狀態(tài)。而短期脅迫的特點是年際間差異較為明顯，往往存在于農作物的某一個或某幾個生長周期中。由此，農作物不同脅迫效應呈現(xiàn)的時間特征差異為區(qū)分農作物重金屬脅迫和其它脅迫提供了可行性，探索一種適合不同脅迫特征分量的分離與提取的時序信號分解算法，是水稻重金屬脅迫準確識別的關鍵。水稻脅迫特征參數通常包含水稻固有生長趨勢，年內長期脅迫和隨機噪聲或年內短期脅迫。時序信號分解可以有效地獲取感興趣的信號，集合經驗模態(tài)分解算法（Ensemble Empirical Mode Decomposition，EEMD）是一種自適應時序信號分解方法，可以很好地揭示信號的非平穩(wěn)性和非線性信息。該方法在信號分析中得到了廣泛地應用[13-15]。

為揭示水稻不同脅迫呈現(xiàn)的時間效應差異，有必要開展水稻脅迫特征參數的時間連續(xù)模擬與計算。WOFOST模型根據氣象和土壤條件以一天時間為間隔[16]模擬作物的葉面積指數（Leaf Area Index，LAI）生成密集的時間序列，準確實現(xiàn)環(huán)境脅迫下作物生長過程的動態(tài)模擬，從而獲取時間連續(xù)的作物脅迫特征參數，為不同脅迫特征分量的分離與提取提供基礎數據。陳艷玲等[17]采集冬小麥多個關鍵生育時期的生理生化、農田環(huán)境、氣象等數據，構建基于WOFOST模型和遙感LAI數據同化的區(qū)域尺度冬小麥單產預測模型。黃健熙等[18]重點優(yōu)化WOFOST模型中與品種相關的積溫參數，該模型在全國尺度取得了較高的模擬精度。以上研究表明了遙感同化WOFOST模型在監(jiān)測作物的長勢等方面具有一定的優(yōu)勢。針對水稻重金屬脅迫具有“年際穩(wěn)定性”特征，因而水稻脅迫特征參數在不同年份存在相似性。動態(tài)時間規(guī)整（Dynamic Time Warping，DTW）方法衡量時間序列曲線的相似度，該算法已經相當成熟，應用廣泛，例如語音識別、信息檢索、土地覆蓋分類等等[19-21]，這些研究成果為重金屬脅迫研究提供了重要的參考。相對于傳統(tǒng)的距離度量方法，DTW更適合度量不同長度以及不同節(jié)奏的時間序列的相似度，這一優(yōu)勢正好可以有效避免不同年份水稻物候期（如分蘗期）的差異性所引起的相似性計算誤差，由此，DTW算法適合度量水稻脅迫特征參數在不同年份的相似度?；谇懊娴姆治觯狙芯恳院鲜≈曛逓樵囼瀰^(qū)域，利用多年Sentinel-2數據同化作物生長模型得到不同年份的LAI時間序列，通過EEMD方法對其進行分解，剝離年內短期脅迫或噪聲，提取可能包含重金屬脅迫在內的年內長期脅迫；在此基礎上，利用DTW計算其年內長期脅迫信號與健康水稻脅迫信號之間的距離，通過相似性度量準確識別水稻重金屬脅迫。

1 材料與方法

1.1 試驗區(qū)概況

試驗區(qū)位于湖南省株洲市（113°05′E～113°15′ E，27°37′N～27°45′N）（圖1）。該地區(qū)屬于亞熱帶季風性濕潤氣候，四季分明，雨量充沛，光照充足，年平均氣溫為16～18 ℃，年降水量大約為1 445 mm，土壤為富含多種營養(yǎng)物質的紅壤，有機質含量（質量分數）2%～3%，適宜水稻生長，是國家重要的商品糧食生產基地。本試驗區(qū)覆蓋一幅Sentinel-2影像的大小為1 000×1 000像元，其范圍約為100 km2。據調查，湘江流域灌溉的土壤中重金屬鎘的含量較高，是主要污染物[22-23]。

1.2 數據獲取與處理

本研究所采用的數據包括遙感數據、作物生長模型所需的輸入數據以及實測的土壤重金屬含量數據。其中作物生長模型所需的數據為氣象數據和作物生長參數，氣象數據來自中國氣象站（http://www.cma.gov.cn），包括逐日最低溫、逐日最高溫以及逐日照時數。而作物參數主要有物候因子、作物初始參數、綠葉面積、CO2同化和干物質分配。土壤重金屬數據針對25個采樣點均選取距土壤表層10 cm深處的土壤，曬干后使用HNO3-HF-HCLO4消毒法在微波消解儀（MARS X，CEM，USA）進行消毒，并采用電感耦合等離子體質譜ICP-MS（7500a，Agilent Technologies，USA）測定Cd、Pb等元素含量。每個土壤樣品平均測定3次，結果以均值表示。遙感數據是2017—2019年的Sentinel-2衛(wèi)星數據（分別為2017年7月17日和24日、2017年8月6日、2017年9月15日；2018年6月17日、2018年7月17日、2018年8月8日、2018年9月10日；2019年6月14日、2019年8月16日、2019年9月15日、2019年9月22日）（https://scihub.copernicus.eu），首先對遙感影像數據進行預處理，包括輻射校正和大氣校正等，然后在遙感影像中選取訓練樣本，利用隨機森林算法提取研究區(qū)域的水稻分布范圍。

1.3 研究方法

本文采用作物生長模型、EEMD與DTW方法有效地甄別并定量分析水稻重金屬脅迫，具體的流程如圖2所示。主要分為4個步驟：1）對數據進行預處理，基于遙感影像數據獲取研究區(qū)域水稻的LAI；2）運用作物生長模型模擬水稻的LAI時間序列；3）通過EEMD方法將原始LAI時間序列分解成不同時間尺度的組分，分析合成含有重金屬脅迫的信號；4）計算分解后含有重金屬脅迫的信號與健康水稻之間的DTW距離并進行歸一化處理，監(jiān)測重金屬脅迫的程度。

1.3.1 基于遙感影像反演水稻LAI

LAI是植被主要的理化參數之一，本文利用LAI作為遙感與作物生長模型之間的同化量。目前，獲取LAI的方法主要有經驗模型、物理模型與混合模型，經驗模型主要是通過建立實地觀測數據與遙感影像數據之間的數量關系來獲取LAI觀測值。統(tǒng)計模型操作簡單，易于實現(xiàn)，精度較高，因此，本文利用實測LAI和Sentinel-2數據獲取的歸一化植被指數（NDVI）值建立適合本研究區(qū)的LAI經驗反演模型。對比不同的擬合方程選取擬合精度最佳的擬合模型[24]。具體的LAI計算式如下：

式中NIR和R分別代表近紅外波段和紅波波段反射率。

1.3.2 基于WOFOST模型獲取長時序LAI

為了將不連續(xù)的遙感信息轉化為時間連續(xù)的作物信息，本文利用遙感與WOFOST作物生長模型同化方法模擬作物的生長過程。該模型基于作物生長發(fā)育的基本過程，解釋了作物的生長，如光合作用以及呼吸作用。它能夠模擬溫度、水分對作物生長發(fā)育形成的脅迫時期和程度[25]。WOFOST模型主要有3種水平的作物生長：潛在生產水平、水分限制生產水平和營養(yǎng)限制生產水平。它的輸入參數有作物、土壤和氣候參數，本文研究綜合影響水稻長勢的各方面因素，包括穩(wěn)定脅迫（如重金屬脅迫等）、短時脅迫（病蟲害等）以及特定時間發(fā)生的災害（干旱等），依據脅迫對作物生長的影響機理，選擇合適的脅迫因子嵌入到日總CO2同化的生產過程。

式中CVFf為脅迫水平下的日總CO2同化物產量，kg/hm2；CVF代表潛在生產水平下的日總CO2同化物產量，kg/hm2；f為脅迫因子，取值范圍為[0,1]，f值越高，說明作物生長過程受到的脅迫越小，反之，脅迫越嚴重。

同化算法是模型同化的關鍵，同化算法的選擇直接影響模型模擬的效率和精度。本文采用原理簡單、易于實現(xiàn)的粒子群優(yōu)化（Particle Swarm Optimization，PSO）算法[26-27]，其主要的思想是通過不斷調整模型參數CVR（Efficiency of Conversion into Roots，干物質轉化為根質量的效率），使得WOFOST模型模擬值與遙感反演實測值之間的差異逐漸減小，收斂條件是全局最優(yōu)解不再變化或達到最大迭代次數，最終輸出CVR的最優(yōu)解即同化結果。本文使用的代價函數為

式中LAIm為實測值序列；LAIs為模擬值序列；N為實測時相數。

通過比較模型模擬的LAI與Sentinel-2數據反演的LAI，得出優(yōu)化后的模擬結果更接近LAI真實值，通過WOFOST模型模擬出連續(xù)的水稻脅迫參數LAI。

1.3.3 基于EEMD分解LAI時間序列

EEMD原理較為簡單，該方法是針對經驗模態(tài)分解（Empirical Mode Decomposition，EMD）方法的不足向原始信號中添加白噪聲，使其信號在不同尺度上具有連續(xù)性，從而避免出現(xiàn)模態(tài)混疊的現(xiàn)象，達到一個更好的分解結果。EEMD分解方法兩個非常重要的參數是白噪聲的幅值系數和總體平均次數，白噪聲的幅值系數需要人為經驗進行確定，添加的白噪聲幅值系數過大，會降低分解的精度，過小，則無法解決EMD中存在的模態(tài)混疊現(xiàn)象。Huang等[28]提出分解以高頻組分為主的原始信號時，白噪聲幅值系數應相對減小，相反，白噪聲幅值系數應相對增大。通常白噪聲幅值系數的取值范圍為0.01～0.40?？傮w平均次數會影響分解效果，為了同時兼顧效率與精度，總體平均次數一般選擇100～300。參考前人的研究結果，結合本論文試驗的需求，將總體平均次數設置為100次，白噪聲幅值系數設置為0.3[10]。

時間序列經過EEMD分解后可以得到不同分量[29]。這些分量由3個波動組分和1個殘差組成，其中3個組分是年內波動組分、年間波動組分和年際波動組分[30-31]。年內波動組分為時間尺度通常小于一年的高頻信號。高頻信號可以檢測信號中存在的異?，F(xiàn)象，其他環(huán)境脅迫和噪聲信息存在于高頻信號中。作物固有的生長趨勢特征與季節(jié)項相似，時間尺度通常為一年，有一個固定的周期[32]。本文將波動周期為一年的年間波動組分表示為作物的固有生長趨勢。去除這些因素的影響，剩余的年際波動組分和殘差可能含有重金屬脅迫信息。本文選取皮爾遜相關系數衡量年際信號之間的相似度。

1.3.4 動態(tài)時間規(guī)整（DTW）

DTW支持時間軸的伸縮和彎曲，能夠對不同長度的時間序列進行相似性度量[33]。它是用滿足一定條件的時間規(guī)整函數描述兩個時間序列之間的對應關系，求出兩個時間序列累加距離最小時對應的規(guī)整函數。其路徑的選擇滿足3個條件：邊界條件、連續(xù)性和單調性。DTW通過尋找到最小的扭曲路徑代價來對齊不同的時間序列。

本文引入DTW方法計算時間序列曲線的距離大小，因為水稻在不同年份、不同地域條件下，其物候期也存在一定的差異，而DTW能夠有效地解決這一類問題。通過計算分解后的時間序列與健康水稻時間序列之間的DTW距離，可以監(jiān)測水稻重金屬脅迫水平。

1.3.5 重金屬脅迫程度的提取

本研究利用WOFOST模型模擬出健康水稻的LAI時間序列，經過EEMD方法分解提取的無脅迫信號作為參考序列，原始LAI時間序列提取的含有重金屬脅迫信號作為查詢序列，計算兩者之間的DTW距離，即得到脅迫指數（Stress Index，SI）。

式中SIi表示水稻像元i的重金屬脅迫指數，Si為水稻像元i當年的原始LAI時間序列分解提取的可能含重金屬脅迫序列，Sh代表當年健康水稻LAI分解提取的無脅迫信號，Ddtw表示兩者之間的DTW距離。

為了更加直觀地了解水稻重金屬脅迫的程度，本文對脅迫指數SI進行歸一化處理，歸一化的公式如下：

式中SIn表示歸一化脅迫指數，SIi0.05和SIi0.95分別表示脅迫指數5%和95%分位數，SIn的值范圍為[0,1]，越接近1代表該水稻像元i受到的重金屬脅迫程度越高。

在本文中，采用標準方差分析[10]將水稻重金屬脅迫程度劃分為輕度、中度和重度3個等級，其SIn值分別介于0～0.40、＞0.40～0.75、＞0.75～1之間。

2 結果與分析

2.1 水稻LAI的分布

通過對Sentinel-2數據的LAI進行反演，得出試驗區(qū)內所有像元的LAI值，根據隨機森林方法提取水稻區(qū)域進而疊加提取水稻像元LAI值（圖3）。圖3中白色表示非水稻的區(qū)域，紅色代表更低的LAI值，綠色代表更高的LAI值。 從圖中可以看出試驗區(qū)西北和東南區(qū)域的LAI值相對較低，偏東北區(qū)域的LAI值相對較大，且分布較為密集。

2.2 EEMD信號分解

基于EEMD將LAI時間序列分解為7個不同波動頻率的IMF分量和1個殘差項，通過不同分量的合成可以得出年內波動組分、年間波動組分以及年際波動組分（圖 4）。圖4中的IMF1、IMF2、IMF3和IMF4的波動周期小于一年，該4個分量之和表示年內波動組分；而IMF5分量的波動周期與原始LAI時間序列的波動周期是具有一致性的，因此將IMF5分量對應于年間波動組分；剩余的IMF6、IMF7對應于年際波動組分，波動時間周期大于一年。

年內波動組分包含的信息具有較高的時間分辨率，時間跨越尺度小于一年，信號波動比較雜亂，其中不僅包含其他短期脅迫還包含自然環(huán)境中的噪聲等信息。可見IMF1～4分量的波動周期均小于一年，環(huán)境中的噪聲信息雖然長期存在并對信號產生一定的波動影響，但是噪聲信號沒有一定的規(guī)律。

通過計算相鄰年份年內波動組分的皮爾遜相關系數可知，2017年和2018年的相關系數為0.13，2018年和2019年的相關系數為0.33，表明相鄰年份年內波動組分之間的相關性很弱，因此，能更好地說明年內波動組分是表示短期脅迫和噪聲等不穩(wěn)定的信號。

IMF5分量的波動周期為一年，與原始水稻LAI具有相同的波動規(guī)律，因此這個分量能夠表示水稻本身固有的生長趨勢。計算年間波動組分與健康水稻LAI之間的相關系數，2017年、2018年和2019年與健康水稻之間的相關系數分別為0.96、0.85、0.90，均表現(xiàn)極強的相關性。說明年間波動組分揭示了水稻固有的生長趨勢。

原始LAI時間序列由持續(xù)的環(huán)境脅迫、作物的固有生長趨勢和其他脅迫信號共同組成，通過對上述兩種組分的剔除，重金屬脅迫信息存在于剩余的分量中，即年際波動組分與殘差之和（IMF6、IMF7和殘差）。

2.3 水稻重金屬脅迫程度

2017—2019年水稻重金屬脅迫的空間分布特征如圖 5所示，研究結果表明，不同年份的水稻重金屬脅迫的空間分布具有很大的相似性，在空間上的分布呈現(xiàn)一種穩(wěn)定的特征。受重金屬脅迫程度較低的區(qū)域大部分位于試驗區(qū)域的中部地區(qū)。而重金屬脅迫程度 較高的地區(qū)主要是集中在試驗區(qū)域西部、東北部以及偏東南地區(qū)。湘江東側和西側的工廠較為密集，重金屬脅迫的空間分布與工廠空間分布具有相似的趨勢，可能是由于早期的污水灌溉導致土壤重金屬含量超標。

進一步對試驗區(qū)域重金屬脅迫程度進行統(tǒng)計分析（表1），由表1可知，2017—2019年輕度重金屬脅迫分布面積比例較大，不同年份間其所占比例相差較小，分別為38%、44%和45%。重度重金屬脅迫分布面積比例約占29%，且不同年份的分布面積比例幾乎一致，2017年與2018年僅相差0.17個百分點，2018年與2019年相差不到1個百分點。從3 a的重金屬脅迫面積統(tǒng)計數據可知土壤中重金屬含量存在時空穩(wěn)定性。

表1 研究區(qū)不同程度重金屬脅迫水稻面積比例Table 1 Percentage of rice area under different heavy metal stress in the study area %

綜上所述，通過計算多年的歸一化脅迫指數SIn可以準確確定不同水稻重金屬脅迫水平的分布情況。

2.4 精度評價

研究利用野外實測的土壤重金屬Cd的含量數據對水稻重金屬脅迫監(jiān)測的結果進行精度評價，通過計算實測的土壤Cd含量與歸一化脅迫指數之間的皮爾遜相關系數進行驗證。兩者之間的相關系數高達0.851（圖6），由此說明歸一化脅迫指數SIn對水稻重金屬脅迫反應比較敏感。

3 討 論

本文借助EEMD方法進行脅迫信號分解，剝離其他脅迫以及噪聲等因素對重金屬脅迫的影響，可以更加準確地研究試驗區(qū)域的水稻重金屬脅迫。由于重金屬脅迫在土壤中可移性小，難以降解，因而土壤重金屬污染具有長期性[34-35]。有研究基于長時序的LAI結合EEMD方法提取水稻重金屬脅迫特征，并進行脅迫特征的穩(wěn)定性分析[8]。目前，越來越多的研究將DTW方法應用于遙感研究領域中，利用動態(tài)時間規(guī)整算法計算得到的距離作為相似性度量的指標來進行物種的識別分類[35]，也有研究通過DTW方法構建指標分析重金屬脅迫[9]，但是該方法沒有考慮水稻同時受到多種脅迫共同影響的情況，因此建立的模型較為簡單，在實際的應用中會存在局限性。因此如何在復雜的生態(tài)環(huán)境中研究重金屬脅迫是一個重要的課題。本研究利用EEMD與DTW相結合的方法先將其他脅迫與噪聲等因素排除，計算與健康水稻之間的距離判別水稻受重金屬污染的程度，在時間和空間上表明重金屬脅迫具有穩(wěn)定的特性，相鄰年份之間的DTW距離體現(xiàn)年際變化特征，可以作為衡量重金屬脅迫穩(wěn)定性的一個重要指標。

本文選擇利用DTW方法研究重金屬脅迫的一個重要的原因是DTW作為一種相似性度量的方法，它可以對時間序列進行拉伸和縮放，將兩個時間序列進行匹配比較，從而使識別效果更佳。湖南由于種植一季稻，插秧時間等的不同會導致水稻不同年份之間物候期有所差異，因此引入DTW方法可以很好地消除物候期不同導致的距離增加，從而提高水稻重金屬脅迫研究的精確度。

SIn是在無其他脅迫和環(huán)境影響因素的情況下計算得出的歸一化脅迫指數，根據其值的大小可以判斷出水稻重金屬污染的程度。試驗結果表明重金屬分布的范圍在3 a間變化較小，在空間上表現(xiàn)出一定的穩(wěn)定特征。因此，引入時序穩(wěn)定檢測方法，對水稻重金屬脅迫進行遙感監(jiān)測具有很高的可信度。

本研究僅采用3 a的數據對水稻重金屬脅迫進行研究，今后可以通過更長時序遙感數據對其進行分析；另外，本研究利用土壤重金屬Cd的含量來驗證水稻重金屬脅迫水平，在今后的研究中將進一步通過測定水稻中的重金屬含量驗證研究結果。

4 結 論

本文利用Sentinel-2數據同化作物生長模型獲取了葉面積指數（Leaf Area Index，LAI）時序數據，利用集合經驗模態(tài)分解（Ensemble Empirical Mode Decomposition，EEMD）方法分解LAI時間序列，提取出含有重金屬脅迫的信號，利用動態(tài)時間規(guī)整（Dynamic Time Warping，DTW）方法計算其脅迫信號與健康水稻LAI之間的距離大小，監(jiān)測水稻重金屬脅迫水平。結果表明：

1）原始LAI時間序列由持續(xù)的環(huán)境脅迫、作物的固有生長趨勢和其他脅迫信號共同組成，經EEMD分解為7個不同的分量和1個殘差項，通過分析得出重金屬脅迫信號存在于IMF6、IMF7和殘差的組合分量中；

2）歸一化脅迫指數是水稻重金屬脅迫敏感的參數，與土壤重金屬含量的相關系數為0.851，水稻受到的脅迫程度越高，歸一化脅迫指數值越大，反之越低；

3）在試驗區(qū)中，水稻重度重金屬脅迫的分布面積比例相對較低，且主要集中在西部、東北部以及偏東南地區(qū)，水稻重金屬脅迫具有一定的時空穩(wěn)定性。