基于 R 语言的 MTM 多窗谱分析

多窗谱分析（Multi-Taper Method，MTM）作为一种低方差、高分辨率的谱分析方法，可用于对气温、降水量等时间序列数据的周期性变化分析，且已在气候研究领域中得到了广泛应用。本文将介绍如何在 R 语言环境下，借助相关 R 包实现 MTM 分析，并导出分析结果。

1 安装及加载 R 包

在开始分析之前，我们需要加载两个关键的 R 包：astrochron 和 openxlsx。在 R 包管理器中，输入以下命令即可完成安装与加载：

# 安装 R 包
install.packages("astrochron") # 主要用于天文年代学研究，可进行 MTM 多窗谱分析
install.packages("openxlsx") # 可用于读写 Excel 文件

# 在 R 语言程序中加载对应 R 包
library(astrochron) 
library(openxlsx) 

2 MTM 详细代码

2.1 读取数据

假设我们的时间序列数据为 Excel 文件 "data.xlsx" 中的 "Sheet1" 工作表，可利用 openxlsx 包中的 read.xlsx 函数读取数据，如下所示：

data <- read.xlsx("data.xlsx", sheet = "Sheet1")
str(data) # 查看数据结构

通过 str() 函数查看数据结构，确保数据正确读取并符合分析要求：

'data.frame':	122 obs. of  2 variables:
 $ Year: num  1901 1902 1903 1904 1905 ...
 $ Data: num  1017 940 1481 1170 1320 ...

上方输出结果表明，data 有 122 条数据，分为 Year 和 Data 两列。其中 Year 为年份，Data 为对应年份的数据。

2.2 执行 MTM 分析

借助 astrochron 包中的 mtm 函数，可轻松进行 MTM 分析。代码如下：

mtm_result <- mtm(data, ar1 = TRUE, output = TRUE)
str(mtm_result)

函数详细说明请参考文末的 astrochron CRAN 网站中的说明文档。

通过 str() 函数查看数据结构，结果如下所示：

'data.frame':	304 obs. of  8 variables:
 $ Frequency   : num  0.00164 0.00328 0.00492 0.00656 0.0082 ...
 $ Power       : num  129.3 107.5 89.9 88.1 102.3 ...
 $ Harmonic_CL : num  5.36 21.59 42.52 49.04 34.42 ...
 $ AR1_CL      : num  52.8 37.4 24.8 23.5 33.6 ...
 $ AR1_fit     : num  134 134 134 134 134 ...
 $ AR1_90_power: num  214 214 214 214 214 ...
 $ AR1_95_power: num  245 245 245 245 245 ...
 $ AR1_99_power: num  311 311 311 311 311 ...

上方输出结果表明，mtm_result 分为八列。其中后续制图中主要会用到频率 (Frequency)、谱值 (Power) 以及 90%、95%、99%检验结果 (AR1_90_power, AR1_95_power, AR1_99_power)。

2.3 将分析结果导出为 Excel 文件

利用 openxlsx 包中的 write.xlsx 函数，将 MTM 分析结果导出为 Excel 文件。

write.xlsx(mtm_result, "mtm_result.xlsx")

输出的表格可通过 Excel 或 Origin 等软件生成对应统计图。Origin 输出的结果图如下所示，其中周期值与 Frequency 数值互为倒数关系。

参考资料：

1. 江志红, 屠其璞, 施能. 多窗谱分析方法及其在全球变暖研究中的应用[J]. 气象学报, 2001(4): 480-490.
2. 许文锋, 张乐满, 范依捷, 等. 1470 年以来长江流域降水重建及其特征分析[J]. 地理科学, 2024, 44(11): 2029-2038.
3. CRAN: Package astrochron
4. CRAN: Package openxlsx