基于PLS模型的高光譜遙感地質巖性成分反演分析

種紹龍

(甘肅省地質礦產勘查開發局測繪勘查院，甘肅 蘭州 730060)

遙感是通過高空傳感器收集地面的電磁波特性，并根據相關特性進行識別、提取和分析地物的技術，已廣泛應用于金屬礦產勘查與預測、礦化蝕變信息提取、地質構造解釋等工程領域中[1-4]。

高光譜遙感誕生于上世紀80年代，具有分辨率高、圖譜合一、信息量大等一系列優勢，在農業生產、災害預報、環境監測以及地質探查領域獲得了較為成熟的應用，在地質勘查方面，可進行巖性識別、蝕變礦物填圖及信息提取、區域性礦產資源調查等多項工作[5-6]。由于高光譜遙感數據具有波段多、信息量大、相關性強、冗余多、信噪比低等特點，在進行計算分析時，需要對數據進行預處理和特征提取，同時盡可能在保證精度的前提下，降低計算維度，實現圖像的簡易精準化處理[7-8]。在巖性定量識別方面，專家學者提出了主成分分析法、多元線性回歸、典型相關分析等單一或者多種判別法組合應用的先例[9-10]，取得了一定的成果，但是上述方法均有各自的局限之處，如多元線性回歸法無法體現復雜變量之間的相關性、典型相關分析僅能描述兩組變量間的線性相關性、主成分分析在數據量不足時的解釋性較差等缺點，因此，有必要對高光譜遙感的降維計算方法進行改進，以解決上述方法的局限性。

本文擬采用偏最小二乘法(PLS)對高光譜遙感數據進行建模和降維計算，通過對研究區玄武巖和花崗巖的成分反演分析，證明該方法的可行性和準確性。

1 研究區概況

研究區位于甘肅酒泉北山-花牛山-柳園鎮一帶，面積約為500 km2，平均海拔1 800 m，屬溫帶大陸性干旱高原氣候，多年平均氣溫0℃，多年平均降雨量50 mm，受構造作用和巖漿活動影響，區內出露大面積為玄武巖和花崗巖，研究區內成礦地質條件優越，包含金、銅、鉛、鋅、鐵、銀、長石等20余種礦產，因此，開展該區域巖性成分的反演分析，有利于成礦帶的預測和確定，研究區地理及地質狀況如圖1所示。

圖1 研究區地理及地質概況

2 研究技術路線

研究主要分為兩條線路，一是通過機載CASI、SASI、TASI等獲得遙感數據，通過影像預處理和數據降維后，再基于PLS建立相應的分析模型，其中影像預處理主要是通過數據重采樣、大氣校正和幾何校正等幾個步驟，數據重采樣采用雙線性內插法，幾何校正采用二次多項式法，大氣校正采用FLAASH輻射傳輸校正模型；數據降維處理主要包括主成分分析和最小噪聲分離變換法，后者相比前者降低了對噪聲的敏感性，提高圖像質量，因此文中選取最小噪聲分離變換法對數據進行降維處理。二是通過野外巖石樣本采集，獲取ASD光譜測試數據和主量成分測試數據，經光譜重采樣后，與PLS建模分析結果疊加綜合分析后得到研究區的巖性成分反演。如圖2所示。

圖2 研究思路

3 偏最小二乘法回歸模型反演分析法

3.1 方法介紹

偏最小二乘法(Partial Least Square，簡稱PLS)是眾多多元統計方法中的一種，可以對多因變量對多自變量的問題進行回歸分析，是集多元線性回歸、典型相關分析和主成分分析三種方法優缺點為一體的方法，通過建立變量數據組之間的回歸模型，實現遙感數據的降維處理，其基本原理為：

假設經ASD巖石樣本數量為n個，巖石主量元素成分數為m，巖石樣品在不同波段處的數據為p，分別構建標準化變量矩陣X和Y：

X={n1……nm}×m

(1)

Y={n1……nm}×p

(2)

求解矩陣X和Y的第一主成分值t1和u1，要求使兩個矩陣的方差和相關度均達到最大值，即有：

Max{Cov(Xt1,Yu1)}=max

(3)

(4)

式中,w1為XTYYTX對應的最大特征值的特征向量；c1為YTXXTY對應最大特征值的特征向量。

求解得到第一主成分t1和u1后，利用偏最小二乘法分別對X～t1和Y～u1進行回歸分析，假設這一回歸方程的相關性達到理想精度，則可以停止計算，若回歸方程的相關性達不到相關精度要求，則進行第二輪主成分的提取和計算，即重復公式(3)～(4)，直至回歸方程的精度達到要求為止。若最終提取的主成分個數為m個，偏最小二乘回歸模型就會通過實施Yk對所有主成分進行回歸分析，然后再反演表達成Yk關于原變量X1,X2,…，XP的相關回歸方程，其中,k=1,2，…，p。

3.2 反演分析流程

在利用偏最小二乘法進行巖性成分反演分析時，主要包括以下三個主要步驟：(1)采用原始反射率法或者去連續統法等方法對光譜數據進行預處理，消除各種非目標因素對于光譜特征的影響；(2)選擇適宜的主成分個數，既要保證足夠的表達精度和預測方程具有較小的殘差平方和，又要盡量降低計算難度，縮短計算過程；(3)建立預測方程并進行對應評價分析；(4)根據模型計算結果對影像圖進行處理，得到成分含量分布圖。

4 反演結果分析

4.1 預處理方法抉擇

由于光譜數據受地形、光照等的影響較為明顯，因此在進行回歸分析之前需要對其進行預處理，分別使用原始反射率法和去連續統法對研究區玄武巖和花崗巖光譜數據進行建模，并選擇不同的主成分個數進行分析，得到的誤差結果如表1所示。從表中可以看到：誤差值隨主成分個數的增加而降低，在玄武巖光譜數據下，采用原始反射率法預處理后的誤差均大于1，而采用去連續統法處理后，當主成分個數為9個時，其誤差已小于1，故玄武巖的光譜數據建議采用去連續統法進行預處理，同時主成分個數為9個即可；在花崗巖光譜數據下，采用原始反射率法和去連續統法的誤差均小于1，但去連續統法的誤差相對更小，故建議也采用去連續統法對數據進行預處理，且主成分個數僅需2個。綜上所述：去連續統法不僅可以減少計算過程，同時還能提升巖性反演精度。

表1 不同處理方法誤差分析結果

4.2 反演預測值與實測值情況

選定主成分個數之后，分別對玄武巖和花崗巖進行PLS建模分析，其中，玄武巖的ASD樣本個數為26個，花崗巖的ASD樣本個數為30個。通過建模分析分別獲得研究區玄武巖和花崗巖的氧化物預測值和實測值的相關關系，如圖3所示。從圖中可以看到：玄武巖和花崗巖巖性成分的實測值和預測值的擬合度較高(擬合度R>0.8)，其中，玄武巖采用PLS模型反演分析的巖性各氧化物成分含量平均擬合度達到0.91，花崗巖采用PLS模型反演分析的巖性各氧化物成分含量平均擬合度達到0.93，表明本文提出的PLS回歸模型反演分析法具有較高的計算精度，可在實際地質勘探工程中予以合理應用；玄武巖的Si成分含量最高，約為50%左右，Ti成分含量最低，約為1%左右，花崗巖的Si成分含量最高，達到70%，Mg和Ti成分含量最低，僅為0.5%左右。

圖3 反演預測值與實測值關系

4.3 巖性成分分布情況

根據玄武巖和花崗巖各氧化物的模型計算結果，對研究區的影像進行疊加計算，然后去除交叉區域，可獲得不同主量礦物元素的分布情況，如圖4所示。從圖中可以看到：研究區AL元素含量約為15%左右，呈均勻分布、幾乎沒有富集區域；Ca元素主要集中在研究區西南地區，其富集程度達到20%以上，東部地區的Ca元素含量相對較少；Fe元素的分布情況與Ca元素相似，但含量相對較低，約為10%左右，中南局部地區出現Fe元素的富集情況；Mg元素也主要集中于西南地區，其含量達到10%以上；Si元素與其它元素分布不同，其主要分布于研究區中東部地區，含量達到70%以上，西南地區的Si元素含量約為40%左右；Ti元素主要集中于研究區西南地區，但其含量僅為3%左右。

圖4 主量元素分布圖

5 結 論

(1)偏最小二乘法結合了主成分分析法、多元線性回歸和典型相關分析法優勢，能夠在保證計算精度的前提下，有效減少主成分分析個數。相比于原始反射率法，采用去連續統法對光譜數據進行預處理，具有更高的精度和更少的計算過程。

(2)通過對研究區玄武巖和花崗巖的巖性成分反演分析，其主要氧化物預測值與實測值的平均擬合度分別達到0.91和0.93。通過疊加計算，分別得到了六種主流元素的分布情況，可為研究區地層巖性和礦產勘查提供借鑒。

(3)雖然基于PLS模型的反演精度較高，但是受地形和風化因素影響仍然較為明顯，且本文僅對該地區花崗巖和玄武巖兩種巖石進行了巖性反演，存在一定的局限性，將在今后做進一步的補充研究。