做气候、海洋和大气资料分析时，我们经常会看到几个很像的说法：
“季节日循环”？
一年里的每一天，都有一个“气候态平均值”：
把这一年 365 天（或 366 天）的平均值连起来，就是所谓：
气候态季节日循环
“去除气候态季节日循环”？
最常见的数学理解其实很简单：
某一天的数据，减去参考时段内“同一天”的多年平均值。
比如你有一套日数据，参考期选的是 1991–2020 年。那么对于 2015 年 8 月 20 日 这一天，你要做的就是：
- 先取出 1991–2020 年所有的 8 月 20 日
- 再用 2015 年 8 月 20 日的实际值减去这个平均值
得到的结果，就是这一天相对于气候态季节背景的异常值，也就是常说的日距平。
“去除年循环”？
问了一下GPT，这里就要特别注意了：
它在不同文章里，可能代表不同处理方式，最常见的三种“去除年循环”含义
1. 减去逐月气候态
这是最常见的一种，尤其是处理月平均资料的时候。
做法是：
这样构造出来的是一个由 12 个月组成的平均年循环。
这种方法得到的是月距平。
2. 减去逐日气候态
有些文章虽然写的是“去除年循环”，但实际做法是：对日数据减去同一历日的多年平均值
也就是前面讲的“去除气候态季节日循环”。
所以很多时候，“去除年循环”和“去除季节日循环”在文献里会被混着说，但它们的时间分辨率其实并不一样。
3. 用谐波拟合一个平滑的年循环，再减去
这在热带波动的信号提取之前，非常之常见，比如说时空频谱分析，往往就会看到数据说明里写到：
seasonal cycle的前三个谐波被删除，这里我理解这个cycle就是指：年内随季节推进而重复出现的气候态变化。
一种比较常见于频谱分析或滤波前处理的方法，是不用直接的多年平均，而是用前几阶正弦、余弦函数去拟合一个平滑的季节周期，比如：
然后再把这个拟合出来的平滑周期减掉。
通过，进行傅立叶转换，可以发现，
从他的频率上也可以看出来，第一个谐波表示一年的周期，第二个谐波表示半年的周期，第三个表示1/3 年的周期，第四个表示1/4 年的周期，依次类推...
一般文献里会说去除年循环，定义为seasonal cycle的前三个谐波，就是去除：
要么将其减掉，要么在傅立叶转换的频率空间赋值为0，再转换回去。
有时候还得看你的数据分辨率，daily、monthly？你的目的是什么？
因为如果你处理的是 daily SST、daily precipitation、daily OLR、daily wind、daily MLD 这类日数据，只减去每个月的平均值通常还不够。
原因很简单：
- 同样是 8 月，8 月 1 日 和 8 月 31 日 的背景态本来就可能不一样。如果你只减去整个 8 月的平均值，那么月内真实存在的季节变化还会残留在结果里。
这时如果你希望更干净地提取：
那么按“同一历日”去除逐日气候态，通常会更合理。反正这很容易实现不是吗？
怎么理解“距平”？
很多人做完去季节循环之后，会得到一个叫 anomaly 的变量。这个 anomaly，其实就是：
原始值减去气候态平均值之后剩下的偏离部分
所以：
我理解就是“相对于平均季节背景的偏差”。
实际分析里最值得注意的几个问题
1. 参考期不同，距平也会不同
比如你选：
算出来的异常值都会有差别。
2. 日数据处理中要特别注意闰年
比如 2 月 29 日怎么处理，就是一个常见问题。
常见做法有两种：
Python 实现代码
在python中，xarray提供了一个分组的函数xarray.DataArray.groupby，可以很方便的实现相关的技计算， including “year”, “month”, “day”, “hour”, “minute”, “second”, “dayofyear”, “week”, “dayofweek”, “weekday” and “quarter”
- groupby('time.dayofyear')
NCL的话也很方便：
甚至xarray的官网提供了相关计算的示例代码，所以说多看官网还是有帮助的
下面直接让GPT生成对应代码
1. 去除气候态季节日循环
import numpy as npimport pandas as pdimport xarray as xrfrom typing import Optional, Tuple, UnionXrObj = Union[xr.DataArray, xr.Dataset]def_check_time_coord(obj: XrObj) -> None:if"time"notin obj.dims and"time"notin obj.coords:raise ValueError("输入数据必须包含 time 维或 time 坐标。")try: pd.DatetimeIndex(obj["time"].values)except Exception as e:raise ValueError("time 坐标必须能转换为 pandas.DatetimeIndex。") from edefremove_daily_climatology( obj: XrObj, ref_start: Optional[str] = None, ref_end: Optional[str] = None, drop_feb29: bool = True,) -> Tuple[XrObj, XrObj]:""" 去除气候态季节日循环（daily climatology）。 参数 ---- obj : xr.DataArray 或 xr.Dataset 含有 time 维的数据。 ref_start, ref_end : str, optional 参考期起止时间，例如 "1991-01-01", "2020-12-31"。 若不提供，则默认用整个时间段作为参考期。 drop_feb29 : bool 是否去除 2 月 29 日。 返回 ---- anomaly : 与原数据同类型 去除逐日气候态后的异常。 climatology : 与原数据同类型 逐日气候态，时间维会被替换为 mmdd。 """ _check_time_coord(obj) obj = obj.sortby("time")if ref_start isNoneor ref_end isNone: ref = objelse: ref = obj.sel(time=slice(ref_start, ref_end))if ref.sizes.get("time", 0) == 0:raise ValueError("参考期为空，请检查 ref_start/ref_end 是否正确。")if drop_feb29: mask_obj = ~((obj.time.dt.month == 2) & (obj.time.dt.day == 29)) mask_ref = ~((ref.time.dt.month == 2) & (ref.time.dt.day == 29)) obj = obj.sel(time=mask_obj) ref = ref.sel(time=mask_ref) obj_mmdd = pd.DatetimeIndex(obj.time.values).strftime("%m-%d") ref_mmdd = pd.DatetimeIndex(ref.time.values).strftime("%m-%d") obj = obj.assign_coords(mmdd=("time", obj_mmdd)) ref = ref.assign_coords(mmdd=("time", ref_mmdd)) climatology = ref.groupby("mmdd").mean("time") anomaly = obj.groupby("mmdd") - climatologyreturn anomaly, climatology
2. 去除月气候态年循环
defremove_monthly_climatology( obj: XrObj, ref_start: Optional[str] = None, ref_end: Optional[str] = None,) -> Tuple[XrObj, XrObj]:""" 去除月气候态年循环（monthly climatology）。 """ _check_time_coord(obj) obj = obj.sortby("time")if ref_start isNoneor ref_end isNone: ref = objelse: ref = obj.sel(time=slice(ref_start, ref_end))if ref.sizes.get("time", 0) == 0:raise ValueError("参考期为空，请检查 ref_start/ref_end 是否正确。") climatology = ref.groupby("time.month").mean("time") anomaly = obj.groupby("time.month") - climatologyreturn anomaly, climatology
3. 用谐波拟合
下面这个是之前我写的用来复习NCL函数，需要对daily数据 先 groupby 得到 climatology，即dayofyear，再用前 3 个谐波重建 climatology，最后如果需要，再从原始数据里减掉它，即去除前三个谐波。
import numpy as npdefextract_low_harmonics(data: xr.DataArray, n_harm: int = 3, dim: str = "dayofyear") -> xr.DataArray:""" Retain the mean and the first n_harm harmonics along a given dimension. """ axis = data.get_axis_num(dim)# FFT z_fft = np.fft.rfft(data.values, axis=axis) z_fft_filtered = z_fft.copy()# Zero out harmonics higher than n_harm slicer = [slice(None)] * z_fft.ndim slicer[axis] = slice(n_harm + 1, None) z_fft_filtered[tuple(slicer)] = 0# Inverse FFT smoothed = np.fft.irfft( z_fft_filtered, n=data.sizes[dim], axis=axis ).real# Return as DataArrayreturn xr.DataArray( smoothed, coords=data.coords, dims=data.dims, attrs={ **data.attrs,"smoothing": f"FFT low-pass: retain mean + first {n_harm} harmonics","long_name": f"{data.attrs.get('long_name', 'data')} smoothed by first {n_harm} harmonics" } )
4. 使用示例
日数据：去除气候态季节日循环
# da 为 daily DataArray，例如 daily SST / precipitation / OLRanom_daily, clim_daily = remove_daily_climatology( da, ref_start="1991-01-01", ref_end="2020-12-31", drop_feb29=True,)print(anom_daily)print(clim_daily)
月数据：去除月气候态年循环
anom_monthly, clim_monthly = remove_monthly_climatology( mon_da, ref_start="1991-01-01", ref_end="2020-12-31",)print(anom_monthly)print(clim_monthly)
日数据去除谐波
clim = da.groupby("time.dayofyear").mean("time")clim.dims# ('dayofyear', 'lat', 'lon')clim_h3 = extract_low_harmonics(clim, n_harm=3, dim="dayofyear")anom = da.groupby("time.dayofyear") - clim_h3
10个月宝宝每天需要喝多少奶粉？