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

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

(1.中國地震局地殼應力研究所地殼動力學重點實驗室，北京 100085； 2.中國國土資源航空物探遙感中心，北京 100083; 3.防災科技學院，三河 065201)

0 引言

黃土地層的劃分對于古地震研究具有重要的意義。特別是在我國地震危險性較強、黃土分布較廣的西部黃土高原區以及華北平原和東北的南部，古地震研究都無法避開黃土。由于黃土的粒度與顏色差別小、古生物化石稀少，沉積學參數提取難，采用傳統的目視分層、生物地層劃分和同位素地層學等方法都很難實現黃土地層的劃分對比。長期以來，黃土地層的劃分問題一直是值得深入研究的課題。

在一定的沉積環境下和等時間段內，沉積物地層化學特征具有某些相似性或具有某些特有的、區別于其他沉積環境和其他時間段內沉積的標志。因此，根據沉積物成分變化對地層剖面進行分層是合理的。磁化率是土壤和沉積物的一個重要參數，能夠很好地指示地層韻律變化[1]。但磁化率樣本的采樣位置、數量等在很大程度上制約了其在地層層序識別的深度，其離散的數據也很難說明復雜的地質現象，特別是在進行黃土剖面的地層結構空間展布特征分析時，會出現以點代面、以偏概全的問題。

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

本研究基于高光譜吸收特征參數，探討光譜和磁化率之間的關系，嘗試使用高光譜磁化率模型對地層剖面進行分層。將該模型應用到高光譜影像上，獲得二維空間上連續的磁化率數據，從而為地層結構空間展布分析提供重要依據。

1 數據和方法

本研究以山西口泉一處黃土剖面為例。通過對剖面進行土壤采樣和高光譜影像采集，在實驗室對土壤樣品進行磁化率和光譜測試。基于光譜特征分析，尋找特征譜段，利用多元逐步線性回歸的方法建立特征譜段與磁化率之間的關系模型。將建立的模型應用到影像上，通過對磁化率強度的分類，分析其在地層劃分的有效性。采樣位置和影像獲取范圍示意圖見圖1。

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

Fig.1Photographofthestudiedsectionwiththeprofileslabeled

1.1 數據采集

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

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

1.2 模型建立與驗證

圖2為光譜反射率和磁化率之間的關系。

(a) 不同磁化率反射率變化 (b) 不同磁化率反射率去連續統變化

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

Fig.2Sepctralfeaturesofsamplesandtheirmagneticsusceptibilityvalues

在400～1 000 nm的波譜范圍內，光譜特征略有差異，采用光譜去連續統的方法凸顯光譜差異性，去連續統后，不同磁化率對應光譜在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)為特征光譜參數，發現三者與磁化率線性關系良好(圖3)，特別是x2和x3相關系數均大于0.96。

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

圖3波段比值參數和磁化率的線性關系模型

Fig.3Linearregressionbetweenbandratiosandmagneticsusceptibility

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

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

(1)

將以上關系模型應用到55個C3樣本。研究采用定性和定量的方法分別對關系模型進行驗證，模型精度在很大程度上決定了應用的有效性和準確性。定性評價主要是對磁化率曲線形態進行簡單評價，包括實測磁化率與反演磁化率明顯的轉折點位置的對應關系，圖4為實例磁化率與模型反演磁化率的對比。圖4顯示反演磁化率和實測磁化率整體趨勢上比較一致，峰谷也比較對應。

圖4 實測磁化率和模型反演磁化率對比

定量評價主要采用波譜角分析方法(spectral angle mapper,SAM)和相關系數法。其中SAM法是光譜分析的一種手段，即用光譜匹配程序對預測磁化率與實測磁化率曲線形態進行定量比較，用以評價波譜質量[15]。該算法是將N個樣本點的磁化率看做N維空間向量，通過計算與實測磁化率曲線之間的夾角判定2個磁化率曲線的相似度，夾角越小，說明越相似。相似度用一個得分來表示，得分越接近于1，說明相似度越高[16]。55個測試樣本點的磁化率曲線與對應的實測磁化率曲線相似度的平均得分為0.896，相關系數R2>0.97，均方根誤差RMSE=4.934 47，證明預測磁化率曲線與實測磁化率曲線匹配度較高，該模型的精度較高，能夠很好地應用于磁化率預測。

2 模型應用

2.1 影像質量評價

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

(a) 點1實測光譜與影像光譜對比(b) 點1實測光譜與影像光譜線性關系

(c) 點2實測光譜與影像光譜對比(d) 點2實測光譜與影像光譜線性關系

(e) 點3實測光譜與影像光譜對比(f) 點3實測光譜與影像光譜線性關系

(g) 點4實測光譜與影像光譜對比(h) 點4實測光譜與影像光譜線性關系

圖5UHD185影像反射率與實測反射率對比

Fig.5ReflectancespectrafromtheUHD185imageandfiledreflectancespectra

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

定量評價主要是采用光譜分析手段——SAM法和相關系數法對影像光譜和實測波譜進行定量比較。由前文分析可知，900～950 nm的影像波譜數據不可信，因此在分析時需要去除該譜段數據，去除后影像光譜曲線和實測光譜曲線的相似度SAM>0.9，相關系數R2>0.995。

磁化率模型中使用波段比值作為參數為了進一步對影像數據進行評估，本研究對模型中用到的波段比值b750/b650(x2)和b880/b810(x3)開展了定量評價。精度用相關系數R2和平均比值差來評價，平均比值差計算公式為

(2)

(3)

2.2 模型應用

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

圖6探槽剖面磁化率強度

Fig.6Magneticsusceptibiltyestimateswiththeregressionmodels

為了精確評價基于影像的磁化率填圖效果，用C2的實測磁化率(35個，描述見2.1)與影像反演的磁化率進行對比(圖7)。

(a) C2實測磁化率與反演磁化率對比 (b) 不同深度實測磁化率與反演磁化率對比

圖7影像反演磁化率與實測磁化率對比

Fig.7ComparisonoftheUHD185imageestimatedandmeasuredmagneticsusceptibilty

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

3 結論與討論

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

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

參考文獻(References)：

[1] 衛蕾華.基于高分辨率粒度、磁化率分析的黃土地層劃分與古地震研究[D].北京:中國地震局地質研究所,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] 季峻峰,陳 駿,劉連文,等.洛川黃土中綠泥石的化學風化與磁化率增強[J].自然科學進展,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] 甘甫平,王潤生.遙感巖礦信息提取基礎與技術方法研究[M].北京:地質出版社,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.