丁松滔，張霞，尚坤，李儒，孫偉超

1.中國科學院空天信息創(chuàng)新研究院,北京 100101;2.中國科學院大學,北京 100049;3.自然資源部國土衛(wèi)星遙感應用中心,北京 100048

1 引言

成像高光譜遙感能夠獲取圖像上連續(xù)且精細的土壤反射光譜，具有低成本大范圍快速監(jiān)測土壤重金屬含量的潛力（Khosravi 等，2018）。然而目前大多數(shù)研究采用實驗室理想條件下采集的光譜，基于高光譜圖像的重金屬反演研究較少（Ding 等，2022），并且由于土壤樣本采集受制于地形、交通可達性等因素，且土壤樣本的重金屬實驗室化驗分析成本高，導致樣本不充足的問題普遍存在，而對于高光譜圖像而言，圖像上的土壤像元數(shù)量與土壤樣本數(shù)相差懸殊，小樣本問題更加突出。光譜特征選擇是解決小樣本問題的有效途徑之一，而服務于光譜特征選擇的特征增強方法也是研究熱點，研究顯示經(jīng)過特征增強的光譜曲線比原始光譜曲線可以更有效地提取出光譜的特征波段（沈強 等，2019）。

微分處理是常用的增強光譜特征的方法，一階微分和二階微分是最常用的光譜微分方式，能夠突出土壤光譜中與重金屬相關的光譜信息，提高反演精度（郭學飛 等，2020）。然而，傳統(tǒng)的整數(shù)階微分（即一階、二階微分）缺乏對可能包含土壤重金屬有益信息的漸變傾斜度或曲率的敏感性（Hong等，2018）。分數(shù)階微分作為整數(shù)階微分的擴展，對光譜微分的處理更加細致，階數(shù)上更精細的量化意味著能更敏感地捕捉各波段的曲率和傾斜度變化，突出光譜特征。已有學者將FOD應用于實驗室光譜（Chen 等，2022），成功地反演了土壤重金屬Cr、Pb、Zn的含量；Meng等（2021）將FOD應用于GF-5圖像，使用FOD處理后的像元光譜計算光譜指數(shù)作為自變量，構建的有機質反演模型精度R2達到0.84，但FOD 應用于高光譜圖像上進行土壤重金屬反演的有效性還有待驗證。

特征選擇算法通過選出特征性較強的波段，提高模型的反演性能?，F(xiàn)有很多研究都已經(jīng)證明，使用優(yōu)選出的特征波段比使用全波段可以獲得更好的反演結果（Tan 等，2020；Zhang 等，2021），這表明波段選擇的重要性。土壤重金屬反演中最為常見的篩選方法是皮爾森相關系數(shù)法，但該方法只考慮單個波段與理化性質間的關系（周冰 等，2021），沒有考慮多個波段間的共線性。GA 是一種以“適者生存”為原則的優(yōu)化算法，以最大化精度為目標，在特征空間中進行啟發(fā)式搜索，可以快速獲取近似最優(yōu)解（柏晗 等，2022），已有將其應用于遙感圖像上反演重金屬含量的先例（Wang等，2022），然而GA的缺點是會過早收斂，導致結果難以跳出局部解（Pavez-Lazo 和Soto-Cartes，2011）。CARS 同樣是遵循“適者生存”的一種篩選算法（Li 等，2009），通過逐步去除冗余和不重要變量來選擇信息，使模型計算效率與穩(wěn)定性提升（Vohland 等，2017），能夠從全波段數(shù)據(jù)中選擇最優(yōu)的波段組合（Vohland 等，2016），Cheng 等（2021）在土壤全氮反演研究中指出CARS 算法在有效變量的選擇上優(yōu)于GA。綜上所述，CARS 算法在這3 種算法中有利于篩選出更有效的波段組合，達到全局最優(yōu)解。

本文以新疆維吾爾自治區(qū)哈密市的黃山南礦區(qū)為研究區(qū)，應用航空高光譜圖像開展對土壤重金屬Pb、Zn、Ni 含量的反演研究，提出了一種基于FOD 估算高光譜圖像上土壤重金屬濃度的反演方法，通過擴充樣本、FOD 增強光譜特征、CARS篩選波段選取優(yōu)光譜特征，探索高光譜圖像反演中小樣本問題的解決途徑，分析確定了研究區(qū)3種重金屬反演的最佳FOD 階數(shù)，并將CARS 與CC、GA 的建模結果進行對比；進一步分析了CARS 的波段篩選結果，得出對3種重金屬反演貢獻最高的波段范圍及其物理意義，分析了重金屬含量分布圖的可靠性。

2 研究區(qū)與研究數(shù)據(jù)

2.1 研究區(qū)概況

研究區(qū)位置如圖1所示，位于新疆維吾爾自治區(qū)哈密市黃山南銅鎳礦區(qū)，是重要的銅、鎳、鐵、鉛、鋅等大型礦床集中區(qū)（王京彬 等，2006），該礦區(qū)于2007年開始勘探，至今仍在開采中。地處新疆維吾爾自治區(qū)東部，屬于典型的溫帶大陸性干旱氣候，干燥少雨，年降水量33.8 mm，年蒸發(fā)量3300 mm，夏季酷熱、蒸發(fā)強。研究區(qū)內(nèi)主要為空曠戈壁，無植被覆蓋，土壤表面多為裸露狀態(tài)。

圖1 研究區(qū)位置示意圖Fig.1 The geographical location of the study area

2.2 數(shù)據(jù)獲取及預處理

航空高光譜圖像獲取于2021 年8 月23 號，傳感器為中國科學院上海技術物理研究所研制的全譜段多模態(tài)成像光譜儀，本次航飛航高3 km，波長范圍覆蓋350—2500 nm，350—1000 nm 光譜區(qū)間的空間分辨率為0.5 m，有251 個波段，1000—2500 nm 波長范圍空間分辨率為1 m，有508 個波段。研究區(qū)由兩條航帶所覆蓋，圖像經(jīng)過輻射定標、大氣校正后得到地表反射率數(shù)據(jù)，對兩航帶圖像進行配準、拼接等處理后，裁剪出研究區(qū)，經(jīng)緯度范圍為：42.197°N—42.216°N，94.625°E—94.684°E。

成像光譜儀在獲取光譜數(shù)據(jù)時，大氣中的水汽對1400 nm 和1900 nm 波長附近的輻射能量存在強吸收，且1900—2500nm 的像元光譜存在較嚴重噪聲。為去除水汽等噪聲的影響，并保留盡可能多的波段，剔除1350—1450 nm 和1800 nm 之后的波段。對剩余的光譜進行SG 濾波去噪處理，由于不同區(qū)間的光譜噪聲不同，對其采用不同的窗口大小進行濾波，350—1000 nm 區(qū)間的噪聲較弱，采用窗口大小為9 的二次多項式，1000—1800 nm區(qū)間為中度噪聲，采用窗口大小為13 的二次多項式。

與航空飛行準同步，2021年7月22日—7月30日開展了土壤樣本地面獲取實驗，由于研究區(qū)屬于干旱區(qū)，且該時段內(nèi)土壤未受降雨影響，兩者獲取數(shù)據(jù)時土壤狀況可認為一致。土壤樣本采集主要針對采礦區(qū)及周邊，沿道路兩旁選擇土壤的顏色和粒徑有明顯差異的樣點，采集0—20 cm 表層土，共采集土壤樣本72 個。圖像及土壤樣本點位置見圖2，根據(jù)土壤樣本的經(jīng)緯度信息，提取對應的圖像像元反射光譜數(shù)據(jù)。

圖2 航空遙感圖像及采樣點位置Fig.2 Aerial remote sensing image and sampling locations

將土樣在實驗室風干研磨后過100目篩制成標準樣，用火焰原子吸收分光光度法測定重金屬含量。圖3 所示為樣本中3 種重金屬含量的直方圖，可見重金屬Pb 和Zn 的含量分布都趨向于正態(tài)分布，而重金屬Ni 含量呈現(xiàn)出明顯的偏分布，本文采用對數(shù)變換對其進行校正，轉換后Ni 含量分布如圖4所示，偏分布現(xiàn)象得到明顯改善。

圖3 土壤樣本的重金屬含量直方圖Fig.3 Histogram of the heavy metal content of soil samples

圖4 對數(shù)變化后Ni含量直方圖Fig.4 Histogram of the Ni content after log change

3 實驗方法

圖5為本文研究技術路線圖，首先根據(jù)土壤樣本的經(jīng)緯度提取對應像元和相似的鄰近像元光譜，對像元光譜數(shù)據(jù)進行SG 濾波及分數(shù)階微分處理，處理后的光譜采用CARS 算法進行特征波段優(yōu)選，選出的波段組合用于偏最小二乘回歸PLSR（Partial Least Squares Regression）建立反演模型，最終將構建的最優(yōu)模型應用于航空高光譜圖像上反演重金屬含量，獲得重金屬含量分布圖。

圖5 技術路線圖Fig.5 Flowchart of the proposed method

3.1 樣本擴充

由于研究區(qū)的地形較平坦，地表一致性較好，可以認為所采集土壤樣本的重金屬含量可以代表該樣本小范圍內(nèi)的土壤重金屬含量，又因為本研究所獲取的高光譜圖像空間分辨率為1 m，可認為鄰近像元所對應的土壤中重金屬含量差異極小。因此，本研究根據(jù)樣本經(jīng)緯度信息定位到圖像上的像元后，以該像元為中心像元，提取中心像元以及其八鄰域的像元光譜，依次計算八個鄰近像元光譜與中心像元光譜的歐氏距離，距離越近則認為其與中心像元的相似度越高，選出歐氏距離最小的兩個像元，將3個像元（中心像元和兩個相似鄰近像元）的光譜共同作為該土壤樣本的反射率光譜信息，達到將土壤樣本的光譜數(shù)據(jù)擴大3 倍的效果。本研究采集的72 個土壤樣本，用此方法進行樣本擴充后，獲得216條像元反射率光譜數(shù)據(jù)。

3.2 分數(shù)階微分（FOD）

分數(shù)階微分是數(shù)學中的重要分支，它將經(jīng)典的整數(shù)階微分推廣為任意階，能更敏感地捕捉光譜反射率細節(jié)的變化。到目前為止，還沒有一個統(tǒng)一的公式來定義分數(shù)階微分。數(shù)學家們從不同的角度分析，得出了分數(shù)階微積分的不同定義。目前，分數(shù)階微分常見的表達形式主要包括Riemann-Liouville（R-L）、Grunwald-Letnikov（G-L）和Caputo。其中最常用的形式為G-L（Wang等，2017），本研究采用G-L形式進行微分，其定義如下：

式中，α為任意階數(shù)；h為微分步長；t與α分別為微分的上、下限；Γ(α)為Gamma 函數(shù)，Γ(α)滿足：

令h=1，式（1）能夠推導出函數(shù)f(x)分數(shù)階微分的表達式為

3.3 競爭自適應重加權采樣法（CARS）

CARS 是一種以“適者生存”為原則的篩選算法，基于與模型性能相關的統(tǒng)計數(shù)據(jù)，通過逐步去除冗余和不重要的變量來選擇信息量大的變量，將其應用于高光譜數(shù)據(jù)進行篩選時，其本質是使用自適應加權采樣技術來保留回歸中絕對系數(shù)較大的光譜波段（Tan等，2021）。

圖6 為CARS 的方法流程圖，CARS 通過N次蒙特卡羅抽樣迭代地選擇N個波段子集，最終目的是選擇一個能使交叉驗證的均方根誤差（RMSEcv）最低的最佳變量子集。在高光譜數(shù)據(jù)上使用CARS篩選的具體步驟如下：

圖6 CARS算法流程圖Fig.6 Flowchart of CARS

（1）使用蒙特卡羅抽樣方法從數(shù)據(jù)集中隨機選擇具有固定比例的樣本，然后用這些樣本建立PLSR 模型。PLSR 模型可用公式表達為y=bX+e，e為常數(shù)項，b為回歸系數(shù)向量，b=［b1，b2，…，bn］，b中第i個元素的絕對值|b|（1≤i≤n）表示第i個波段對y的貢獻，波段對目標變量的貢獻越大，該波段就越重要。為評價每個波長的重要性，定義權重Wi：

通過CARS 算法去掉的變量，其權重Wi均設為0；

（2）利用指數(shù)遞減函數(shù)EDF（Exponentially Decreasing Function）和自適應重加權采樣ARS（Adaptive Reweighted Sampling）分別強行和競爭性地消除權重低的波段；

（3）重復步驟（1）—（2），直到達到N次采樣運行，最后選擇RMSEcv 最低的波段子集作為最優(yōu)波段組合。

由于蒙特卡洛抽樣每次迭代抽取固定比例的樣本參與運算，不使用所有樣本，因此選出的波段組合具有更好的適應性；步驟（2）中使用指數(shù)遞減函數(shù)去除波段，能在迭代前期篩選掉大量重要性低的波段，使篩選過程的運算量顯著下降。根據(jù)多次實驗結果，將蒙特卡羅抽樣次數(shù)設置為50次，每次迭代抽取90%的樣本用于運算。

3.4 模型構建與評價指標

PLSR 是土壤光譜分析中最為常用的反演模型，在處理共線性強、計算復雜度高的問題上具有天然優(yōu)勢（Rossel 和Behrens，2010），能夠應對自變量多于樣本個數(shù)的情況，因此其非常適合處理高光譜數(shù)據(jù)，并且已被廣泛的應用于土壤重金屬反演。本文采用CARS+PLSR 的模型構建方法，先通過CARS篩選出特征波段，再將特征波段輸入PLSR 建立反演模型。在PLSR 建模中，通過留一交叉驗證的最小均方根誤差確定PLSR 的潛變量個數(shù)。由于CARS的蒙特卡羅抽樣具有隨機性，為了避免實驗的偶然性影響，本文對每組反演實驗均進行5 次，以5 次反演實驗中的精度最優(yōu)值作為該模型精度將精度最優(yōu)的模型應用于高光譜圖像，得到土壤重金屬的含量分布圖。

精度評定采用預測均方根誤差RMSEP（Root Mean Square Error of Prediction）、相對分析誤差RPD（Ratio of Prediction to Deviation）和決定系數(shù)R2（Coefficient of Determination）3 個評價指標。3 個評價指標，RMSEP 值越小，RPD 值越大，R2值越接近1，說明反演模型的精度越高；反之，RMSEP值越大，RPD 值越小，R2值越小，說明反演模型的精度越低，本文模型優(yōu)劣參考現(xiàn)有的土壤屬性含量高光譜反演的評價標準（Wang 等，2014）：出色模型，R2>0.9；良好模型，0.9>R2>0.8；近似模型，0.8>R2>0.65；具有一定反演能力，0.65>R2>0.50；不具備反演能力，0.50>R2。

4 實驗結果與分析

反演中將216個樣本按照2∶1的比例劃分為訓練集和測試集。由于樣本擴充出的3個樣本重金屬含量相同，如果將樣本重金屬含量排序后采用分層抽樣方法，按2∶1進行劃分，會使訓練集和測試集的樣本含量分布完全一樣，從而導致反演精度偏高。為了解決該問題，選出具有代表性的訓練集，采用Kennard-Stone（KS）算法（Zhang 等，2017）來劃分樣本集，具體的實施步驟如下：

（1）首先計算兩兩樣本之間的光譜歐氏距離，選擇光譜距離最大的兩個樣本；

（2）然后分別計算剩余樣本與已選兩樣本之間的光譜歐氏距離；

（3）對于每個剩余樣本而言，計算其與已選各樣本之間的最短光譜距離，選擇這些最短光譜距離中相對最大的距離所對應的樣本，作為新入選的樣本；

（4）重復步驟（3），直至所選樣本的個數(shù)等于事先設定的數(shù)目為止。

通過KS 算法選出144 個樣本作為訓練集，剩下的72 個樣本作為測試集，樣本集劃分結果通過直方圖展示于圖7，3 種重金屬的訓練集、測試集和總樣本集的含量直方圖分布一致且有一定差異性，表明該方法劃分出的訓練集和測試集具有較好的代表性。

圖7 各樣本集的重金屬含量直方圖Fig.7 Histogram of heavy metal content for each sample set

4.1 分數(shù)階微分的最佳階數(shù)確定

由于分數(shù)階微分需使用相鄰的波段進行計算，理想條件下的土壤反射率光譜應是一條連續(xù)的光滑曲線，相鄰波段的反射率不應出現(xiàn)較大變化。本研究使用的航空高光譜圖像采集350—1000 nm和1000—2500 nm 波長范圍的傳感器不同，導致光譜在1000 nm 波長處的反射率曲線不夠平滑。此外，由于1350—1450 nm 的水汽吸收帶被剔除后，導致該光譜區(qū)間兩端波長相差100 nm 的兩個波段變?yōu)橄噜彶ǘ巍榱舜_保微分時相鄰波段反射率的差異在合理范圍內(nèi)，將光譜分為350—1000 nm、1000—1350 nm 和1450—1800 nm 這3 個區(qū)間來分別進行分數(shù)階微分處理。

從圖8中可以看出，隨著階數(shù)的增加，光譜的特征變的越來越明顯，峰谷差異越發(fā)增大；0—1階的微分結果顯示出，階數(shù)增大的過程中，平緩區(qū)間的波段反射率值越發(fā)趨近于0，而波峰波谷區(qū)間的波段反射率值被逐漸放大；1—2 階的微分結果中，反射率曲線的波形已較為相似（由于1—2階內(nèi)的微分曲線較為相似，圖8中僅展示1階、1.5階和2 階的微分結果），均與原始光譜曲線差異較大，特征放大的效果明顯，隨著階數(shù)的增加，變化主要體現(xiàn)光譜曲線的極值不斷增大。

圖8 不同階數(shù)微分后的像元光譜反射曲線Fig.8 Image spectral curves with different orders of differential

在不同階數(shù)的微分處理后，反射率曲線中的光譜信息有明顯差異，為了選出適合于反演各重金屬的最佳微分階數(shù)，在0—2階的范圍內(nèi)，以0.1階為間隔，將不同階數(shù)微分后的光譜輸入CARS+PLSR構建反演模型，以測試集精度為選擇依據(jù)。為了避免實驗的偶然性影響，確保選出的階數(shù)能使反演精度達到最優(yōu)，每組階數(shù)下均進行了5次反演建模。

圖9顯示了在不同階數(shù)微分條件下，CARS+PLSR反演3種重金屬的測試集R2最大值和平均值。對于Pb，階數(shù)為1.2 時測試集精度的最大值和平均值都達到了峰值（分別為0.7974 和0.7431）；對于Zn，當階數(shù)為0.8 時，R2最大值曲線達到峰值0.8690，此時平均值曲線的值為0.8096，與平均值曲線的峰值0.8126 相差極小，因此將階數(shù)0.8 認為是反演重金屬Zn 的最佳微分階數(shù)；對于Ni，當階數(shù)為0.3 時最大值曲線和平均值曲線均達到峰值（分別為0.8303 和0.7681）。同時考慮精度的最優(yōu)值和平均值，保證選出的最佳階數(shù)能實現(xiàn)反演精度最優(yōu)且模型穩(wěn)定性好。

圖9 不同階數(shù)下反演3種重金屬的測試集精度Fig.9 Estimation accuracy of three heavy metals at different orders

每種重金屬反演時光譜區(qū)間內(nèi)對不同種類重金屬敏感有效的波段不同，即各重金屬的特征波段不同，不同階數(shù)的微分處理后對各區(qū)間波段的特征化效果有所差異，因此各重金屬的最佳微分階數(shù)不盡相同。圖9 中當階數(shù)為0 時，代表不做微分處理使用原始像元光譜，當階數(shù)為1 和2 時相當于使用一階微分和二階微分，3種重金屬的反演結果都顯示出使用最佳階數(shù)分數(shù)階微分的精度高于原始像元光譜、一階微分和二階微分光譜的反演精度，證明了使用分數(shù)階微分能夠有效增強光譜特征，提高重金屬反演精度，并且分數(shù)階微分比整數(shù)階微分能更敏感地突出對土壤重金屬反演有益的光譜信息。根據(jù)本節(jié)的實驗結果，本文在反演重金屬Pb、Zn、Ni的各項實驗中，分別選定1.2、0.8和0.3階的分數(shù)階微分對像元光譜進行處理。

4.2 樣本擴充有效性分析

為驗證樣本擴充的有效性，每種重金屬均進行兩組反演實驗，一組對樣本進行擴充，使用216個樣本參與建模，另一組不擴充樣本，使用72 個樣本，反演每種重金屬時，都以其最佳階數(shù)對像元光譜進行微分處理后，建立CARS+PLSR模型。

表1 中展示了各重金屬的兩組對照實驗結果，在沒有進行樣本擴充前，小樣本問題引起的過擬合現(xiàn)象明顯，3 種重金屬的反演模型訓練集精度R2都大于0.98，然而只有Zn 的測試集精度為0.8178，Pb和Ni的測試集精度都較低，R2小于0.7；對樣本集進行擴充后建立的反演模型，過擬合現(xiàn)象得到了很好的緩解，3種重金屬的訓練集精度和測試集精度的R2差距小于0.05。樣本擴充后測試集的精度都得到了明顯的提升，Pb 的測試集R2從0.6128提升到0.7974，Zn 的精度從0.8178 提升到0.8690，Ni 的R2從0.6969 提升到0.8303。由此可見，樣本擴充不僅緩和了模型的過擬合現(xiàn)象，還有效地提升了3種重金屬的反演精度。

表1 樣本擴充前后的3種重金屬反演精度（ 和 分別代表訓練集和測試集的反演精度R2）Table 1 The estimation accuracies of three heavy metals before and after sample expansion（ and represent the inversion accuracies R2 of the training and test sets，respectively）

圖10 展示了樣本擴充前后重金屬Ni 的反演結果散點圖，樣本中Ni含量大于2000 mg/kg的高含量樣本僅有兩個，從圖10（a）（b）可以看出，兩個高含量樣本均被選入訓練集，測試集中未包含高含量樣本，難以通過測試集精度評估該模型應用于圖像后在高含量區(qū)域的反演能力。從圖10（c）（d）可以看出，樣本擴充后訓練集和測試集都包括了高含量的樣本，因此測試集精度可以代表在高含量區(qū)域的反演效果。在使用3.1 節(jié)所述的樣本擴充條件下，使用KS算法劃分樣本集能夠確保高含量樣本存在于訓練集和測試集，使本研究方法在反演含量偏分布的重金屬時，能夠得到可靠性較高的反演模型。

圖10 樣本擴充前后重金屬Ni反演散點圖Fig.10 Scatter plot of Ni estimation results before and after sample expansion

4.3 波段選擇算法的建模精度對比分析

GA 在PLSR 建模中被認為是一種有效的波段選擇算法（Leardi 和González，1998）。此外，相關系數(shù)法在土壤重金屬反演中被廣泛使用，因此本研究將GA+PLSR、CC+PLSR 兩種建模方法與CARS+PLSR 進行對比分析。參考已有的GA+PLSR研究（Sun 等，2022），GA 的參數(shù)設置為染色體個數(shù)20，迭代次數(shù)150，代際間隙90%，基因變異概率10%。CC 的相關性篩選閾值通過在0—1 范圍內(nèi)，以0.1 為間隔設置閾值進行3 種重金屬的反演實驗，以CC+PLSR 的反演精度最優(yōu)為原則設置閾值的參數(shù)。

通過分數(shù)階微分對擴充后的216條像元光譜進行處理后，分別用CC、GA 和CARS 等3 種算法篩選波段，選出的3 組波段組合分別構建PLSR 反演模型。

表2 展示了3 種重金屬在使用不同波段選擇算法下的反演精度，表中的T（s）代表單次模型構建所花費的秒數(shù)，分別代表訓練集和測試集的反演精度R2。3種建模方法中，CC+PLSR 構建的模型測試集精度最低，甚至在反演Pb 時R2<0.5，模型沒有估算能力，并且在3種重金屬反演中均出現(xiàn)過擬合現(xiàn)象，反演Pb 和Zn 時過擬合最為顯著，訓練集精度遠高于測試集精度。

表2 不同波段選擇算法的重金屬反演精度Table 2 Heavy metal estimation accuracy of different band selection algorithms

GA+PLSR 的模型精度較高，3 種重金屬反演結果僅略低于CARS+PLSR 模型，反演Zn時的測試集精度R2為0.8119，屬于良好模型，反演Pb 和Ni時測試集精度R2也均大于0.78，模型沒有顯示出過擬合現(xiàn)象。但由于GA 在篩選波段時，需要多次迭代選擇最優(yōu)波段組合，每次迭代中所有波段均參與運算，并且由于GA 算法中有種群機制的設置，在本研究的參數(shù)設置下，每次迭代需生成20種波段組合，意味著有20倍的總波段數(shù)參與運算，導致整個模型的構建時間較長，3種重金屬的GA+PLSR反演模型構建時間都在300 s 以上，耗費時間遠遠多于CARS和CC兩種篩選方法。

CARS+PLSR 在3 種建模方法中獲得了最高的反演精度，Zn 和Ni 的反演精度R2都高于0.8，Pb的反演精度R2也達到0.7974，并且模型沒有呈現(xiàn)出明顯的過擬合現(xiàn)象，訓練集和測試集都保持著較高的精度。此外，雖然CARS 和GA 算法都需要多次迭代尋找最優(yōu)波段組合，但由于CARS 內(nèi)的EDF 和ARS 算法在迭代初期就快速、強力的去除重要性低的波段，因此迭代中后期僅有少數(shù)波段參與運算，所需的建模時間顯著縮短。因此CARS是3種方法中最優(yōu)的波段選擇算法，能夠更快篩選出最優(yōu)波段組合，使模型具有更好的反演能力。

4.4 CARS的波段選擇結果分析

圖11 展示了CARS 選出的最優(yōu)波段組合的直方圖，直方圖內(nèi)柱形圖越高代表該波長范圍內(nèi)選擇的波段越多，3 種重金屬的CARS 篩選結果在不同波段區(qū)間選擇的波段數(shù)量不同，說明各波長范圍對不同重金屬反演的貢獻不同。土壤反射光譜分析表明土壤在VNIR-SWIR 區(qū)間的吸收特征主要由土壤有機質、鐵氧化物等土壤光譜活性物質引起（Kooistra 等，2003）。由于土壤光譜活性物質在土壤中對重金屬的吸附具有主導作用，可據(jù)此推算土壤重金屬濃度，這是通過反射光譜間接反演重金屬的主要機理（Rathod 等，2013）。然而有機質、鐵氧化物等土壤光譜活性物質對于不同的重金屬吸附強度不同，因此反演不同重金屬時起主導作用的土壤組分不同（Covelo 等，2007），現(xiàn)有研究表明土壤鐵氧化物對重金屬Pb 的吸附作用是土壤反射光譜反演Pb 含量的主要機理，土壤有機質對重金屬Ni、Zn 的吸附作用是土壤反射光譜反演Ni、Zn 的主要機理，鐵氧化物的吸收特征在500 nm（Wu 等，2007），而600—800 nm 附近的吸收峰被認為是土壤有機質的吸收特征（徐彬彬 等，1991）。已有研究證明500 nm 的鐵氧化物特征波段對反演重金屬Pb 是有效的（張霞 等，2022），Sun 等（Sun 和Zhang，2017；Sun 等，2018）證明了使用600—800 nm的有機質特征譜段對反演Zn和Ni是有效的。根據(jù)圖11所示，在反演Pb時，CARS篩選結果的直方圖中500 nm 范圍內(nèi)的柱形圖最高；在反演Zn 和Ni 時，直方圖內(nèi)最高的柱形圖落在600—800 nm區(qū)間和1600 nm范圍內(nèi)，Li等（2022）指出有機質的官能團對1600 nm 光譜影響顯著，證明1600 nm 同樣是對有機質敏感的特征波段。綜上所述CARS的波段選擇結果與已有反演機理研究保持一致，證明CARS能夠篩選出對反演重金屬有益且合理的特征波段。

圖11 3種重金屬反演模型的CARS優(yōu)選波段結果Fig.11 Results of CARS band selection for three heavy metal estimation models

4.5 重金屬含量分布圖分析

為了保證模型應用到高光譜圖像上時反演結果的可靠性，通過光譜角匹配法（Ramirez-Lopez 等，2013）篩選出圖像上與土壤樣本像元光譜相似度高的區(qū)域，將3種重金屬的最優(yōu)模型分別應用于該區(qū)域進行重金屬含量反演，圖12（b）（c）（d）展示了依據(jù)此方法制作出的3 種重金屬含量分布圖。由于Ni 的含量區(qū)間跨度較大，連續(xù)的圖例難以體現(xiàn)出含量差距，圖12（d）的Ni含量分布圖中采用了分級圖例進行展示，根據(jù)2018 年發(fā)布的《土壤環(huán)境質量 建設用地土壤污染風險管控標準》（GB 36600-2018）相關規(guī)定，礦區(qū)作為二類用地，Ni 含量的篩選值（可能存在風險）為900 mg/kg，因此設置900—2500 mg/kg 的級別代表超過篩選值，用于表示對人體健康可能存在風險的區(qū)域，900 mg/kg 以下采用分位數(shù)法進行分級，確保每一級別內(nèi)所含的像元個數(shù)大致相等。

圖12 重金屬含量分布圖Fig.12 Heavy metal content distribution map

研究區(qū)屬于鎳銅礦山區(qū)域，可以看出Ni 在整個研究區(qū)范圍內(nèi)的含量偏高，但幾乎沒有超過篩選值的情況，大部分區(qū)域的Ni 含量明顯高于Pb 和Zn 的含量值，彭再華和蔣素芳（2018）研究也表明在哈密黃山南礦山中最主要的有價元素為銅和鎳，其他金屬元素含量較低。重金屬Pb 和Zn 具有相似的分布趨勢，李玲等（2020）在新疆礦區(qū)的研究也指出Pb、Zn 含量表現(xiàn)出顯著相關性。3 種重金屬分布的共同點為：在北部以及南部區(qū)域的重金屬含量較高，在靠近居住區(qū)的東部區(qū)域含量較低。

圖13展示了研究區(qū)的高程分布，可以看出研究區(qū)的地勢呈現(xiàn)出東北高，西南低。經(jīng)實地調(diào)研得知圖12（a）中的A 區(qū)域用于堆放采礦渣，圖12（a）中的B區(qū)域是采礦區(qū)，采礦被認為是重金屬污染的最重要因素之一，重金屬會伴隨尾礦渣，廢水等進入土壤（王海洋 等，2022），該采礦區(qū)于2007年開始勘探，已在此地進行了多年的采礦活動，由于研究區(qū)西南地勢較低，重金屬在重力作用下不斷發(fā)生遷移，長時間作用下導致研究區(qū)南部、西南部區(qū)域受到污染，因此3 種重金屬在南部的含量較高。

圖13 研究區(qū)高程圖Fig.13 Elevation map of the study area

研究區(qū)北部含量較高的區(qū)域對應圖13 中高程出現(xiàn)顯著變化的部分，圖12（a）中可看出該區(qū)域圖像與周圍有明顯差異，根據(jù)實地調(diào)研得知該區(qū)域為山體斜坡，坡上的土壤粒徑極小呈細砂狀，坡面平滑。在坡面上采集的4個土壤樣本重金屬含量均不低，其中3 個樣本的Pb、Zn 含量高于平均值，2 個樣本的Ni含量遠高于平均值，且其中1 個達1220 mg/kg，表明該區(qū)域土壤的重金屬含量偏高，反演結果與采樣分析結果一致。

5 結論

目前針對高光譜圖像的重金屬反演研究較少，本文以3 種重金屬Pb、Zn、Ni 為例，開展面向高光譜圖像的反演研究。由于土壤樣本點與圖像土壤像元數(shù)差距懸殊，通過提取采樣點鄰近像元的方法擴充樣本數(shù)據(jù)，同時增加樣本的光譜多樣性；采用FOD 增強對重金屬反演有益的光譜信息，提高反演精度；由于高光譜數(shù)據(jù)的波段眾多，使用CARS 算法選出特征波段以建立反演模型，從降維角度減弱小樣本的影響。以新疆哈密黃山南銅鎳礦區(qū)為研究區(qū)，分析了該反演方法的有效性。

研究結果表明，樣本擴充有效地提升了3種重金屬的測試集反演精度，緩解了模型的過擬合問題，擴充前3種重金屬的訓練集精度R2均在0.98以上，遠高于測試集精度，模型明顯過擬合，而樣本擴充后，3 種重金屬的測試集反演精度都得到了提升，同時，訓練集R2也在0.8 以上，與測試集精度相近，過擬合問題得到顯著改善。

3 種重金屬在最佳階數(shù)分數(shù)階微分下的模型反演精度均高于使用整數(shù)階微分的反演精度，說明分數(shù)階微分能夠更加有效地突出對土壤重金屬反演有益的光譜信息。在最佳階數(shù)微分的基礎上，將CARS與常用的CC和GA兩種波段選擇算法進行了對比，CARS+PLSR 模型的精度最高，并且建模耗時短，本研究認為CARS 是3 種算法中最優(yōu)的波段選擇策略。

對CARS的波段選擇結果進行分析，反演Pb時，選出的波段大多位于500 nm 范圍內(nèi)，該區(qū)間是對Pb 敏感有效的鐵氧化物的特征波段，反演Zn、Ni時，選出的波段大多位于600—800 nm 和1600 nm附近，該波長范圍屬于對Zn、Ni 敏感的有機質特征譜段，與現(xiàn)有反演機理研究相一致，證明CARS算法能夠有效地篩選出光譜中對反演重金屬有益的特征波段。

本研究方法應用于航空高光譜圖像上反演土壤重金屬Pb、Zn、Ni 的含量，反演精度較高，并且反演出的3種重金屬分布圖與實際相符，表明該方法有很好的魯棒性，具有反演多種土壤重金屬的能力。但該方法尚需要在其它礦區(qū)、農(nóng)作區(qū)驗證其適用性，且土壤類型的影響也需要進一步探討。