基于ASAMC算法的地质层中氧元素含量的变点问题

A Change-Point Problem of Oxygen Element in Geological Strata Is Based on ASAMC Algorithm

摘要: 变点是地质学的热门问题，多用于微量元素检测，具有很强的应用价值。由Kim，J. and Cheon，S.在2010年提出的退化后的随机逼近蒙特卡罗算法(Annealing Stochastic Approximation Monte Carlo，简称ASAMC)；ASAMC算法既可以检测变点个数，又可以检测变点所在的位置。本篇文章基于ASAMC算法，首先对地质层中氧元素含量进行正态性检验，随后，研究不同温度、不同季节下，土壤中氧元素含量的变化情况。最后，利用R软件，我们发现高温多雨的天气会较明显的影响土壤中氧元素的含量。

Abstract: The change-point is a hot issue in geology, and it is widely used in tracing element detection, and it has strong application value. The annealing stochastic approximation Monte Carlo (ASAMC) algo-rithm proposed by Kim, J. and Cheon, S. in 2010 can detect both the number of change points and the location of change points. In this paper, based on the ASAMC algorithm, the normality test of oxygen content in the geological layer is carried out firstly. Then, the change of oxygen content in the soil under different temperatures and seasons is studied. Last, using R software, we find that the high temperature and rainy weather will obviously affect the oxygen content in the soil.

1. 引言

在大数据背景下，各种分布的极限分布均为正态分布。因此，考察正态分布均值变点模型，对于研究大数据下各种分布的均值变点模型来说都有很大的应用价值，正态分布均值变点已成为解决其他分布均值变点较为普遍的方法。这篇文章基于ASAMC算法，首先对地质层中氧元素含量进行正态性检验，随后，利用R软件研究不同温度、不同季节下，土壤中氧元素含量的变化情况，我们发现高温多雨的天气会较明显的影响土壤中氧元素的含量。

ASAMC算法是通过SAMC算法演变而来，这种方法最早出现在Kim，J. and Cheon，S (2010) [1] 的文章中，随后，许欢在2016年将ASAMC算法用于研究气象变化规律 [2]，成守尧等在2022 [3] 年采用Wilcoxon秩检验法处理变点数据，程夏2022 [4] 年将变点理论与混合分布模型相结合处理数据，文章 [5] [6] [7] 将变点理论与地质数据相结合，利用变点理论的处理地形、水文等自然现象。本篇利用ASAMC算法与气象数据相结合，找出数据链中的变点。本篇文章将ASAMC算法用于探寻地质层中氧元素含量变化的成因，利用统计方法处理地质数据，减少了地质调研的时间并结合乌鲁木齐市水磨沟区地质特点，分析氧元素含量与哪些因素有关。

2. 正态分布均值变点模型

变点理论是上世纪五十年代提出的，具有很强的应用价值。变点模型在工业、金融方面应用较广，变点主要分为以下三种，参数变点、概率变点、位置变点，其中位置变点定义为：

假设存在一个数据集，每个数据观测值相互独立，如果在某一时刻，模型中的某个或某些变量突然发生了变化，及存在一个时间点，在该点之前，数据集符合一个分布，在该点之后，数据集符合另一个分布，则该点为该数据集的位置变点 [8]。

2.1. 正态分布均值变点模型

假设 $\left\{Y_i,i=\left(1,2,\cdots ,n\right)\right\}$ 是服从正态分布的随机变量序列。其中，随机变量序列存在k个均值变点，即含有k个变点的正态分布均值变点模型如下：

$Y_i~\{\begin{array}{l}N\left({\mu }_1,{\sigma }^2\right),i=c_0,2,\cdots ,c_1\\ N\left({\mu }_2,{\sigma }^2\right),i=c_1+1,c_1+2,\cdots ,c_2\\ ⋮\\ N\left({\mu }_k,{\sigma }^2\right),i=c_{k-1}+1,c_{k-1}+2,\cdots ,c_k\\ N\left({\mu }_{k+1},{\sigma }^2\right),i=c_k+1,c_k+2,\cdots ,c_{k+1}\end{array}$

其中 $c_i\left(i=1,2,\cdots ,n\right)$ 表示变点位置，其中 $c_0=1$， $c_{k+1}=n$， $\sigma$ 为已知数，记 $\theta =\left(k,c_0,\cdots ,c_{k+1},{\mu }_1,\cdots ,{\mu }_{k+1}\right)$ 为该模型参数向量，则正态分布均值变点序列似然函数为

$L\left(Y|\theta \right)={\displaystyle \underset{j=c_0}{\overset{c_1}{\prod }}{f}_1\left(y_j|{\mu }_1,{\sigma }^2\right)}\cdot {\displaystyle \underset{j=c_1}{\overset{c_2}{\prod }}{f}_2\left(y_j|{\mu }_2,{\sigma }^2\right)}\cdots {\displaystyle \underset{j=c_k+1}{\overset{c_{k+1}}{\prod }}{f}_{k+1}\left(y_j|{\mu }_{k+1},{\sigma }^2\right)}$

其中 $f_i=\frac{1}{\sqrt{2\pi }\sigma }\exp \left\{\frac{-{\left(y-{\mu }_i\right)}^2}{2{\sigma }^2}\right\},i=1,\cdots ,k+1$。

2.2. 各种分布与正态分布之间的联系

在大数据背景下，样本均值服从的分布可以转成正态分布，因此，正态分布均值变点在大数据背景下有较强的应用价值。

辛钦大数定理 [9]：设 $X_1,X_2,\cdots$ 是相互独立，服从同一分布的随机变量序列，且具有数学期望

$E\left(X_k\right)=\mu \left(k=1,2,\cdots \right)$。作前n个变量的算术平均 $\frac{1}{n}{\displaystyle \underset{k=1}{\overset{n}{\sum }}{X}_k}$，则对于任意 $\epsilon >0$，有

$\underset{n\to \infty }{lim}P\left\{\left|\frac{1}{n}{\displaystyle \underset{k=1}{\overset{n}{\sum }}{X}_k}-\mu \right|<\epsilon \right\}=1$。

即序列 $\bar{X}=\frac{1}{n}{\displaystyle \underset{k=1}{\overset{n}{\sum }}{X}_k}$ 依概率收敛于 $\mu$， $\bar{X}\stackrel{P}{\to }\mu$。(n趋向于正无穷)

独立同分布的中心极限定理 [9]：设随机变量 $X_1,X_2,\cdots ,X_n,\cdots$ 相互独立，服从同一分布，其具有数学期望和方差： $E\left(X_k\right)=\mu ,D\left(X_k\right)={\sigma }^2>0\left(k=1,2,\cdots \right)$，则随机变量之和的标准变化量为

$Y_n=\frac{{\displaystyle \underset{k=1}{\overset{n}{\sum }}{X}_k}-E\left({\displaystyle \underset{k=1}{\overset{n}{\sum }}{X}_k}\right)}{\sqrt{D\left({\displaystyle \underset{k=1}{\overset{n}{\sum }}{X}_k}\right)}}=\frac{{\displaystyle \underset{k=1}{\overset{n}{\sum }}{X}_k}-n\mu }{\sqrt{n}\sigma }$

的分布函数 $F_n\left(x\right)$ 对于任意x满足

$\underset{n\to \infty }{\lim }{F}_n\left(x\right)=\underset{n\to \infty }{\lim }P\left\{\frac{{\displaystyle \underset{k=1}{\overset{n}{\sum }}{X}_k}-n\mu }{\sqrt{n}\sigma }\le x\right\}={\displaystyle {\int }_{-\infty }^x\frac{1}{\sqrt{2\pi }}{\text{e}}^{-\frac{{{\displaystyle t}}^2}{2}}\text{d}t}=\Phi \left(x\right).$

3. 实例分析

3.1. ASAMC算法

ASAMC算法既可以检验正态分布均值变点位置又可以检验变点个数，运用ASAMC算法，一般先将数据正态化，将服从其他分布的随机变量序列转化为服从正态分布的随机序列；再利用ASAMC算法检测正态分布均值变点的位置和个数。

3.2. ASAMC算法步骤 [1]

(1) 样本数据正态化

随机产生一组样本数据，记作 $Z=\left(z_1,z_2,\cdots ,z_n\right)$，检验样本数据 $Z=\left(z_1,z_2,\cdots ,z_n\right)$ 是否服从正态分布，具体做法如下：画出样本数据的散点图，观察这些点的分布特点，是否在一条直线附近波动，若存在这样的直线，我们就认为该样本数据服从正态分布；反之，样本数据不服从正态分布，要将样本数据进行正态化处理。

(2) 设定一个新的变量 $X_1=\left(x_1,x_2,\cdots ,x_n\right)$，随机产生k个变点，( $k\le \left[\frac{n}{10}\right]-1$ )，从1到n中随机选取

k个位置记为 $c_1,c_2,\cdots ,c_k$，其中 $1<c_1<c_2<\cdots <c_k<n$，假定，处于变点位置处的 $x_i=1,i=c_1,c_2,\cdots ,c_k$，其余位置为0。

(3) 定义相对样本频率 $\frac{n_i}{n}\cdot 100\%$ 来代替概率，迭代过程具体分为以下步骤：

(3.1) 根据样本数据的对数后验分布概率密度函数 $U\left(x\right)=\log p\left(x^k|Z\right)$，我们将样本空间 $\chi$ 进行划分，

分成m份，即 $\begin{array}{l}{E}_1=\left\{x\in \chi :U\left(x\right)<u_1\right\}\\ E_2=\left\{x\in \chi :u_1\le U\left(x\right)<u_2\right\}\\ ⋮\\ E_m=\left\{x:U\left(x\right)\ge u_{m-1}\right\}\end{array}$，第t次迭代时样本空间记为： ${\chi }^{\left(t\right)}={\displaystyle \underset{i=1}{\overset{\varpi \left(U_{\min }^{\left(t\right)}+\Delta \right)}{\cup }}{E}_i}$

(3.2) 根据贝叶斯分析理论，从建议分布 $q\left(x,y\right)$ 中随机抽取样本y

(3.3) 计算接受概率 $\min \left(1,r\right)$，其中

$r=C\frac{\varphi \left(y\right)q\left(y,x_t\right)}{\varphi \left(x_t\right)q\left(x_t,y\right)}$

其中 $q\left(x,y\right)$ 是定义在样本空间 $\chi$ 上的函数，C是与参数无关的常系数

(3.4) 利用R软件在 $\left(0,1\right)$ 上的均匀分布中抽取随机数u，

$x_{t+1}=\{\begin{array}{l}y,u\le \min \left(1,r\right)\\ x_t,u>\min \left(1,r\right)\end{array}$

(3.5) 多次重复(3.1)~(3.4)，找出不同的U和x的值。

(4) 若满足 $U\left(x_{t+1}\right)<U_{\max }^{\left(t\right)}$，则 $U_{\max }^{\left(t+1\right)}=U\left(x_{t+1}\right),x_{\max }^{\left(t+1\right)}=x_{t+1}$

(5) 输出 $U_{\max }^{\left(t+1\right)},x_{\max }^{\left(t+1\right)}$，做出相关数据的频数直方图。

(6) 观察频数直方图，若各组间有显著差别，则认为样本序列存在变点，反之不存在变点。

(7) 若模型存在变点，找出最大的 $U_{\max }^{\left(t+1\right)},x_{\max }^{\left(t+1\right)}$， $x_{\max }^{\left(t+1\right)}$ 的值表示变点的位置。

3.3. 运用ASAMC算法检测土壤中氧元素含量

变点理论在地质学中应用广泛。ASAMC算法既可以检测出正态分布均值变点个数，又可以检测出变点所在的位置。这篇文章将ASAMC算法应用于地质层中氧元素含量的变化上，探究地质层中元素含量与哪些因素有关。下面我们以乌鲁木齐水磨沟区新疆师范大学温泉校区的土地为研究对象(数据来源于新疆师范大学地理科学实验室)。

为了对已有数据进行分析处理 [10]，我们进行散点图处理，分析温度，土壤中氧元素含量的变化情况。图1表示从2014年1月到2017年5月新疆师范大学温泉校区温度变化情况；图2表示近年的温度变化导致地质层中氧元素含量变化情况，如下图所示：

Figure 1. Change of temperature

图1. 温度变化情况

Figure 2. Change of oxygen content

图2. 氧元素变化情况

观察图1、图2，我们发现地质层中氧元素含量大致呈现正态化，并且地质层中氧元素含量与2014年1月到2017年5月温度变化呈正相关，我们以地面为参考平面，测量土壤所在位置的距离(即该土壤层距离地面的最短垂直距离)，随着距离的加深，氧元素的含量在降低。考察温度和土壤中氧元素含量的关系见下表1：

Table 1. Change of oxygen content

表1. 氧元素变化情况

通过对表1的研究，地质层中氧元素含量与月平均温度的关系，通过对数据的研究，我们发现数据中存在变点，大多数据中只含有一个变点，因此，我们发现变点位置大致为24，即年，即2016年6月月平均温度对土壤中氧元素的含量的影响发生明显的变化，高温会使得土壤中氧元素含量增加。

下面研究季节对土壤中氧元素含量的影响，我们将研究重点放在春夏两个季节，见图3。

Figure 3. Changes of oxygen content in soil in spring and summer

图3. 春夏两季土壤中氧元素含量变化情况

根据图3中的最大对数概率值可以得到土壤中氧元素含量变化情况，其中夏季氧元素含量明显增加，即变点一般产生在夏季，由于夏季高温炎热，降水量大。因此，我们可以发现温度和季节都会影响土壤中氧元素的含量。

这篇文章通过ASAMC方法检验出的乌鲁木齐水磨沟区的新疆师范大学温泉校区的土地中氧元素含量变化情况。

新疆乌鲁木齐的夏季高温、多雨，使得土壤中氧元素含量升高，符合新疆气候特点，分析图2、图3我们发现，影响氧元素含量变化与新疆夏季气候有关，因此利用ASAMC算法处理数据具有较强的实用价值。

基金项目

石河子大学自主立项科研项目(ZZZC202032B)。