煤礦沉陷水域重金屬含量高光譜反演

徐良驥，劉曙光，孟雪瑩，韋 任

(1.深部煤礦采動(dòng)響應(yīng)與災(zāi)害防控國(guó)家重點(diǎn)實(shí)驗(yàn)，安徽 淮南 232001; 2.安徽理工大學(xué) 測(cè)繪學(xué)院，安徽 淮南 232001)

近年來(lái)礦區(qū)沉陷水域的研究主要集中在水體的氮、磷含量及分布特征[1-4]、富營(yíng)養(yǎng)化[1,5]、水化學(xué)特征及影響[6-7]、微動(dòng)物群落結(jié)構(gòu)及影響因子[8-9]、水質(zhì)[10]、水環(huán)境評(píng)價(jià)和綜合利用分析[11-12]等;隨著對(duì)重金屬污染的重視，很多學(xué)者研究了海洋[13]、湖泊[14]、礦區(qū)附近河水[15-19]、沉陷水域[20-21]的重金屬含量及分布特征、重金屬污染評(píng)估、有害元素的健康風(fēng)險(xiǎn)評(píng)估[22]、重金屬對(duì)微生物群落及生態(tài)環(huán)境的影響[9]。但關(guān)于沉陷水域中重金屬含量反演模型的建立很少。本工作利用地物光譜儀采集沉陷水域的光譜反射率并對(duì)光譜數(shù)據(jù)進(jìn)行微分變換，采用單波段回歸法、逐步回歸法和波段深度與偏最小二乘結(jié)合3種方法分別建立各重金屬元素含量的預(yù)測(cè)模型，并對(duì)3種回歸模型進(jìn)行精度評(píng)定，確定各重金屬含量的最佳預(yù)測(cè)模型。

1 研究區(qū)概況

潘一礦沉陷水域位于安徽省淮南市潘集區(qū)(地理坐標(biāo)為北緯32.80°～32.83°，東經(jīng)116.78°～116.82°)，該沉陷水域附近存在煤矸石堆、煤矸石充填復(fù)墾的試驗(yàn)田和動(dòng)物養(yǎng)殖，該研究區(qū)包括兩塊沉陷水域，左側(cè)沉陷水域目前用于水產(chǎn)養(yǎng)殖，右側(cè)水域暫無(wú)工農(nóng)業(yè)用途。研究區(qū)水域面積約為0.35 km2，共有43個(gè)采樣點(diǎn)，采樣點(diǎn)均勻分布，各采樣點(diǎn)間距約為150 m(圖1)。

圖1 采樣點(diǎn)分布示意

2 樣品采集與分析

2.1 水體樣品

研究水域內(nèi)采集43個(gè)0～10 cm表層水體樣品，每個(gè)樣品取約2 L。采用原子吸收分光光度法測(cè)定水樣中Cu,Zn,Pb,Cd和Cr五種重金屬元素的含量，采用原子熒光光度法測(cè)定水樣中As元素的含量，測(cè)定結(jié)果見表1。

2.2 光譜數(shù)據(jù)的采集

采用美國(guó)ASD公司生產(chǎn)的地物光譜儀FieldSpec 4在沉陷水域進(jìn)行光譜采集，該儀器的光譜數(shù)據(jù)波段范圍是350～2 500 nm，有兩種采樣間分別為1.4 nm(350～1 000 nm)、2 nm(1 000～2 500 nm);該實(shí)驗(yàn)的重采樣的間隔為1 nm[23]。光譜數(shù)據(jù)為2017-07-27T10:00—12:30 在潘一沉陷水域進(jìn)行采集，每個(gè)采樣點(diǎn)采集20條光譜曲線，并將利用ViewSpecPro軟件剔除掉異常曲線后的反射率的算數(shù)平均值作為樣本的原始反射率光譜值[24]。

2.3 光譜數(shù)據(jù)預(yù)處

水體光譜數(shù)據(jù)采集過(guò)程中受隨機(jī)誤差影響，且重金屬光譜響應(yīng)信號(hào)微弱、原始光譜數(shù)據(jù)難以直接反應(yīng)特征波段的特點(diǎn)，對(duì)原始光譜數(shù)據(jù)進(jìn)行一階微分變換、二階微分變換和倒數(shù)對(duì)數(shù)變換，增強(qiáng)光譜響應(yīng)特征波段，變換后的光譜曲線如圖2所示。微分變換是常見的光譜數(shù)據(jù)處理方法。可在某種程度上消除、削弱光譜數(shù)據(jù)中的噪聲，放大光譜信息，改善多重共線性，提高標(biāo)定方程的性能[25]。計(jì)算公式為

(1)

(2)

表1 水樣的重金屬元素含量統(tǒng)計(jì)特征

Table 1 Statistical characteristics of heavy metal elements in water samples

編號(hào)元素含量/(mg·L-1)Cu2+Pb2+Zn2+AsCd2+Cr6+PB10.00500.00700.05400.00400.00060.0040PB20.02400.01700.11700.06700.00150.0050PB30.01200.01500.15700.02900.00120.0060PB40.00300.00400.14600.01200.00010.0040PB50.00800.01000.06900.04500.00010.0040PB60.00400.02200.01100.00400.00130.0100PB70.00100.01200.00600.03100.00120.0040PB80.00100.00400.05200.01600.00100.0030PB90.02000.01700.06200.03500.00100.0090PB100.00200.00100.04100.02000.00100.0010PB110.00600.00100.10500.00700.00100.0020PB120.00300.00700.03600.04100.00100.0010PB130.00100.02100.01700.00100.00100.0050PB140.00100.00700.00600.02400.00100.0010PB150.00500.00700.05400.00400.00060.0040PB160.02400.01700.11700.06700.00150.0050PB170.01200.01500.15700.02900.00120.0060PB180.00300.00400.14600.01200.00010.0040PB190.00800.01000.06900.04500.00010.0040PB200.00400.02200.01100.00400.00130.0100PB210.00100.01200.00600.03100.00120.0040PB220.00100.00400.05200.01600.00100.0030PB230.00790.01230.07510.03160.00060.0055PB240.00860.01050.08540.01930.00080.0049PB250.00600.00100.10500.00700.00100.0021PB260.00410.00870.09230.01850.00050.0050PB270.00740.01400.05930.01940.00090.0062PB280.01020.01430.07190.03410.00090.0058PB290.02000.01700.06200.03510.00100.0091PB300.00370.01070.05730.02710.00080.0049PB310.00210.01140.03440.02890.00100.0042PB320.00200.00100.04100.02000.00100.0010PB330.00820.00890.06000.02330.00100.0049PB340.00290.00700.03600.04110.00100.0010PB350.00100.02100.01700.00100.00100.0051PB360.00350.00910.08070.02090.00060.0046PB370.00520.00640.06030.01660.00090.0034PB380.00660.00670.08260.01760.00100.0033PB390.00960.01130.08510.02650.00110.0040PB400.00620.01160.04710.03000.00110.0032PB410.00640.00710.05830.02440.00100.0032PB420.00410.01380.05050.01610.00100.0041PB430.00100.00700.00600.02400.00100.0010

圖2 沉陷水域光譜曲線

在對(duì)光譜研究時(shí)，將光譜原始的反射率進(jìn)行倒數(shù)對(duì)數(shù)變換是常用的處理方法，該變換形式可以增強(qiáng)相似光譜間差異，減少因光照條件、環(huán)境變化等因素的影響，計(jì)算公式為

(3)

式中，λi+1,λi,λi-1為相鄰的波長(zhǎng);R(λi+1),R(λi),R(λi-1)為相應(yīng)波長(zhǎng)的原始反射率;R′(λi),R″(λi),Rlgλi分別為波長(zhǎng)λi對(duì)應(yīng)的一階微分、二階微分、倒數(shù)對(duì)數(shù)光譜反射率值。

2.4 統(tǒng)計(jì)分析

采用決定系數(shù)(Coefficient of Determination,R2)、均方根誤差(Root Mean Square Error，RMSE)對(duì)模型的可靠性以及它的預(yù)測(cè)能力進(jìn)行評(píng)價(jià)[26-27]。R2越接近1，表示模型的預(yù)測(cè)值越接近于實(shí)測(cè)值，模擬精度越高。RMSE越接近0，說(shuō)明擬合能力好。

3 結(jié)果與分析

3.1 重金屬含量與水體光譜反射率的相關(guān)性

建立沉陷水域重金屬含量反演模型的首要步驟是統(tǒng)計(jì)分析水體中重金屬元素含量與水體光譜反射率數(shù)據(jù)的相關(guān)性。分析兩類數(shù)值變量的相關(guān)性的常用指標(biāo)為Pearson相關(guān)系數(shù)，計(jì)算公式為

(4)

式中，ri為光譜指標(biāo)反射率與水體中重金屬元素含量的相關(guān)系數(shù);i為波段序號(hào);cov(Xi,Y)為第i個(gè)波段水體光譜指標(biāo)反射率和重金屬含量間協(xié)方差;var(Xi)為第i個(gè)波段光譜指標(biāo)反射率值的方差;var(Y)為水體樣本重金屬含量的方差[26]。

將6種重金屬(Cu,Pb,Zn,As,Cd和Cr)含量分別與原始反射光譜(REF)、一階微分光譜(FDR)、二階微分光譜(SDR)、倒數(shù)對(duì)數(shù)光譜(lg(1/R))進(jìn)行相關(guān)性分析，相關(guān)性分析結(jié)果如圖3所示。由圖3可知:

(1)Pb,As,Cd元素與原始光譜相關(guān)性較低，相關(guān)系數(shù)<0.3;且Cd元素與原始光譜在400～977 nm呈負(fù)相關(guān)。Cr相關(guān)性最大值在1 071 nm處，最大相關(guān)系數(shù)為0.43，為低相關(guān)。Cu,Zn元素含量與原始光譜相關(guān)性達(dá)到中度相關(guān)，分別在波段1 077,1 085 nm處出現(xiàn)最大相關(guān)系數(shù)為0.615 4,0.53。

(2)Cu元素、Cr元素與一階微分光譜的大多數(shù)波段含量呈顯著相關(guān)，相關(guān)系數(shù)波動(dòng)范圍為-0.768～0.874;Zn元素含量與部分波段達(dá)到中度相關(guān)，相關(guān)系數(shù)波動(dòng)范圍為-0.746～0.750;Pb元素、As元素含量與一階微分光譜呈中度相關(guān)，相關(guān)系數(shù)波動(dòng)范圍為-0.444～0.566;Cd元素含量與一階微分光譜呈低度相關(guān)，相關(guān)系數(shù)波動(dòng)范圍為-0.49～0.47。相較于原始光譜的相關(guān)系數(shù)，都有不同程度的提高。

(3)二階微分光譜絕大多數(shù)波段與Cu元素、Zn元素、Cr元素含量達(dá)到高度相關(guān)，相關(guān)系數(shù)在-0.885～0.824;部分波段與Pb元素、As元素、Cd元素含量呈中度相關(guān)，相關(guān)系數(shù)在-0.648～0.764。

(4)倒數(shù)對(duì)數(shù)光譜與各金屬元素相關(guān)系數(shù)較小，與原始光譜的相關(guān)性曲線接近且近似呈鏡面對(duì)稱分布。

3.2 基于單波段的水體重金屬含量估算

選出29個(gè)水樣進(jìn)行建模，14個(gè)水樣進(jìn)行驗(yàn)證。選取400～1 200 nm波段范圍內(nèi)的光譜(REF,FDR,SDR, lg(1/R))與重金屬含量相關(guān)系數(shù)最大的波段建立重金屬含量反演模型，各重金屬的最優(yōu)擬合模型方程(表2)。

圖3 光譜指標(biāo)反射率與重金屬元素含量相關(guān)性

表2 特征單波段與重金屬含量最優(yōu)擬合模型

注:擬合模型方程中，Y為水體重金屬元素含量;x為水體光譜擬合波段的波長(zhǎng)位置的反射率值(下同)。

相比于REF,FDR,lg(1/R)，SDR與各重金屬含量的擬合模型方程效果較好。由表2可知，SDR與Cu,Zn,As和Cr元素含量的二次方程的R2均>0.5;該擬合方程擬合效果較好;SDR與重金屬Cu含量的二次擬合模型精度最高，建模R2為0.796，REMS為0.000 28 μg/L;模型驗(yàn)證R2為0.823，REMS為9.66×10-6mg/L。Pb和Cd元素含量的建模R2均<0.5，驗(yàn)證R2也不高，表明該擬合方程擬合效果不好。

3.3 基于SMLR的水體重金屬含量反演回歸分析

將400～1 200 nm波段內(nèi)的波峰、波谷和(在 0.05 檢驗(yàn)水平上)與水體重金屬含量成顯著相關(guān)的光譜波段選為預(yù)測(cè)模型的備選波段。利用 SPSS 軟件建立擬合模型，得到的擬合模型方程見表3。

由表3可知:同一金屬元素的REF,FDR,SDR,lg(1/R)指標(biāo)下，SDR與重金屬(Cu,Pb,Zn,As和Cr)含量的擬合模型方程的建模R2均>0.7，表明SDR的特征波段建立的擬合方程能較好的擬合光譜數(shù)據(jù)與重金屬含量之間的關(guān)系。Pb,Zn,As和Cd四種元素只有二階微分光譜能與之建立較好的擬合方程。

表3中方差擴(kuò)大因子VIF均<10，表明各多元方程的特征波段之間不存在多重共線性。

表3 水體重金屬含量逐步回歸模型

Table 3 SMLR model for heavy metal content in water

光譜指標(biāo)光譜波段/nm擬合模型方程建模R2RMSE/(mg·L-1)R2驗(yàn)證RMSE/(mg·L-1)共線性統(tǒng)計(jì)容差VIFCuREFFDRSDRlg(1/R)116110771158672719785100511671076Y=-0.070x1161+0.091x1077-0.0030.5058.97×10-60.5431.44×10-5Y=-11.653x1158+0.0010.1201.64×10-50.1372.81×10-5Y=107.217x719-35.718x672-45.799x785-107.217x1005+0.0140.7361.22×10-40.7701.36×10-4Y=0.011x1167-0.035x1076+0.0230.5119.08×10-60.5191.57×10-50.1536.5550.1536.5551.0001.0000.4822.0740.4332.3100.6061.6510.8341.1990.1516.6330.1516.633PbFDRSDR9985658641045Y=7.32x998+0.0050.1732.16×10-50.2103.70×10-5Y=-146.102x565-198.881x864-277.202x1045+0.0010.7456.6×10-60.7721.07×10-51.0001.0000.3313.0240.2264.4250.4502.224ZnREFFDRSDRlg(1/R)10851153687845101311654794836831074Y=0.31x1085+0.0180.2629.46×10-40.2811.72×10-3Y=0.162x1153+0.0560.1910.0013640.3932.42×10-3Y=-5.901x479-501.493x483+116.69x683+433.8486x687-1270.793x845+384.686x1013+636.293x1165+0.1020.7752.89×10-40.7745.37×10-4Y=-0.101x1074+0.1550.2279.91×10-40.2201.86×10-31.0001.0001.0001.0000.2434.1180.3442.9040.4832.0680.5861.7070.2813.5580.2743.6450.2833.5361.0001.000AsFDRSDR418411704Y=27.977x418+0.0270.2861.26×10-40.3202.26×10-4Y=79.1021x411+751.033x704+0.0160.7095.10×10-50.7249.12×10-51.0001.0001.0001.0001.0001.000CdFDRSDR42842298710511138Y=-0.783x428+0.0010.2901.60×10-70.2411.97×10-7Y=0.735x422-2.424x987-8.974x1051-5.251x1138+0.0010.5654.05×10-80.5646.70×10-81.0001.0000.5291.8900.5571.7950.9551.0480.5921.689CrFDRSDR604657566Y=5.433x604+8.025x657+0.0030.7669.91×10-70.7671.68×10-6Y=-28.710x566+0.0030.7811.07×10-60.9021.67×10-60.5961.6770.5961.6771.0001.000

3.4 波段深度與PLSR結(jié)合的水體重金屬含量反演回歸分析

連續(xù)統(tǒng)去除后計(jì)算如下光譜吸收特征:BD(波段深度)、BDR(波段深度比)、NBD(歸一化面積波段深度)、BNA(歸一化面積波段指數(shù))，表達(dá)式為

BD=1-R

(5)

BDR=BD/BDmax

(6)

NBDI=(BD-BDmax)/(BD+BDmax)

(7)

BNA=BD/BDarea

(8)

R(連續(xù)統(tǒng)去除反射率)是光譜反射曲線上的波長(zhǎng)反射率與對(duì)應(yīng)連續(xù)統(tǒng)線上數(shù)值的比值。選取光譜的波段深度作為基礎(chǔ)數(shù)據(jù)，計(jì)算出該光譜的特征指數(shù)，通過(guò)偏最小二乘回歸方法構(gòu)建與水體重金屬含量的預(yù)測(cè)模型，模型精度見表4。

由表4可知:重金屬Cu,Cr元素與BD參數(shù)結(jié)合PLSR建立的模型的決定系數(shù)均大于0.5，表明重金屬Cu,Cr元素含量與該模型預(yù)測(cè)值具有較好的的線性關(guān)系。其中，BNA結(jié)合PLSR的模型的預(yù)測(cè)值與Cu,Cr元素含量的線性關(guān)系最好，決定系數(shù)分別為0.783,0.758，均方根誤差為5×10-7,1.7×10-6mg/L。NBDI,BNA結(jié)合PLSR能建立Zn元素含量反演效果較好的模型;BDR,NBDI,BNA能建立As,Cd元素反演效果較好的模型;NBDI能建立Pb元素含量反演效果較好的模型。

4 結(jié) 論

(1)通過(guò)對(duì)沉陷水域的光譜數(shù)據(jù)進(jìn)行不同形式的變換可知，一階微分變換、二階微分變換、倒數(shù)對(duì)數(shù)變換對(duì)光譜信息都進(jìn)行不同程度的放大。① 從相關(guān)系數(shù)來(lái)看，一階微分變換、二階微分變換效果較為明顯;微分變換后的光譜波段與重金屬含量達(dá)到高相關(guān)度;Cu元素與一階微分光譜的相關(guān)系數(shù)最大，為0.874;Pb元素、Zn元素、As元素、Cd元素和Cr元素與二階微分光譜的相關(guān)系數(shù)最大，相關(guān)系數(shù)分別為-0.648,0.824,0.764,0.636,-0.885。② 根據(jù)單波段回歸和多元逐步回歸建模的決定系數(shù)來(lái)看，二階微分光譜對(duì)6種重金屬(Cu,Pb,Zn,As,Cd和Cr)元素含量的擬合效果效果最好;即二階微分變換的預(yù)處理方法比較適合。

表4 波段深度結(jié)合偏最小二乘反演重金屬結(jié)果

Table 4 BDA combined with PLSR inversion of heavy metal content model accuracy

指標(biāo)BDBDRNBDIBNACuR20.6150.6920.7680.783RMSE/(mg·L-1)4.30×10-61.56×10-51.45×10-45.00×10-7PbR20.1560.3460.5090.465RMSE/(mg·L-1)5.80×10-64.50×10-64.56×10-57.00×10-7ZnR20.4860.2130.5580.643RMSE/(mg·L-1)3.70×10-64.30×10-66.90×10-65.00×10-7AsR20.2330.5590.6780.511RMSE/(mg·L-1)1.12×10-53.64×10-54.50×10-61.40×10-6CdR20.1590.5100.5120.499RMSE/(mg·L-1)1.55×10-54.53×10-51.85×10-57.00×10-7CrR20.5980.6100.7200.758RMSE/(mg·L-1)2.56×10-53.69×10-51.23×10-51.70×10-6

(2)根據(jù)單波段回歸，多元逐步回歸和波段深度與PLSR結(jié)合3種模型的R2和RMSE大小可知，二元微分的單波段回歸對(duì)Cu,Cr元素含量的擬合效果最好;Cu元素建模R2為0.796，RMSE為2.8×10-7mg/L，模型驗(yàn)證的R2為0.823，RMSE為9.66×10-6mg/L;Cr元素建模R2為0.794，RMSE為5.63×10-7mg/L，模型驗(yàn)證的R2為0.806，RMSE為1.65×10-6mg/L;二元微分的多波段逐步回歸對(duì)Pb,Zn,As和Cd四種元素的擬合效果最好，建模R2分別為0.745,0.775,0.709,0.565，RMSE分別為6.6×10-6,2.89×10-4,5.1×10-5,4.05×10-8mg/L;模型驗(yàn)證R2分別為0.772,0.774,0.724,0.564，RMSE分別為1.07×10-5，0.537，9.12×10-5，6.7×10-8mg/L。