基于光譜磁化率模型的黃土剖面地層劃分

崔 靜， 董新豐， 丁 銳， 張世民， 王琮禾， 魯恒新， 孫艷云

(1.中國(guó)地震局地殼應(yīng)力研究所地殼動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室，北京 100085； 2.中國(guó)國(guó)土資源航空物探遙感中心，北京 100083; 3.防災(zāi)科技學(xué)院，三河 065201)

0 引言

黃土地層的劃分對(duì)于古地震研究具有重要的意義。特別是在我國(guó)地震危險(xiǎn)性較強(qiáng)、黃土分布較廣的西部黃土高原區(qū)以及華北平原和東北的南部，古地震研究都無(wú)法避開(kāi)黃土。由于黃土的粒度與顏色差別小、古生物化石稀少，沉積學(xué)參數(shù)提取難，采用傳統(tǒng)的目視分層、生物地層劃分和同位素地層學(xué)等方法都很難實(shí)現(xiàn)黃土地層的劃分對(duì)比。長(zhǎng)期以來(lái)，黃土地層的劃分問(wèn)題一直是值得深入研究的課題。

在一定的沉積環(huán)境下和等時(shí)間段內(nèi)，沉積物地層化學(xué)特征具有某些相似性或具有某些特有的、區(qū)別于其他沉積環(huán)境和其他時(shí)間段內(nèi)沉積的標(biāo)志。因此，根據(jù)沉積物成分變化對(duì)地層剖面進(jìn)行分層是合理的。磁化率是土壤和沉積物的一個(gè)重要參數(shù)，能夠很好地指示地層韻律變化[1]。但磁化率樣本的采樣位置、數(shù)量等在很大程度上制約了其在地層層序識(shí)別的深度，其離散的數(shù)據(jù)也很難說(shuō)明復(fù)雜的地質(zhì)現(xiàn)象，特別是在進(jìn)行黃土剖面的地層結(jié)構(gòu)空間展布特征分析時(shí)，會(huì)出現(xiàn)以點(diǎn)代面、以偏概全的問(wèn)題。

磁化率的高低主要與鐵磁性礦物有關(guān)，磁赤鐵礦、磁鐵礦以及風(fēng)化成壤過(guò)程中的一些含鐵的硅酸鹽礦物(如綠泥石等)都與土壤磁化率強(qiáng)度有關(guān)[2-8]。而鐵氧化物和氫氧化物在可見(jiàn)光/近紅外譜段具有診斷性光譜特征[9]，且被用來(lái)識(shí)別土壤和沉積物中的鐵氧化物，并對(duì)其含量進(jìn)行估算[10-12]。由此推測(cè)磁化率和反射光譜之間可能存在一定關(guān)系。Smith等首次分析了洛川黃土剖面磁化率和光譜特征參量之間的關(guān)系，結(jié)果表明磁化率和光譜的反射率、一階導(dǎo)數(shù)、吸收深度、吸收深度面積具有較高的相關(guān)系數(shù)，并提出具有圖譜合一特性的高光譜影像是未來(lái)開(kāi)展黃土堆積區(qū)地層劃分的一個(gè)研究方向[13]； 但并沒(méi)有利用關(guān)系模型進(jìn)行磁化率反演或應(yīng)用到高光譜影像上驗(yàn)證其可行性。

本研究基于高光譜吸收特征參數(shù)，探討光譜和磁化率之間的關(guān)系，嘗試使用高光譜磁化率模型對(duì)地層剖面進(jìn)行分層。將該模型應(yīng)用到高光譜影像上，獲得二維空間上連續(xù)的磁化率數(shù)據(jù)，從而為地層結(jié)構(gòu)空間展布分析提供重要依據(jù)。

1 數(shù)據(jù)和方法

本研究以山西口泉一處黃土剖面為例。通過(guò)對(duì)剖面進(jìn)行土壤采樣和高光譜影像采集，在實(shí)驗(yàn)室對(duì)土壤樣品進(jìn)行磁化率和光譜測(cè)試?；诠庾V特征分析，尋找特征譜段，利用多元逐步線性回歸的方法建立特征譜段與磁化率之間的關(guān)系模型。將建立的模型應(yīng)用到影像上，通過(guò)對(duì)磁化率強(qiáng)度的分類(lèi)，分析其在地層劃分的有效性。采樣位置和影像獲取范圍示意圖見(jiàn)圖1。

圖1探槽剖面與采樣點(diǎn)位置(綠色點(diǎn)為影像光譜驗(yàn)證點(diǎn))

Fig.1Photographofthestudiedsectionwiththeprofileslabeled

1.1 數(shù)據(jù)采集

該剖面位于山西省懷仁縣西北部一處黃土臺(tái)地上，剖面中心位置經(jīng)緯度坐標(biāo)為113° 1′13.10″E，39° 51′33.04″N。該剖面自上而下分別為耕植土(L0)、全新世灰黑色黑壚土(S0)、上更新統(tǒng)馬蘭黃土(L1)和礫石層。從古土壤頂層開(kāi)始，設(shè)計(jì)3條測(cè)線自上而下每隔10 cm進(jìn)行樣品采集，3條測(cè)線分別命名為C1,C2和C3。其中C1和C2的樣品用來(lái)建模，樣品數(shù)分別為52個(gè)和35個(gè)，總數(shù)為87個(gè)； C3的樣品用來(lái)做模型驗(yàn)證，樣品總數(shù)為55個(gè)。低頻磁化率采用美國(guó)AGICPOI公司生產(chǎn)的Kappabridge MFK1-FA各向磁化率儀進(jìn)行樣品測(cè)試，選取的頻率為976 Hz，得到的磁化率單位為10-11m3·kg-1。反射率光譜采用美國(guó)ASD公司生產(chǎn)的FieldSpec光譜儀對(duì)樣品進(jìn)行測(cè)試。野外波譜測(cè)試時(shí)，要求每個(gè)測(cè)試點(diǎn)地物盡量均一，面積大于等于17.8 cm2，每個(gè)測(cè)點(diǎn)都選取5個(gè)位置進(jìn)行測(cè)量。為保證儀器的穩(wěn)定性和盡可能多地去除儀器噪聲，每次測(cè)量記錄30條連續(xù)波譜，求取30條波譜的均值作為該樣品的測(cè)量值。在剖面上選取4驗(yàn)證個(gè)點(diǎn)進(jìn)行高光譜影像的光譜驗(yàn)證，4個(gè)驗(yàn)證點(diǎn)分別從暗色地物向亮色地物過(guò)渡(圖1)。

2015年5月20日，天氣晴朗無(wú)風(fēng)，北京時(shí)間11：00—14：00，采用德國(guó)Cubert公司生產(chǎn)的UHD185機(jī)載光譜成像儀進(jìn)行光譜影像采集。視場(chǎng)角范圍為27°，光譜范圍為450～950 nm，光譜分辨率為4 nm。傳感器垂直探槽剖面，距離為9 m，空間分辨率為4.22 cm。

1.2 模型建立與驗(yàn)證

圖2為光譜反射率和磁化率之間的關(guān)系。

(a) 不同磁化率反射率變化 (b) 不同磁化率反射率去連續(xù)統(tǒng)變化

圖2不同磁化率光譜反射率變化

Fig.2Sepctralfeaturesofsamplesandtheirmagneticsusceptibilityvalues

在400～1 000 nm的波譜范圍內(nèi)，光譜特征略有差異，采用光譜去連續(xù)統(tǒng)的方法凸顯光譜差異性，去連續(xù)統(tǒng)后，不同磁化率對(duì)應(yīng)光譜在650～750 nm和810～880 nm之間斜率明顯不同，隨著磁化率的增加，b750/b650和b880/b810逐漸增加，500～600 nm之間光譜也有略微的差異特征，但不是很明顯。這些特征是由于土壤中鐵氧化物的含量差異造成的[9,11,14]。

選擇b600/b500(x1)、b750/b650(x2)、 b880/b810(x3)為特征光譜參數(shù)，發(fā)現(xiàn)三者與磁化率線性關(guān)系良好(圖3)，特別是x2和x3相關(guān)系數(shù)均大于0.96。

(a) b600/b500與磁化率關(guān)系 (b) b750/b650與磁化率關(guān)系 (c) b880/b810與磁化率關(guān)系

圖3波段比值參數(shù)和磁化率的線性關(guān)系模型

Fig.3Linearregressionbetweenbandratiosandmagneticsusceptibility

因此將這2個(gè)特征參數(shù)作為自變量，采用多元逐步線性回歸的方法建立模型，模型表達(dá)式為

y=156.031x2+611.195x3-720.957 78，R2=0.984,RMSE=3.875 09。

(1)

將以上關(guān)系模型應(yīng)用到55個(gè)C3樣本。研究采用定性和定量的方法分別對(duì)關(guān)系模型進(jìn)行驗(yàn)證，模型精度在很大程度上決定了應(yīng)用的有效性和準(zhǔn)確性。定性評(píng)價(jià)主要是對(duì)磁化率曲線形態(tài)進(jìn)行簡(jiǎn)單評(píng)價(jià)，包括實(shí)測(cè)磁化率與反演磁化率明顯的轉(zhuǎn)折點(diǎn)位置的對(duì)應(yīng)關(guān)系，圖4為實(shí)例磁化率與模型反演磁化率的對(duì)比。圖4顯示反演磁化率和實(shí)測(cè)磁化率整體趨勢(shì)上比較一致，峰谷也比較對(duì)應(yīng)。

圖4 實(shí)測(cè)磁化率和模型反演磁化率對(duì)比

定量評(píng)價(jià)主要采用波譜角分析方法(spectral angle mapper,SAM)和相關(guān)系數(shù)法。其中SAM法是光譜分析的一種手段，即用光譜匹配程序?qū)︻A(yù)測(cè)磁化率與實(shí)測(cè)磁化率曲線形態(tài)進(jìn)行定量比較，用以評(píng)價(jià)波譜質(zhì)量[15]。該算法是將N個(gè)樣本點(diǎn)的磁化率看做N維空間向量，通過(guò)計(jì)算與實(shí)測(cè)磁化率曲線之間的夾角判定2個(gè)磁化率曲線的相似度，夾角越小，說(shuō)明越相似。相似度用一個(gè)得分來(lái)表示，得分越接近于1，說(shuō)明相似度越高[16]。55個(gè)測(cè)試樣本點(diǎn)的磁化率曲線與對(duì)應(yīng)的實(shí)測(cè)磁化率曲線相似度的平均得分為0.896，相關(guān)系數(shù)R2>0.97，均方根誤差RMSE=4.934 47，證明預(yù)測(cè)磁化率曲線與實(shí)測(cè)磁化率曲線匹配度較高，該模型的精度較高，能夠很好地應(yīng)用于磁化率預(yù)測(cè)。

2 模型應(yīng)用

2.1 影像質(zhì)量評(píng)價(jià)

影像數(shù)據(jù)的準(zhǔn)確性直接影響著應(yīng)用的效果，所以首先需要對(duì)數(shù)據(jù)質(zhì)量進(jìn)行評(píng)價(jià)。在剖面上選取4個(gè)點(diǎn)進(jìn)行高光譜影像的光譜驗(yàn)證，驗(yàn)證點(diǎn)位置見(jiàn)圖1(綠色點(diǎn))。本次研究同樣分別采用定性和定量的方法對(duì)UHD185反射率數(shù)據(jù)進(jìn)行質(zhì)量評(píng)價(jià)，其中定性評(píng)價(jià)主要是從反射率譜形上對(duì)其進(jìn)行簡(jiǎn)單評(píng)價(jià)，包括影像波譜與實(shí)測(cè)波譜明顯吸收位置的對(duì)應(yīng)關(guān)系(圖5)。

(a) 點(diǎn)1實(shí)測(cè)光譜與影像光譜對(duì)比(b) 點(diǎn)1實(shí)測(cè)光譜與影像光譜線性關(guān)系

(c) 點(diǎn)2實(shí)測(cè)光譜與影像光譜對(duì)比(d) 點(diǎn)2實(shí)測(cè)光譜與影像光譜線性關(guān)系

(e) 點(diǎn)3實(shí)測(cè)光譜與影像光譜對(duì)比(f) 點(diǎn)3實(shí)測(cè)光譜與影像光譜線性關(guān)系

(g) 點(diǎn)4實(shí)測(cè)光譜與影像光譜對(duì)比(h) 點(diǎn)4實(shí)測(cè)光譜與影像光譜線性關(guān)系

圖5UHD185影像反射率與實(shí)測(cè)反射率對(duì)比

Fig.5ReflectancespectrafromtheUHD185imageandfiledreflectancespectra

從圖5中可以明顯看出ASD實(shí)測(cè)波譜與UHD185影像波譜整體趨勢(shì)上比較一致，吸收位置也比較對(duì)應(yīng)。但是900～950 nm波段范圍的反射率整體是下降的，與實(shí)際不符，說(shuō)明這些通道的數(shù)據(jù)不可信。

定量評(píng)價(jià)主要是采用光譜分析手段——SAM法和相關(guān)系數(shù)法對(duì)影像光譜和實(shí)測(cè)波譜進(jìn)行定量比較。由前文分析可知，900～950 nm的影像波譜數(shù)據(jù)不可信，因此在分析時(shí)需要去除該譜段數(shù)據(jù)，去除后影像光譜曲線和實(shí)測(cè)光譜曲線的相似度SAM>0.9，相關(guān)系數(shù)R2>0.995。

磁化率模型中使用波段比值作為參數(shù)為了進(jìn)一步對(duì)影像數(shù)據(jù)進(jìn)行評(píng)估，本研究對(duì)模型中用到的波段比值b750/b650(x2)和b880/b810(x3)開(kāi)展了定量評(píng)價(jià)。精度用相關(guān)系數(shù)R2和平均比值差來(lái)評(píng)價(jià)，平均比值差計(jì)算公式為

(2)

(3)

2.2 模型應(yīng)用

將模型應(yīng)用到UHD185影像上，得到了磁化率強(qiáng)度(圖6)。為了便于肉眼識(shí)別，本研究使用ArcMap自帶的標(biāo)準(zhǔn)差法選擇1倍標(biāo)準(zhǔn)差間隔進(jìn)行分級(jí)顯示。磁化率強(qiáng)度可以明顯地將地層分為6層，自上而下表現(xiàn)為： 藍(lán)綠混合層、橙色層、黃色層、綠色層、藍(lán)色層和綠色層。將肉眼劃分的剖面層序分界點(diǎn)(圖6粉色虛線和黑色虛線)與磁化率強(qiáng)度分級(jí)點(diǎn)(圖6藍(lán)色實(shí)線)顯示的層序?qū)Ρ?，可以看出磁化率?qiáng)度分級(jí)將S0和L1之間劃分為2層，而且這2層的分界點(diǎn)對(duì)應(yīng)于磁化率曲線的轉(zhuǎn)折點(diǎn)。這表明磁化率強(qiáng)度在縱向上的波動(dòng)特征與地層層序的旋回性有較好的對(duì)應(yīng)關(guān)系。相較于單測(cè)線測(cè)試磁化率進(jìn)行地層劃分，二維磁化率強(qiáng)度可以直觀地展現(xiàn)地層結(jié)構(gòu)空間展布特征，避免局部由于生物效應(yīng)等引起的磁化率變化異常，造成誤判。例如圖6中A和B部分磁化率曲線有所不同，若利用其對(duì)應(yīng)的磁化率劃分地層，A和B的結(jié)果將會(huì)不同。但是從磁化率強(qiáng)度上，可以宏觀地識(shí)別出A和B為同一地層。需要注意的是，本次研究中C處有一個(gè)樹(shù)根，磁化率強(qiáng)度呈現(xiàn)出局部高值，并不能代表整個(gè)地層的特征。

圖6探槽剖面磁化率強(qiáng)度

Fig.6Magneticsusceptibiltyestimateswiththeregressionmodels

為了精確評(píng)價(jià)基于影像的磁化率填圖效果，用C2的實(shí)測(cè)磁化率(35個(gè)，描述見(jiàn)2.1)與影像反演的磁化率進(jìn)行對(duì)比(圖7)。

(a) C2實(shí)測(cè)磁化率與反演磁化率對(duì)比 (b) 不同深度實(shí)測(cè)磁化率與反演磁化率對(duì)比

圖7影像反演磁化率與實(shí)測(cè)磁化率對(duì)比

Fig.7ComparisonoftheUHD185imageestimatedandmeasuredmagneticsusceptibilty

圖7顯示影像預(yù)測(cè)的磁化率整體曲線特征與實(shí)測(cè)一致，且模型反演的影像磁化率與實(shí)測(cè)磁化率具有良好的線性關(guān)系，相關(guān)系數(shù)均大于0.95。但是影像反演的磁化率和實(shí)測(cè)磁化率存在一定的偏移，模型反演的值表現(xiàn)為高估實(shí)際的磁化率。一種可能是對(duì)目標(biāo)物觀測(cè)的方式不同，造成了這些差異，但是仍然需要對(duì)其原因做進(jìn)一步的研究。整體上本文的磁化率光譜模型表現(xiàn)出具有相對(duì)較高的精度。

3 結(jié)論與討論

本文基于高光譜磁化率模型對(duì)黃土剖面地層劃分進(jìn)行了探索,建立了光譜與指示地層韻律變化的磁化率之間的關(guān)系模型。該模型可以有效地反演磁化率，具有較高的精度： 實(shí)驗(yàn)室測(cè)試樣品估計(jì)磁化率與實(shí)測(cè)磁化率的相關(guān)系數(shù)R2>0.97，剖面高光譜影像估計(jì)磁化率和實(shí)測(cè)值相關(guān)系數(shù)R2>0.95。反演得到的磁化率強(qiáng)度分布不僅可以識(shí)別出肉眼識(shí)別的地層，還能將肉眼無(wú)法識(shí)別的黃土與古土壤的過(guò)渡層識(shí)別出來(lái)，其在縱向上的波動(dòng)特征與地層層序的旋回性有較好的對(duì)應(yīng)關(guān)系。因此，基于光譜磁化率模型的高光譜影像有利于實(shí)現(xiàn)黃土地層的精細(xì)劃分。

通常情況下，古地震事件分析對(duì)應(yīng)的剖面需要較高的空間分辨率，剖面影像需要多幅影像拼接而成。然而影像獲取的條件很難保證完全一致，高光譜影像多個(gè)波段的鑲嵌融合是首先要解決的問(wèn)題。利用高光譜影像反演得到的磁化率是單波段影像，鑲嵌融合技術(shù)要求較低，分類(lèi)方法也較成熟。因此，本研究并沒(méi)有直接使用反射波譜進(jìn)行地層的直接劃分。同時(shí)不同時(shí)間尺度的地層層序尺度不同，即還有分層閾值的設(shè)定問(wèn)題，本研究并沒(méi)有將地層層序和時(shí)間關(guān)聯(lián)起來(lái)，只是探索了利用高光譜進(jìn)行地層劃分的可行性。與時(shí)間關(guān)聯(lián)的層序劃分對(duì)于古地震定量分析具有重要作用，因此后續(xù)研究可以根據(jù)不同的應(yīng)用需求，開(kāi)展層序序列分析。

參考文獻(xiàn)(References)：

[1] 衛(wèi)蕾華.基于高分辨率粒度、磁化率分析的黃土地層劃分與古地震研究[D].北京:中國(guó)地震局地質(zhì)研究所,2015.

Wei L H.Loess Stratigraphy and Paleo-Earthquake Identification Based on High-Resolution Analysis of Granularity and Magnetic Susceptibility[D].Beijing:Institute of Geology,China Earthquake Administration,2015.

[2] Balsam W,Ji J F,Chen J.Climatic interpretation of the Luochuan and Lingtai loess sections,China,based on changing iron oxide mineralogy and magnetic susceptibility[J].Earth and Planetary Science Letters,2004,223(3/4):335-348.

[3] Chen J,Ji J F,Balsam W,et al.Characterization of the Chinese loess-paleosol stratigraphy by whiteness measurement[J].Palaeogeography,Palaeoclimatology,Palaeoecology,2002,183(3/4):287-297.

[4] Liu X M,Hesse P,Rolph T.Origin of maghaemite in Chinese loess deposits:Aeolian or pedogenic?[J].Physics of the Earth and Planetary Interiors,1999,112(3/4):191-201.

[5] Maher B A,Thompson R.Mineral magnetic record of the Chinese loess and paleosols[J].Geology,1991,19(1):3-6.

[6] Maher B A.Magnetic properties of modern soils and Quaternary loessic paleosols:Paleoclimatic implications[J].Palaeogeography,Palaeoclimatology,Palaeoecology,1998,137(1/2):25-54.

[7] Zhou L P,Oldfield F,Wintle A G,et al.Partly pedogenic origin of magnetic variations in Chinese loess[J].Nature,1990,346(6286):737-739.

[8] 季峻峰,陳 駿,劉連文,等.洛川黃土中綠泥石的化學(xué)風(fēng)化與磁化率增強(qiáng)[J].自然科學(xué)進(jìn)展,1999,9(7):619-623.

Ji J F,Chen J,Liu L W,et al.Chemical weathering of chlorite and enhancement of magnetic susceptibility in Luochuan Loess[J].Progress in Natural Science,1999,9(7):619-623.

[9] 甘甫平,王潤(rùn)生.遙感巖礦信息提取基礎(chǔ)與技術(shù)方法研究[M].北京:地質(zhì)出版社,2004:43-47.

Gan F P,Wang R S.Basis Theory and Technical Methods Study of Remote Sensing Rock and Mineral Information Extraction[M].Beijing:Geological Publishing House,2004:43-47.

[10] Deaton B C,Balsam W L.Visible spectroscopy：A rapid method for determining hematite and goethite concentration in geological materials[J].Journal of Sedimentary Petorology,1991,61(4):628-632.

[11] Ji J F,Balsam W,Chen J,et al.Rapid and quantitative measurement of hematite and goethite in the Chinese loess-epaleosol sequence by diffuse reflectance spectroscopy[J].Clays and Clay Minerals,2002,50(2):208-216.

[12] Scheinost A C,Chavernas A,Barrón V,et al.Use and limitation of second-derivative diffuse reflectance spectroscopy in the visible to near-infrared range to identify and quantify Fe oxide minerals in soils[J].Clays and Clay Minerals,1998,46(5):528-536.

[13] Smith M J,Stevens T,MacArthur A,et al.Characterising Chinese loess stratigraphy and past monsoon variation using field spectroscopy[J].Quaternary International,2011,234(1/2):146-158

[14] Hunt G R,Salisbury J W,Lenhoff C J.Visible and near-infrared spectra of minerals and rocks:III.Oxides and hydroxides[J].Modern Geology,1971,2:195-205.

[15] Cui J,Yan B K,Wang R S,et al.Regional-scale mineral mapping using ASTER VNIR/SWIR data and validation of reflectance and mineral map products using airborne hyperspectral CASI/SASI data[J].International Journal of Applied Earth Observation and Geoinformation,2014,33:127-141.

[16] Kruse F A,Lefkoff A B,Boardman J W,et al.The spectral image processing system(SIPS)-interactive visualization and analysis of imaging spectrometer data[J].Remote Sensing of Environment,1993,44(2/3):145-163.