丁松滔,張霞,尚坤,李儒,孫偉超
1.中國(guó)科學(xué)院空天信息創(chuàng)新研究院,北京 100101;2.中國(guó)科學(xué)院大學(xué),北京 100049;3.自然資源部國(guó)土衛(wèi)星遙感應(yīng)用中心,北京 100048
成像高光譜遙感能夠獲取圖像上連續(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)最高的波段范圍及其物理意義,分析了重金屬含量分布圖的可靠性。
研究區(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
航空高光譜圖像獲取于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
圖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
由于研究區(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ù)。
分?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á)式為
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)算。
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。
反演中將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
由于分?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)行處理。
為驗(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
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)波段組合,使模型具有更好的反演能力。
圖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
為了保證模型應(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é)果一致。
目前針對(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)一步探討。