基于分?jǐn)?shù)階微分的土壤重金屬高光譜遙感圖像反演

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

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

1 引言

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

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

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

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

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

2.1 研究區(qū)概況

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

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

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

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

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

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

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

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

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

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

3 實(shí)驗(yàn)方法

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

圖5 技術(shù)路線圖Fig.5 Flowchart of the proposed method

3.1 樣本擴(kuò)充

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

3.2 分?jǐn)?shù)階微分（FOD）

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

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

令h=1，式（1）能夠推導(dǎo)出函數(shù)f(x)分?jǐn)?shù)階微分的表達(dá)式為

3.3 競(jìng)爭(zhēng)自適應(yīng)重加權(quán)采樣法（CARS）

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

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

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

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

通過(guò)CARS 算法去掉的變量，其權(quán)重Wi均設(shè)為0；

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

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

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

3.4 模型構(gòu)建與評(píng)價(jià)指標(biāo)

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

精度評(píng)定采用預(yù)測(cè)均方根誤差RMSEP（Root Mean Square Error of Prediction）、相對(duì)分析誤差RPD（Ratio of Prediction to Deviation）和決定系數(shù)R2（Coefficient of Determination）3 個(gè)評(píng)價(jià)指標(biāo)。3 個(gè)評(píng)價(jià)指標(biāo)，RMSEP 值越小，RPD 值越大，R2值越接近1，說(shuō)明反演模型的精度越高；反之，RMSEP值越大，RPD 值越小，R2值越小，說(shuō)明反演模型的精度越低，本文模型優(yōu)劣參考現(xiàn)有的土壤屬性含量高光譜反演的評(píng)價(jià)標(biāo)準(zhǔ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 實(shí)驗(yàn)結(jié)果與分析

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

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

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

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

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

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

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

4.1 分?jǐn)?shù)階微分的最佳階數(shù)確定

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

從圖8中可以看出，隨著階數(shù)的增加，光譜的特征變的越來(lái)越明顯，峰谷差異越發(fā)增大；0—1階的微分結(jié)果顯示出，階數(shù)增大的過(guò)程中，平緩區(qū)間的波段反射率值越發(fā)趨近于0，而波峰波谷區(qū)間的波段反射率值被逐漸放大；1—2 階的微分結(jié)果中，反射率曲線的波形已較為相似（由于1—2階內(nèi)的微分曲線較為相似，圖8中僅展示1階、1.5階和2 階的微分結(jié)果），均與原始光譜曲線差異較大，特征放大的效果明顯，隨著階數(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構(gòu)建反演模型，以測(cè)試集精度為選擇依據(jù)。為了避免實(shí)驗(yàn)的偶然性影響，確保選出的階數(shù)能使反演精度達(dá)到最優(yōu)，每組階數(shù)下均進(jìn)行了5次反演建模。

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

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

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

4.2 樣本擴(kuò)充有效性分析

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

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

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

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

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

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

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

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

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

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

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

4.4 CARS的波段選擇結(jié)果分析

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

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

4.5 重金屬含量分布圖分析

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

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

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

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

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

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

5 結(jié)論

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

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

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

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

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