1 案例（一）

定量地层学-Graphic correlation图形相关：在使用图解的、半定量的方法对成对井进行人工对比的基础上，对许多井或剖面中事件进行地层对比。

1.1 构造数据

Oslo <- c(0, 1, 0.9, 3.9, 4, 6, 0.3, 4)
SE_Scania <- c(0.7, 6.2, 5, 5.2, 10.5, 12.5, 4.5, 13.1)
Bornholm <- c(1, 2, 1, 1.1, 2.5, 3.1, 0, 4.1)
Speci <- c('M.simon F', 'M.simon L', 'D.acut F', 'D.acut L', 'I.sarsi F', 'I.sarsi L', 'Cyrt F', 'Cyrt L')
data <- data.frame(Speci = Speci, Oslo = Oslo, SE_Scania = SE_Scania, Bornholm = Bornholm) #生成数据框
head(data)

##       Speci Oslo SE_Scania Bornholm
## 1 M.simon F  0.0       0.7      1.0
## 2 M.simon L  1.0       6.2      2.0
## 3  D.acut F  0.9       5.0      1.0
## 4  D.acut L  3.9       5.2      1.1
## 5 I.sarsi F  4.0      10.5      2.5
## 6 I.sarsi L  6.0      12.5      3.1

数据为Nielsen(1995)的奥陶纪三叶虫数据。他在三个地方记录的四种物种的首次和最后一次出现:Bornholm(丹麦);东南斯堪尼亚(瑞典南部);和Slemmestad(挪威奥斯陆地区)。向量Speici表示的是四个不同的物种，F和L表示的是该物种在剖面中延续的深度最小和最大值。

1.2 画图

par(mfrow = c(2, 2), mar = c(4, 4, 1, 1))
plot(data$Oslo, data$SE_Scania, xlab = 'Oslo', ylab = 'SE Scania', ylim = c(0, 14), xlim = c(-3, 7), pch = 19, bty = 'l')
abline(lm(data$SE_Scania ~ data$Oslo), col = 'black', lwd = 2) #添加趋势线
text(data$Oslo + 0.5, data$SE_Scania - 0.5, data$Speci, cex = 0.8) #添加文字
text(-2.75, 13.8, '(a)', cex = 1.5, font = 2) #左上角添加'(a)'
plot(data$Oslo, data$SE_Scania, xlab = 'Oslo', ylab = 'SE Scania', ylim = c(0, 14), xlim = c(-3, 7), pch = 19, bty = 'l')
abline(lm(data$SE_Scania ~ data$Oslo), col = 'black', lwd = 2) #添加趋势线
text(data$Oslo + 0.5, data$SE_Scania - 0.5, data$Speci, cex = 0.8) #添加文字
arrows(x0 = 0, y0 = 0.7, x1 = -1.4, y1 = 0.7, lwd = 2, angle = 20) ; arrows(x0 = 3.9, y0 = 5.2, x1 = 3.9, y1 = 9.4, lwd = 2, angle = 20) ; arrows(x0 = 4, y0 = 13.1, x1 = 6, y1 = 13.1, lwd = 2, angle = 20) #添加箭头
text(-2.75, 13.8, '(b)', cex = 1.5, font = 2) #左上角添加'(b)'
#计算Composite Oslo-Scania
lm(data$SE_Scania ~ data$Oslo) # y = 3.092 + 1.64 * x (其中y为变量SE Scania,x为Oslo)

## 
## Call:
## lm(formula = data$SE_Scania ~ data$Oslo)
## 
## Coefficients:
## (Intercept)    data$Oslo  
##       3.092        1.640

Composite_Scania <- (data$SE_Scania - 3.092) / 1.64 
Composite_Oslo_Scania <- c(Composite_Scania[1 : 3], Oslo[4 : 5], Composite_Scania[6 : 8]) #在图b中箭头为'↑'或'↓'值为Oslo(原始值),箭头为'←'或'→'值为Composite_Scania
plot(Composite_Oslo_Scania, data$Bornholm, xlab = 'Composite Oslo-Scania', ylab = 'Bornholm', ylim = c(0, 4.5), xlim = c(-3, 8), pch = 19, bty = 'l')
abline(lm(data$Bornholm ~ Composite_Oslo_Scania), col = 'black', lwd = 2) #添加趋势线
text(Composite_Oslo_Scania + 0.5, data$Bornholm + 0.2, data$Speci, cex = 0.8) #添加文字
text(-2.75, 4.3, '(c)', cex = 1.5, font = 2) #左上角添加'(c)'
plot(Composite_Oslo_Scania, data$Bornholm, xlab = 'Composite Oslo-Scania', ylab = 'Bornholm', ylim = c(0, 4.5), xlim = c(-3, 8), pch = 19, bty = 'l')
abline(lm(data$Bornholm ~ Composite_Oslo_Scania), col = 'black', lwd = 2) #添加趋势线
text(Composite_Oslo_Scania + 0.5, data$Bornholm + 0.2, data$Speci, cex = 0.8) #添加文字
arrows(x0 = -1.46, y0 = 1, x1 = -1.46, y1 = 0.15, lwd = 2, angle = 20) ; arrows(x0 = 0.86, y0 = 0, x1 = -1.65, y1 = 0, lwd = 2, angle = 20) ; arrows(x0 = 3.9, y0 = 1.1, x1 = 3.9, y1 = 2.3, lwd = 2, angle = 20) ; arrows(x0 = 1.9, y0 = 2, x1 = 3.1, y1 = 2, lwd = 2, angle = 20) ; arrows(x0 = 6.1, y0 = 4.1, x1 = 8.2, y1 = 4.1, lwd = 2, angle = 20)#添加箭头
text(-2.75, 4.3, '(d)', cex = 1.5, font = 2) #左上角添加'(d)'

图形相关过程首先选择一个高质量的部分作为生成复合标准参考的起点。我们选择了Oslo部分。图(a)为奥斯陆和斯堪尼亚地区首次和最后一次出现事件的地层位置，单位为米，我们选择使用一条直线相关(LOC)。然后将事件转移到LOC上图(b)，始终通过扩展范围，从不缩短范围，事件在LOC上的位置构成了奥斯陆-斯堪尼亚组合的序列。接下来，以相同的方式将Bornholm局部集成到复合材料中图(c)和图(d)。理想情况下，我们应该多次重复这个过程，从添加Oslo局部开始，直到序列稳定。

2 案例（二）

定量地层学-Constrained optimization约束最优：在许多井或剖面中地层事件的对比，以及这些事件在时间尺度或沉积物厚度上的可选定位。

图形相关法有许多基本问题和局限性。它基本上是一种手工方法，是一项耗时的工作，充满了关于相关线构建的主观选择。此外，位置一次一个地集成到解决方案(组合部分)中，它们被引入的顺序会严重影响结果。

约束优化可以看作是克服这些困难的一种方法。首先，所有的部分都被平等地同时对待。在图形相关中，两个部分同时相关，其中一个通常是组合部分。这是通过在由两个部分中事件的位置所跨越的二维空间中建立一条相关线(LOC)来实现的。相比之下，约束优化(restricted optimization, CONOP)在J维空间中构造了一个LOC，其中J是section的数量。因此，将各部分按顺序集成到解决方案中所造成的偏差消除了。

conop_data <- read.csv('C:\\Users\\dell\\Desktop\\北大R语言\\约束最优与图形相关\\conop_data.csv') #读取数据
head(conop_data[, 1 : 7]) #查看前7列前6行数据

##     name WangjiaWan.North.FAD WangjiaWan.North.LAD WangjiaWan.South.FAD
## 1 name_1                242.5                273.5                   NA
## 2 name_2                   NA                   NA                   NA
## 3 name_3                200.0                312.0                  100
## 4 name_4                207.0                222.5                  142
## 5 name_5                   NA                   NA                   NA
## 6 name_6                   NA                   NA                   NA
##   WangjiaWan.South.LAD Fenxiang.FAD Fenxiang.LAD
## 1                   NA          158          235
## 2                   NA           NA           NA
## 3                  159          100          331
## 4                  159          158          276
## 5                   NA          227          331
## 6                   NA           NA           NA

第一列name表示一系列物种（FAD指第一个出现，LAD指最后一次出现，顶上是不同的剖面）

min <- apply(conop_data[, -1], 1, function(x) min(x, na.rm = TRUE)) #选出每行最小值
max <- apply(conop_data[, -1], 1, function(x) max(x, na.rm = TRUE)) #选出每行最大值
mm <- max - min #计算增量(每行最大值-每行最小值)
min_mm <- matrix(c(min, mm),nrow=2, ncol=31, byrow=TRUE) #生产2行31列矩阵
barplot(min_mm, xlab='', ylab='Height In Composite Section', col=c('white', 'gray'), border=NA, ylim = c(0, 2000))
text(seq(0.8,37,1.2), max + 25, 1:31, cex = 0.6)

该图反映了不同类别的存在时间，如图中类别1的存在时间大概为150-750之间。

plot(conop_data$WangjiaWan.North.FAD, conop_data$WangjiaWan.South.FAD, xlim = c(200, 400), ylim = c(100, 250), pch = 19, xlab = 'WangjiaWan North', ylab = 'WangjiaWan South')
par(new = TRUE)
plot(conop_data$WangjiaWan.North.LAD, conop_data$WangjiaWan.South.LAD, xlim = c(200, 400), ylim = c(100, 250), pch = 3, xlab = '', ylab = '')
x <- c(conop_data$WangjiaWan.North.FAD, conop_data$WangjiaWan.North.LAD)
y <- c(conop_data$WangjiaWan.South.FAD, conop_data$WangjiaWan.South.LAD)
dd <- data.frame(x = x, y = y) #生产数据框
dd <- dd[complete.cases(dd), ] #删除有缺失值的行
lines(lowess(dd$x, dd$y), col = 'red', lwd = 2) #添加趋势线

横坐标是WangjiaWan North，纵坐标是WangjiaWan South。其中圆圈表示的是WangjiaWan.North.FAD数据，加号表示的是WangjiaWan.North.LAD数据。

3 本章汇总

参数	类别	功能
abline	函数	图中添加回归线
text	函数	图中添加文字
arrows	函数	图中添加箭头
lm	函数	建立回归方程
apply	函数	对行或列进行迭代操作
barplot	函数	绘制柱图

4 参考文献

[1] Westrop, S R.Trilobite systematics, biostratigraphy and palaeoecology of the lower ordovician komstad limestone and huk formations, southern Scandinavia[J].Earth-Science Reviews, 03/1997, 卷 42, 期 1-2

[2]Daniel Goldman,Stig M,Bergström,H. David Sheets,Carolyn Pantle.A CONOP9 composite taxon range chart for Ordovician conodonts from Baltoscandia: a framework for biostratigraphic correlation and maximum-likelihood biodiversity analyses[J].GFF:Volume 136, 2014 - Issue 2

[3]Junxuan Fan a,Qing Chen a, b, c,Michael J.Melchin b,H. David Sheets d.Quantitative stratigraphy of the Wufeng and Lungmachi black shales and graptolite evolution during and after the Late Ordovician mass extinction[J].Palaeogeography, Palaeoclimatology, Palaeoecology:Volume 389, 1 November 2013, Pages 96-114

源文件下载在线调试

返回：数据分析教程案例, R 语言