基于 R 語(yǔ)言的 MTM 多窗譜分析

多窗譜分析（Multi-Taper Method，MTM）作為一種低方差、高分辨率的譜分析方法，可用于對(duì)氣溫、降水量等時(shí)間序列數(shù)據(jù)的周期性變化分析，且已在氣候研究領(lǐng)域中得到了廣泛應(yīng)用。本文將介紹如何在 R 語(yǔ)言環(huán)境下，借助相關(guān) R 包實(shí)現(xiàn) MTM 分析，并導(dǎo)出分析結(jié)果。

多窗譜分析（Multi-Taper Method，MTM）作為一種低方差、高分辨率的譜分析方法，可用于對(duì)氣溫、降水量等時(shí)間序列數(shù)據(jù)的周期性變化分析，且已在氣候研究領(lǐng)域中得到了廣泛應(yīng)用。本文將介紹如何在 R 語(yǔ)言環(huán)境下，借助相關(guān) R 包實(shí)現(xiàn) MTM 分析，并導(dǎo)出分析結(jié)果。

1 安裝及加載 R 包

在開(kāi)始分析之前，我們需要加載兩個(gè)關(guān)鍵的 R 包：astrochron 和 openxlsx。在 R 包管理器中，輸入以下命令即可完成安裝與加載：

# 安裝 R 包
install.packages("astrochron") # 主要用于天文年代學(xué)研究，可進(jìn)行 MTM 多窗譜分析
install.packages("openxlsx") # 可用于讀寫 Excel 文件

# 在 R 語(yǔ)言程序中加載對(duì)應(yīng) R 包
library(astrochron) 
library(openxlsx) 

2 MTM 詳細(xì)代碼

2.1 讀取數(shù)據(jù)

假設(shè)我們的時(shí)間序列數(shù)據(jù)為 Excel 文件 "data.xlsx" 中的 "Sheet1" 工作表，可利用 openxlsx 包中的 read.xlsx 函數(shù)讀取數(shù)據(jù)，如下所示：

data <- read.xlsx("data.xlsx", sheet = "Sheet1")
str(data) # 查看數(shù)據(jù)結(jié)構(gòu)

通過(guò) str() 函數(shù)查看數(shù)據(jù)結(jié)構(gòu)，確保數(shù)據(jù)正確讀取并符合分析要求：

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

上方輸出結(jié)果表明，data 有 122 條數(shù)據(jù)，分為 Year 和 Data 兩列。其中 Year 為年份，Data 為對(duì)應(yīng)年份的數(shù)據(jù)。

2.2 執(zhí)行 MTM 分析

借助 astrochron 包中的 mtm 函數(shù)，可輕松進(jìn)行 MTM 分析。代碼如下：

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

函數(shù)詳細(xì)說(shuō)明請(qǐng)參考文末的 astrochron CRAN 網(wǎng)站中的說(shuō)明文檔。

通過(guò) str() 函數(shù)查看數(shù)據(jù)結(jié)構(gòu)，結(jié)果如下所示：

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

上方輸出結(jié)果表明，mtm_result 分為八列。其中后續(xù)制圖中主要會(huì)用到頻率 (Frequency)、譜值 (Power) 以及 90%、95%、99%檢驗(yàn)結(jié)果 (AR1_90_power, AR1_95_power, AR1_99_power)。

2.3 將分析結(jié)果導(dǎo)出為 Excel 文件

利用 openxlsx 包中的 write.xlsx 函數(shù)，將 MTM 分析結(jié)果導(dǎo)出為 Excel 文件。

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

輸出的表格可通過(guò) Excel 或 Origin 等軟件生成對(duì)應(yīng)統(tǒng)計(jì)圖。Origin 輸出的結(jié)果圖如下所示，其中周期值與 Frequency 數(shù)值互為倒數(shù)關(guān)系。

參考資料：

1. 江志紅, 屠其璞, 施能. 多窗譜分析方法及其在全球變暖研究中的應(yīng)用[J]. 氣象學(xué)報(bào), 2001(4): 480-490.
2. 許文鋒, 張樂(lè)滿, 范依捷, 等. 1470 年以來(lái)長(zhǎng)江流域降水重建及其特征分析[J]. 地理科學(xué), 2024, 44(11): 2029-2038.
3. CRAN: Package astrochron
4. CRAN: Package openxlsx