基于光谱磁化率模型的黄土剖面地层划分

崔 静， 董新丰， 丁 锐， 张世民， 王琮禾， 鲁恒新， 孙艳云

(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.