不同降噪方式下基于高分五號影像的土壤有機質反演

劉煥軍，鮑依臨，孟祥添，崔 楊，張艾明，劉云超，王丹丹

（1. 赤峰學院資源環境與建筑工程學院，赤峰 024000；2. 東北農業大學公共管理與法學院，哈爾濱 150030；3. 中國科學院東北地理與農業生態研究所，長春 130012）

0 引 言

土壤有機質（Soil Organic Matter，SOM）含量是土壤肥力和土壤安全的重要指標，實現快速、精確的SOM空間制圖可以為精準農業和區域土地可持續發展提供監測依據[1-2]。高光譜衛星系統被設計成通過電磁光譜的可見光-近紅外區域，將表面光學性質分離成數百個光譜分辨率小于10 nm的光譜[3]，根據其在400～2 500 nm波長范圍內的反射率，可以為預測土壤的物理、化學和生物學特性提供解決方案[4-5]。然而，衛星傳感器的周期性偏移，載荷元器件間的電磁干擾以及波段間相互作用產生的噪聲降低了圖像的質量，甚至掩蓋數字圖像中真正的輻射信息[6-7]；同時，相比于植被、道路、建成區等地物，裸土的反射率低，光譜更易受到噪聲影響，因此，對高光譜衛星數據的校正及噪聲過濾是獲取精準土壤光譜反射率的前提與保證[8]。

高光譜圖像具有高分辨率和光譜間高相關性的特性，基于低秩的降噪方法成為研究熱點之一[9]。鑒于噪聲同時存在于空間域和光譜域，孟令軍[10]提出基于高光譜圖像像元光譜間相關性的光譜域降噪算法，分類精度提高了近8%。然而，其研究數據為合成高光譜圖像，包含Cuprite Nevada的部分數據和Indian Pines，這種數據類型并不常見，其結果在指導類似研究時具有一定的限制。黃冬梅等[11]提出基于多云協同的遙感圖像降噪方案，將遙感圖像共享給多個云服務器，云服務器根據改進的非局部均值降噪算法在密文域中完成圖像降噪。胡立栓[12]基于信息量和類間可分性的評價標準，分別提出了最優指數因子法與二進制粒子群算法，設計了基于該算法的特征選擇和分類器參數同步尋優策略。證實了降噪算法在降低處理成本、提升分類精度的有效性。Meng等[13]以反射率一階微分為基礎，結合離散小波算法對高分五號（GF-5）高光譜影像進行降噪處理，使土壤有機碳的預測精度提升了約10%。現有的降噪方法可分為空間域、光譜域和空譜結合 3種。基于信號頻率成分的奇異值分解方法隸屬于光譜域降噪方法；離散小波變換可以分離出高頻噪聲信息并保留光譜中的基本成分，屬于空譜結合的降噪方式[14]。中值濾波是基于排序統計理論的非線性信號處理技術，屬于典型的空間域降噪方法，尤其適用于對椒鹽噪聲的處理[15-16]。

上述降噪方法在圖像及信號處理中得到了廣泛的應用，然而，大部分研究僅應用于低維（圖片、音頻、多光譜影像）數據類型中。GF-5高光譜衛星影像是國內光譜分辨率最高的衛星，也是國際上首次實現對大氣和陸地進行綜合觀測的全譜段高光譜衛星，它由數以百計的波段信息組成，高維的空間信息中夾雜著較多的噪聲；此外，同類像素間的光譜反射率具有較強的一致性，導致同類地物之間的光譜特征相似。由于黑龍江省西部明水縣位于黑土帶中心處，SOM含量空間差異明顯，更具有代表性，因此本文選取明水縣為研究區，利用多種降噪方式對GF-5影像進行降噪處理，評估不同方式下反射率與SOM含量的相關性，提取二維光譜指數，運用降維方法對輸入量進行篩選并利用隨機森林算法構建SOM反演模型，以提升 SOM預測精度，并實現高精度的SOM土壤制圖，以期為實時定量監測土壤肥力變化提供依據。

1 研究區概況及數據來源

1.1 研究區概況

明水縣（124°18′～125°21′E，46°44′～47°29′N）位于黑龍江省西南部，松嫩平原西北部，地處高緯度地帶，平均海拔249.2 m，地勢由中部向東西兩側逐漸降低，屬中溫帶亞濕潤氣候，年平均氣溫約2.9 ℃，年平均降水量約480 mm。明水縣主要包括黑土、黑鈣土、風沙土，其中黑土和黑鈣土的面積占總面積 60%以上。黑中含有較多腐殖質，SOM含量高；風沙土中黏粒含量較少，保水保肥能力較差，SOM含量較低。縣域總面積為2 400 km2，其中耕地面積為1 087 km2，草原面積為467 km2，林地面積為360 km2。圖1為明水縣GF-5高光譜影像及采樣點空間位置分布。

圖1 研究區位置及樣點布置Fig.1 Location of study fields and sampling point distribution

1.2 數據獲取與處理

1.2.1 土壤數據獲取

2016年5月，在研究區范圍內0～20 cm深度選取采樣點共38個（圖1），利用GPS記錄采樣點的空間坐標。主要選取原則為沿公路兩側，兼顧土壤類型及物理性狀選取開闊平坦的非均質性區域作為樣地，樣本間隔約10 km。每個樣本是1個復合樣本，由5～6個子樣本組成，子樣本在30 m×30 m大小的網格范圍內隨機選取，以確保SOM含量能夠反映該區域平均水平。將獲取的土壤樣本置于布制土壤袋中存放，帶回實驗室進行風干處理，并去除結石，雜草根和其他雜質，研磨、干燥篩分至≤2 mm，采用重鉻酸鉀法測定樣品的 SOM 含量。分析38個采樣點 SOM 含量數據，結果表明，SOM的質量分數介于3.45%～5.53%之間，平均值為4.44%，樣本標準差為4.28%，變異系數9.6。

1.2.2 高光譜衛星影像數據獲取及預處理

研究區耕地范圍內采用傳統播種方式，作物生長期一般為一年一熟制，從 4月中旬陸續進行播種，此時土壤直接暴露于表面，地表無作物植被及秸稈殘留，無積雪覆蓋，這一時期被定義為“裸土窗口期”[17]。獲取裸土時期（2019年 4月 24日）GF-5影像（http://data.cresda.com:90/#/home）。影像獲取日期與土壤樣本采集日期接近，SOM含量不會發生明顯變化。GF-5影像數據的空間分辨率為30 m，地面覆蓋寬度60 km，包含 330個波段，其中在可見光近紅外（Visible and Near-infrared，VNIR）范圍內共包含150個波段，光譜分辨率為4.28 nm；在短波紅外（Shortwave Infrared，SWIR）范圍內共包含180個波段，光譜分辨率為8.42 nm。由于傳感器在400～430、900～1 050、2 450～2 500 nm波譜范圍內信噪比較低，且在1 350～1 451、1 771～1 982 nm范圍內受到大氣中水蒸氣吸收的影響，造成光譜數據不連續。因此，本文選取430～900、1 050～1 350、1 451～1 771、1 982～2 450 nm作為本次研究的光譜區間。使用ENVI 5.1輻射校準工具和GF-5數據頭文件中的增益和偏差比對影像數據進行輻射校準，利用FLAASH大氣校正模型對遙感影像進行大氣校正。

2 研究方法

2.1 降噪方法

2.1.1 奇異值分解

奇異值分解（Singular Value Decomposition，SVD）可以被看作一個能夠反映矩陣信息的矩陣代表值。利用奇異值分解的方法確定光譜反射信號的Hankle矩陣，對矩陣進行奇異值分解，構造逼近矩陣對含噪干涉信號進行降噪處理，奇異值越大時，它代表的信息越多，這一原理與主成分分析相似。假設Hm是m×z階矩陣，其元素全部屬于域Hm，則存在一個分解使得Hm=U∑V＊，式中U是m×m階酉矩陣；V＊即V的共軛轉置，是z×z階酉矩陣；Σ是半正定m×z階對角矩陣。當信號中存在噪聲時，Hm為滿秩矩陣，使該分解的奇異值之間存在如下關系：

式（1）中σk是SVD獲得的第k個奇異值，為噪聲信號中的分界點。因此，k個奇異值包含信號能量，隨后的奇異值對應于信號中的噪聲分量[18]，奇異值越大所包含的信息越多，取k前面若干個最大的奇異值和對應的左右奇異向量，可以基本還原出數據本身，因此可用于數據降維與壓縮。該降噪方式在Matlab2016a中實現。

2.1.2 離散小波分解

離散小波分解（Discrete Wavelet Transform，DWT）是以 2的整數次冪為底對連續小波的尺度和平移參數采樣，包含了信號分解與信號重構 2個過程。依據信號的長度和小波基長度將原始信號分解為高頻信號和低頻信號，其中高頻信號主要由噪聲信息構成，低頻信號就是小波系數，其后的每一步均分解上一層的低頻信號，迭代進行直至分解的最大尺度。在每個分解層j，建立 1個第j層的高頻信號（Dj），因此原始高光譜信號（S）就等于最終分解層的低頻信號（Aj）與所有分解層Dj之和，即

小波重構過程是將每層的低頻信息進行重構，從而得到各層的特征光譜。常用于對高頻噪聲的過濾，且對于光譜曲線中的“毛刺”部分的平滑效果顯著[13]。本研究在 Matlab2016a中測試了大量的母函數，最終選擇‘db4’小波母函數。

2.1.3 中值濾波

中值濾波（Median Filtering，MF）是基于排序統計理論的一種廣泛應用的平滑噪聲和保持圖像邊緣的非線性濾波方法，其基本原理為：用像素點鄰域灰度值的中值代替該像素點的灰度值，使周圍的像素值接近真實值從而消除孤立的噪聲點。對于給定的圖像I(i,j)，中值濾波過程如式（3）所示：

式中(u,v)∈(1,2,…,h)×(1,2,…,l)，Imf(u,v)為處理后得到的圖像；h和l是圖像的寬度和高度，r和s代表窗口位移大小，W是方形窗口內的坐標集。該中值濾波法在取出脈沖噪聲、椒鹽噪聲的同時能保留圖像的邊緣細節，在光學測量條紋圖像的相位分析處理方法中有特殊作用，是經典的圖像平滑噪聲的方法。本研究在Matlab2016a中測試了不同窗口濾波器，最終選擇 5×5方形窗的濾波器對影像進行降噪處理。

2.2 光譜指數的構建

光譜指數是不同波段反射率的線性或非線性組合形式，光譜指數的構建有助于提高遙感數據預測SOM含量的精度[19]。本研究探討差值指數（Difference Index，DI）、比值指數（Ratio Index，RI）、歸一化指數（Normalization Index，NDI）[20-21]與SOM之間的關系，計算過程如下：

式中Ro和Rp為波段o和p獲取的光譜反射率。

在900～1 050 nm受到傳感器自身缺陷的影響、且在1 350～1 451、1 771～1 982 nm受到大氣中水蒸氣吸收的影響，光譜數據不連續，因此，SOM估計的最佳光譜波段從 430～900、1 050～1 350、1 451～1 771、1 982～2 450 nm共計237個波段內選擇。波段選擇在Matlab2016a中實現。

2.3 主成分分析

主成分分析（Principal Component Analysis，PCA）的實質是基于數據統計特征的多維正交線性變化分析，將原有的高度相關變量的數據集轉換為一組不相關變量，各變量間相互獨立，前幾個主成分信息可以最大限度地反映原始多個變量的信息[22]。在本研究中，保留了總光譜方法解釋量＞85%的主成分信息作為本研究的輸入量。這一過程是在Matlab2016a軟件中通過princomp函數實現的。

2.4 RF預測模型

RF（Random Forest，RF）預測模型是一種集成模型，預測變量集在每個分裂點中都是隨機限制的[23]，避免了計算中樹與樹之間的相關性問題，從而實現高精度的預測。RF的性能取決于以下3個參數的設置：生成樹的數量（ntree）、葉片最小數量（nodesize）和每個節點處用于分割節點的預測變量數（mtry）。Ramón等[24]使用更多的樹的數量（ntree）可以獲得更穩定的估計變量重要性的結果，因此使用ntree =500。mtry默認值是預測因子總數的1/3，隨后，使用迭代方法以最高預測精度確定最佳nodesize值。

2.5 SOM預測模型精度驗證

研究采用建模集與驗證集之比為 2:1的方式建立模型，建模集26個，驗證集12個，校準和驗證樣本在研究區域內空間分散。為了檢驗不同降噪方法的效果，模型的性能通過精度度量決定系數（R2）、均方根誤差（Root-Mean-Square-Error，RMSE）、測定值標準偏差與標準預測誤差的比值（Residual Predictive Deviation，RPD）進行評估。R2越大，RMSE越小，證明預測精度越高；當 1.5＜RPD＜2時表明模型可以對樣品含量進行粗略估測，當 RPD＞2時表明模型具有較好的定量預測能力[25]。

3 結果與分析

3.1 DWT最優分解尺度的確定

圖2為9個DWT分解層（L1～L9）的光譜反射率曲線。由L1～L4可見，DWT降低了光譜曲線中的噪聲，使曲線變得更加平滑，特別是在2 200～2 450 nm之間較為明顯。整體來看，L1～L4可以較好地保留原始光譜曲線的特征。在L5～L9中，DWT消除了光譜曲線中大部分的細節信息，使得不同SOM含量的土壤除反射率高低的差異外，無其他明顯差異，因此，過高的分解層數不利于光譜細節的保留以及SOM預測精度的提升。原始光譜曲線中，在2 200～2 450 nm處存在“小毛刺”形狀，這是由于原始反射率噪聲傳遞導致的，DWT可以將這些高頻信號進一步去除，減弱噪聲傳遞現象。將整個波段范圍反射率與SOM做相關分析，由于波段信息豐富，表1列出了相關系數最高的3個波段，這3個波段主要集中在600～700 nm之間。

圖2 不同離散小波分解（DWT）分解層數下不同SOM含量土壤反射光譜曲線Fig.2 Reflectance spectral curves of soils with different SOM contents under different Discrete Wavelet Transform (DWT) decomposition layers

表1 不同DWT分解層光譜反射率與SOM相關性最強的前3名Table 1 Top three based on correlation between spectral reflectance of different DWT decomposition layers and SOM

由表1可知，L1～L3均可提升光譜曲線與SOM含量的相關性，且當分解尺度為1層時，反射率與SOM含量的相關性最高。同時，這些相關性高的波段出現的位置與原始反射率（給出值）最接近，因此，將高光譜衛星影像進行1層DWT時，可以實現較高SOM預測精度，下文將選擇1層DWT進行SOM反演。

3.2 不同降噪方式下土壤光譜曲線與SOM含量相關性

土壤反射率隨SOM含量的增加而降低（圖3），且不同 SOM 含量的土壤反射光譜曲線形狀大體一致。與OR相比，不同的降噪方法均有效降低了光譜曲線中的噪聲，使曲線變得更加平滑。

圖3 不同降噪方法下的土壤反射光譜曲線Fig.3 Soil reflection spectral curves under different noise reduction methods

在600～1 300 nm范圍內，反射率隨著波長增加而增加，光譜曲線的形狀差異比較明顯，其中，在 900～1 050 nm間出現1個反射峰，這個反射峰是由于2個傳感器連接處的噪聲較大造成的。在1 450～1 700 nm范圍內，反射率隨著波長增加呈現先降低，后升高的趨勢，且能夠較好地保留光譜細節；在1 950～2 450 nm內，反射率隨著波長增加先升高后降低。對比來看，SVD及DWT均對光譜曲線的不同位置進行了平滑，有效地實現了對高光譜衛星影像的降噪；然而，通過 MF降噪下的光譜曲線可知，SOM 質量分數處于 3.5%～4.0%與＞4.0%～4.5 %的光譜曲線不符合反射率隨SOM含量增加而降低的規律，此外，該降噪方式縮減了不同SOM含量間光譜曲線的差異，不利于SOM含量的預測。

降噪方法可有效地提高反射率與 SOM 之間的相關性，如SVD及DWT（圖4）。在800 nm附近，不同降噪條件的相關性由高到低依次為DWT、SVD、OR、MF，與 SOM 含量的相關系數依次為 0.625、0.557、0.581、0.481；在1 000～1 750 nm內，不同降噪方法的差別進一步擴大；在2 000～2 450 nm內，SVD、OR與SOM之間的相關系數大致相等。

圖4 不同降噪方法下光譜反射率與SOM含量的相關性Fig.4 Correlation between spectral reflectance and SOM content under different noise reduction methods

3.3 光譜指數的選擇

圖5展示了不同降噪處理方式下，所選取的光譜指數（DI，RI，NDI）與SOM之間的相關關系。將與SOM相關性最高的波段組合列于表2，整體而言，與OR相比，基于DWT降噪方式下的光譜指數的相關性有所提升，SVD方式下僅RI相關性有所提升，MF處理下的相關性均低于OR反射率的相關性，不具備SOM反演的優勢。

圖5 土壤有機質含量與二維光譜指數相關圖Fig.5 Correlation diagram of SOM content and 2-D spectral index

表2 土壤有機質含量與二維光譜指數相關性Table 2 Correlation between SOM and 2-D spectral index

3.4 不同降噪方式下的SOM預測結果

每種降噪方式包含237波段信息及3個二維光譜指數，共計240個輸入變量。利用PCA方法提取貢獻率大于 85%的信息進行預測，各種降噪處理下的主成分信息見表3。由表3可知，不同降噪處理下，累計貢獻率均高于95%。可見，不同降噪處理后利用PCA進行波段降維，仍能夠保留約95%的原始信息。

表3 基于主成分分析法的輸入量統計表Table 3 Statistics of inputs based on principal component analysis method

建模集SOM預測精度R2為0.62～0.80，RMSE為1.80%～2.35%，RPD為1.63～2.23，其中以DWT降噪方式下模型模擬精度最高（圖6）。驗證集SOM的預測精度R2約為0.49～0.69，RMSE為2.26%～2.65%，RPD為1.41～1.80。結合 1:1線比較 SOM 的預測值與實測值，當以OR和SVD作為輸入量時，SOM預測值分布比較集中；以MF作為輸入量時，無論SOM含量處于較低或較高時，其預測值與實測值相差較大，說明該方式未能對數據進行有效的降噪。而使用DWT降噪方法時，其SOM預測值擬合曲線與1:1線更為接近，模擬精度更高，誤差更小，且DWT方法與其他降噪方式相比具有更好的降噪效果，能夠實現較高的SOM反演精度。

3.5 基于GF-5高光譜衛星數據反演圖

明水縣土壤分布以黑土為主，土壤肥力較高。中部黑鈣土形成于低地貌環境，受低平原地下水控制，有利于鈣質成分聚集，而形成鈣積層，表層土被有充分的成壤條件，使腐殖質能夠保存并得以逐漸加厚。西部地區地勢較低，以草甸土為主，草甸土通過土壤侵蝕等作用沉積了較多的土壤養分（圖7）。整體而言，中部漫川漫崗地形，SOM含量高且相對均一；東南部侵蝕嚴重，影響了SOM含量的積累。

圖6 基于RF預測模型的SOM預測結果Fig.6 Results of SOM inversion based on RF model

圖7 基于GF-5的RF模型明水縣SOM分布Fig.7 Distribution of GF-5-based SOM predicted by RF model in Mingshui County

4 討 論

高光譜衛星影像中包含豐富的波段信息，但由于受多種因素影響，影像中包含著大量噪聲信息，不利于進行SOM的反演研究。因此，本文選取常見的降噪方法對高光譜衛星數據進行處理，同時，判斷降噪后的影像信息是否可以提升SOM預測的精度。本文從空間域、光譜域和空譜結合 3種降噪方法出發，分別選取不同降噪原理的方式進行SOM預測。SVD方法通過保留影像中權重大的矩陣信息，實現光譜數據的降維；DWT是基于不同的分解尺度對高頻信息進行過濾[26]；MF則采用象元平均化的方式以“抹平”噪聲。相比于土壤本身的反射率，影像中的噪聲更偏重于高頻噪聲。SVD和DWT降噪方法可以消除光譜數據中的冗余信息，增大不同SOM含量之間光譜曲線的差異，有助于實現高精度土壤屬性的預測。對比來看，DWT將光譜曲線進行不同尺度的分解，得到不同尺度下土壤的光譜特征，使光譜與SOM含量的相關性得以提升。相比之下，MF這種象元平均化的方式未能有效的消除土壤光譜中的噪聲，且降噪后的光譜與原始反射率相比，其與SOM含量的相關性有所降低。

利用高光譜遙感進行SOM預測研究的數據降維主要通過兩方面實現：輸入量和模型[27-28]。偏最小二乘回歸模型適合處理自變量大于樣點數的數據，適用于服從于正態分布的數據集，因此不具備普適性。這項研究中，每種降噪方式下均包含大量的波段信息，致使數據冗余，不利于模型的計算及分析。因此，本文采用PCA方法，將高維信息投射到更小維度的子空間中，保留原始數據大部分的信息，以縮減計算機運行時間，提升SOM預測精度[29]。本文利用不同的降噪方法對GF-5高光譜衛星數據進行降噪處理，為降低高光譜衛星影像中的噪聲提供思路，然而，研究未能完全識別并消除影像中存在的噪聲，在今后的研究中，應該與室內實測數據結合，充分校正數據中潛在的誤差信息，將其用作開發其他領域（例如植被和水體）的參考方法。

5 結 論

考慮到高光譜衛星數據中的噪聲不利于獲取真實的地物反射光譜曲線，本文基于空間域、光譜域和空譜結合3種降噪方法，篩選用于GF-5高光譜數據的最佳降噪方式，主要結論如下：

1）奇異值分解（Singular Value Decomposition，SVD），離散小波變換（Discrete Wavelet Transform，DWT）能夠有效降低高光譜數據中的噪聲，與未降噪處理的光譜曲線相比，這 2種方法有效地減少了曲線中的“毛刺”部分，使光譜曲線更加平滑。

2）光譜指數可以提升波段信息與土壤有機質（Soil Organic Matter，SOM）的相關性，不同降噪方法對于SOM預測的潛力由高到低依次為 DWT、SVD、原始反射率（Original Reflectance，OR）、中值濾波（Median Filtering，MF），與 SOM 含量的相關性由高到低依次為 0.625、0.557、0.581、0.481。

3）降噪后，采用組成分分析降維，基于隨機森林構建SOM模型。比較各降噪方式，最佳降噪方法為DWT，其驗證集R2為0.69，RMSE為2.26%，RPD為1.80。基于DWT降噪方式，采用隨機森林模型預測的明水縣SOM分布與實際相符，進一步證實了降噪方法應用于處理高光譜衛星數據的可行性，在未來的研究中，建議使用室內實測光譜曲線作為高光譜衛星影像方法的補充，以更好地模擬不同光譜區間的噪聲特征，拓展高光譜衛星在其他領域中的應用。