基于R语言的光释光年代学数据处理程序包numOSL程序设计及应用实例分析

中国沙漠 (Journal of Desert Research) .2015.

基于R语言的光释光年代学数据处理程序包numOSL程序

设计及应用实例分析

彭俊1，董治宝1，韩凤清2，隆浩3，刘向军2

(1.中国科学院寒区旱区环境与工程研究所 中国科学院沙漠与沙漠化重点实验室, 甘肃 兰州 730000; 2.中国科学院青海盐湖研究所, 青海 西宁 810008；3.中国科学院南京地理

与湖泊研究所，江苏 南京 210008)

摘 要: 光释光测年是测定第四纪沉积物年代的常用方法，年代数据分析与可视化是释光测年过程的重要环节。本文使用开源R程序语言开发了光释光年代学数据分析程序包 numOSL，该程序包是发布于R软件官方网站的标准程序包之一，其设计旨在提供一款可用于释光年代学数据分析的实用软件包。该程序包以R软件为运行平台并以Fortran语言编写核心运算程序，在设计上同时兼顾了程序的流动性、算法稳定性与实用性，其内容涵盖了光释光年代学数据常规数值分析的各个方面，如生长曲线非线性优化、等效剂量计算及蒙托卡罗误差估计、衰减信号曲线分解、快组分等效剂量计算、统计学年龄模型的极大似然估计与马尔科夫链蒙托卡罗采样等。本文对光释光年代学数据分析的常用算法及R程序包numOSL的主要数值分析功能做了简要介绍，并通过具体应用实例展示了该程序包的用法。

关键词：光释光测年; 数据分析; R程序包numOSL; 等效剂量计算; 衰减曲线分解; 快组分等效剂量；年龄模型优化

1 引言

释光测年[1]是测定第四纪沉积物常用的常用测年法。近年来，释光测年在测量仪器 (如Risø释光仪) 与技术方法 (如单片再生剂量法) 上的迅猛发展使其得到普遍应用，被广泛用于各种类型 (如沙漠、河流、湖泊、冰川等) 沉积物年代测定[2-4]。除了仪器与技术的发展，目前通用的年代学数据处理软件 (如软件Viewer、Analyst等) 的开发实现了测量数据的可视化与批量处理，大大提高了数据分析的效率。首先，可靠的数据分析结果是对年代学数据进行合理解释的基础，数值分析方法如非线性参数优化与极大似然估计是光释光数据处理中常用的方法。如指数生长曲线的建立与衰减信号曲线的分解都要用到非线性参数优化技术，统计学年龄模型优化常用到极大似然估计法等，这类算法的主要缺点在于其结果对参数的初始值较为敏感，常常需要使用不同的初始参数反复运算才能确定最佳参数组合，而对于如何确定最佳优化参数，目前的通用软件 (Analyst) 没有做太多考虑，这是本文程序包要解决的问题之一。其次，随机模拟技术近年来也被广泛应用于年代学数据分析，如等效剂量的蒙托卡罗误差评价与年龄模型的马尔科夫链蒙托卡罗采样等都需要进行大量的模拟计算，因此实用高效的数值程序对于提高优化效率与结果的可靠性极为重要，本程序包在设计上通过混合编程极大提高了这些算法的运算效率。此外，软件Analyst 主要用于基本的数据操作处理如数据的导入与导出，等效剂量计算与绘图等，诸如衰减信号曲线分解及统计学年龄模型估计等问题尚需通过其他独立软件 (如SigmaPlot, Matlab等) 实现。本程序包的设计旨在提供一款涵盖释光年代学数据常规数值分析的各个方面的实用程序包，为释光年代学数据分析提供便利。

R语言是一款开源程序语言[5]，被广泛应用于统计模拟、数据挖掘等。R程序语言的第一个显著特点是开源并且完全免费，广泛的用于自然与社会科学的各个分支；R语言的第二个重要特点是支持向量化运算及面向对象程序设计，因此程序的流动性很强；此外，R还提供了与其它程序语言如Fortran、C++ 之间连接的接口，易于实现与其他程序语言间的交互编程。R程序包 numOSL[6] 在设计上兼顾了程序的流动性与算法实用性，旨在提供一款可进行光释光年代数据分析的实用数值程序包。自2013年12月numOSL版本1.1发布于CRAN (The Comprehensive R Archive Network) 以来，笔者对该程序包进行了持续的维护与更新，目前为版本1.7[7]。表1给出了程序包numOSL中的主要函数及用途概述。以下对R程序包numOSL的主要数值计算功能做了简要介绍，并通过实例展示了该程序包的数据优化、分析结果可视化、及数据自动导出功能。

资助项目：国家重大科学研究计划项目 (2013CB956000)；国家973前期计划项目 (2012CB426501)。

作者简介: 彭俊(1987-), 男, 湖北红安人, 博士, 研究方向为风沙地貌学、气候变化与光释光年代学。

E-mail:

表1 R程序包numOSL 版本1.7中的函数

Table 1 A summary of functions in version 1.7 of R package numOSL

函数名称 用途

根据多组衰减信号曲线数据执行常规De计算，功能类似软件Analyst

calED() 给定曲线数据进行De计算，误差估计方法包括简单转换与蒙托卡罗模拟

dbED() De分布的统计描述，包括权重De值、偏度、峰度、百分比分位数

decomp() 基于差分进化与非线性阻尼最小二乘法的衰减信号曲线组分分解

fastED() 根据多组衰减信号曲线数据执行快组分De计算

fitGrowth() 生长曲线优化，包括线性、饱和指数、线性+饱和指数，双饱和指数模型

mcFMM() 有限组分混合年龄模型的马尔科夫链蒙托卡罗采样算法

mcMAM() 最小年龄模型的马尔科夫链蒙托卡罗采样算法

RadialPlotter() 统计学年龄模型的极大似然优化

reportSAM() 统计学年龄模型的马尔科夫链蒙托卡罗采样结果统计汇报

2 等效剂量计算与蒙托卡罗误差模拟

单片再生剂量法 (SAR)[8] 对样品同一单片进行重复预热、人工辐照、光激发等会引起单片感量发生变化，对这种变化若不加以校正，将会给结果带来巨大误差。感量校正是通过紧邻再生剂量之后的恒定实验剂量 (Test Dose) 的释光信号将再生剂量的释光信号进行正规化来实现的。这种标准化的再生信号值随人工再生剂量的变化曲线即为生长曲线。生长曲线通常用指数形式描述：

f(x)?a(1?e?bx)?c (2.1)

式 (2.1) 中a、b、c都是待估参数。a为生长曲线的饱和水平，b为饱和剂量之倒数，c为一个补偿参数。确定了生长曲线的参数之后便可以通过插值或求反解的方法计算自然标准化信号值对应的等效剂量 (De) 值。为了估计生长曲线在误差范围内形状的变化对De值的影响从而更加全面的描述De计算过程中的误差传递，蒙托卡罗 (Monte Carlo) 误差估计法常被用来估计De的误差 [9-10]。该方法需要利用随机数发生器根据测量数据及其误差生成一系列随机自然标准化信号值及人工再生标准化信号值，再利用生成的随机数据重复建立生长曲线并模拟De值计算过程，最后以获得的大批量随机De值的标准偏差作为实际De的误差估计。由于此方法需要重复估计大量随机生成的生长曲线的参数，而生长曲线优化过程中用到的非线性阻尼最小二乘法对初值敏感，优化结果容易呈病态特征，因此对于这种需要进行自动大批量运算的数值优化过程，高效稳定的程序模块显得极为重要。R程序包 numOSL 通过降维处理及随机初始化的方法[6,11]使这种蒙托卡罗误差模拟过程极为稳定，提高了批量数据自动化处理的能力。本程序包中的De计算函数 calED() 不仅能像Analyst 软件那样使用单条生长曲线计算De值 (函数 calED() 与Analyst 软件的单条生长曲线De计算结果非常一致，其对比见文献[6])，还可以根据多个单片的一系列生长曲线建立一条通用生长曲线 (SGC) 并计算一系列De值[12-13]。以下以两组生长曲线数据说明SGC De计算及蒙托卡罗误差模拟函数 calED() 的用法。

图1 函数 calED() 根据两组生长曲线数据计算的8组De值及其误差估计. 图中生长曲线下方的灰色区域

为蒙托卡罗误差传递法模拟的De值的密度分布

Fig. 1 The resulting 8 De values and their standard errors estimated from two growth curves using function calED(). The gray regions under the growth curves denote the kernel density plots of the De values simulated

using the Monte Carlo protocol

library(numOSL) #1

data1<-cbind(c(0, 18, 36, 54, 72, 0, 18), #2

c(0.026, 1.55, 2.39, 3.46, 4.13, 0.023, 1.61), #3

c(0.005, 0.11, 0.27, 0.22, 0.20, 0.008, 0.24)) #4

data2<-cbind(c(0, 18, 36, 54, 72, 0, 18), #5

c(0.021, 1.47, 2.51, 3.37, 4.31, 0.028, 1.52), #6

c(0.003, 0.09, 0.35, 0.21, 0.29, 0.01, 0.17)) #7

data <- rbind(data1, data2) #8

Ltx<-cbind(c(0.5, 1.0, 1.8, 2.3, 2.8, 3.1, 3.6, 4.0), #9

c(0.02, 0.08, 0.11, 0.12, 0.23, 0.31, 0.33, 0.51)) #10

calED(Curvedata=data, Ltx=Ltx, #11

model= "exp", origin=FALSE, #12

nsim=1000, weight=TRUE, outfile="out1") #13 在R平台中输入以上代码则会输出形如图1的图形。第1行代码加载R程序包numOSL； 2-4行代码利用函数analyst() (表1) 或Analyst软件计算的第1个单片建立生长曲线所需3列数据 (第1列为人工再生剂量值、第2列为人工标准化再生信号值、第3列为人工标准化再生信号值的误差估计) 创建一个矩阵；类似的，第5-7行代码利用第2个单片的生长曲线数据创建另一个数据矩阵；第8行代码将两个单片的生长曲线数据合并到一个矩阵中；9-10行代码利用一系列自然标准化OSL信号值生成一个含两列数据的矩阵 (第1列为一系列自

然标准化信号值、第2列为不同自然标准化信号值对应的误差估计)，此处若仅需计算单个自然标准化信号的De值，则第9-10行代码可相应更替为如Ltx<-c(0.5,0.02) 的形式，表示所需估计的De值的自然标准化信号值为0.5，误差为0.02 (数据可通过函数analyst() 或Analyst软件计算获得)；11-13行代码调用函数 calED() 根据给定的自然标准化信号值及生长曲线数据计算一系列De并估计其误差，拟合模型为饱和指数类型 (model="exp")，生长曲线为不过原点形式 (origin=FALSE), 蒙托卡罗误差模拟的最大次数为1,000 (nsim=1000)，生长曲线的优化过程采用权重最小二乘估计 (weight=TRUE), 最后将模拟的8组共8,000个De值保存到当前工作目录下名为 "out1.csv" 的文件中 (outfile="out1")。上例中生长曲线的建立用到了两个不同单片的数据，用户可以通过类似的方式使用多个单片的数据建立一条SGC生长曲线并进行De计算。

3 衰减信号曲线分解及快组分De计算

释光年代样品的光释光信号在持续的光激发下表现为一条信号量随激发时间变化的曲线，称为衰减信号曲线。常见的衰减信号曲线依据激发过程中激发光的强度变化形式不同分为两种：一种为激发过程所使用的激发光强度不随时间变化，衰减曲线称为CW-OSL衰减曲线；另外一种为激发过程中，激发光强度呈线性增加，衰减曲线称为LM-OSL衰减曲线。通常衰减信号曲线可以分解为一系列衰减速率不同的可以由一阶指数形式描述的信号组分, 一般认为这些不同信号组分来自于石英晶格中的不同“陷阱”[14]。依衰减速率不同可将常见的信号组分分为快、中、慢等组分。描述CW-OSL及LM-OSL衰减曲线的数学形式分别为：

ICW(t)?a1b1e?b1t?...?akbke?bkt (3.1)

ILM(t)?a1b1(t/P)e?b1t2/(2P)?...?akbk(t/P)e?bkt2/(2P) (3.2)

式 (3.1) 与 (3.2) 中，I为衰减信号随激发时间t的变化值，a、b分别为某信号组分的陷阱电子数和衰减常数，t为激发时间，P为LM-OSL衰减曲线的最大激发时间; k为总信号组分数。沉积样品的自然衰减曲线中的信号组分数最多可达6~7个。对于组分数为k的指数衰减信号曲线模型，至少有2k个参数需要优化，随着组分数的增加模型的病态越来越严重，对初值及数据结构的要求也越来越高。文献[15] 给出了估计衰减信号曲线组分数k的方法。R程序包numOSL对衰减信号曲线的分解采用了差分进化 (Differential evolution) 与非线性阻尼最小二乘法 (Levenberg Marquardt) 相结合的一种算法[16]，这种算法能有效的对参数值进行自动初始化，进而确定参数的全局最优解，适用于对批量衰减信号曲线数据的自动化分解。参数的误差是通过有限差分接近 (Finite difference approximation) 来估计的。以下代码使用来自文献[17] 的实验室实测 LM-OSL衰减信号曲线数据来说明函数decomp() 在衰减信号曲线分解中的应用，输出的图形见图2。

data(Signaldata) #1

decomp(Sigdata=Signaldata$lm, #2

ncomp=4, constant=TRUE, #3

typ="lm", weight=TRUE, outfile="out2") #4

上述第1行代码加载程序包内置衰减信号曲线数据集， 该LM-OSL衰减信号由两列数据组成 (第1列为激发时间，第2列为激发时间对应的信号序列)； 第2-4行代码调用函数 decomp() 对LM-OSL衰减信号曲线进行分解，曲线被分解为4个组分 (ncomp=4)，分解过

程中自动扣除常量信号组分 (constant=TRUE), 所分解的信号曲线的类型为 LM (typ="lm")，优化过程采用权重优化算法 (weight=TRUE)，最后将分解出的4组信号值输出保存到当前工作目录下名为 "out2.csv" 的文件中 (outfile="out2")。通常衰减信号曲线所含的数据量较大，而此处用于示例的 LM-OSL 衰减曲线已内置于程序包，故不存在数据导入的问题。但用户若想对自己实验获得的衰减信号数据进行组分分解计算，则需自行导入数据。数据导入常用的 R 函数为 read.table(), Windows用户可以使用 Viewer 软件获取衰减信号曲线数据，然后将数据保存到当前工作目录 (一般是 “我的文档”，用户也可自行指定工作目录) 下的 ".txt" 格式文本文件中。例如，若用户希望对获得的CW-OSL衰减曲线进行组分数为3 (不作常量背景信号扣除) 的非权重信号分解，则可以将需要分解的两列数据保存到当前工作目录下名为 "decay" (用户可以自行任意命名) 的 ".txt" 文件中，然后使用read.table() 函数将信号曲线数据读入并执行分解运算，执行代码可相应地更改为: decay <- read.table(file="decay.txt") #1

decomp(Sigdata=decay, #2

ncomp=3, constant=FALSE, #3

typ="cw", weight=FALSE) #4

图2 函数 decomp() 分解的4组分+背景LM-OSL曲线. 右上角为不同信号组分占总信号量的百分比

Fig. 2 The resulting 4 decaying components plus a constant background component for a LM-OSL decay curve decomposed using function decomp(). Proportions of different decaying components are shown in the right

upper region

SAR是目前光释光定年中最常用的一种方法, 该方法计算De值的两个基本前提是：(1) 样品的释光衰减信号的初始部分以快组分信号为主体[17]；(2) 根据实验室人工辐照剂量建立的生长曲线能很好的描述自然条件下样品的信号随辐射剂量的变化[8]。因为自然条件下快组分信号的衰减速率非常快 (一般2s持续光照即可完全晒退)，所以在沉积物埋藏前残余快组分信号很容易被清空，故实验激发过程中产生的快组分信号可以认为完全由埋藏过程中接受的累积辐射剂量所致，快组分信号由此成为释光测年地质计时的理想载体。但通常有两种情况会显著影响到SAR测量结果的可靠性。第一，年轻沉积物样品可能出现不完全晒退现象

[18]，造成沉积前释光信号残余，若残余释光信号占衰减信号曲线的初始部分的比重较大则会造成De值的显著高估[19]。第二，沉积样品的释光衰减信号中可能存在不稳定的信号组分，