matlab 使用LOWESS滤波估计海洋热含量(OHC)的趋势曲线

[matlab 使用LOWESS滤波估计海洋热含量(OHC)的趋势曲线]
matlab 使用LOWESS滤波估计海洋热含量(OHC)的趋势曲线
估计海洋热含量（OHC）的长期趋势及其不确定性
结果展示以及解释：
中央趋势（红线）

趋势的 不确定性区间（灰带）

可用于：气候变化信号的显著性分析、热含量变化的统计置信范围估计

图片
画图思路解释：

用 LOWESS 平滑提取长期趋势
对原始 OHC 序列应用 LOWESS 滤波器，获得平滑的趋势曲线。

IAP_OHC_0_700m_LOWESS = smooth(IAP_OHC_year,IAP_OHC_0_700m,window_size,‘lowess’,2);

smooth(…, ‘lowess’) 是 局部加权回归（locally weighted scatterplot smoothing）

对每个时间点，利用它周围 window_size 年的观测点进行局部回归，得到平滑曲线

作用：消除高频波动，提取长期趋势（如气候变化趋势）

计算残差序列（原始值 - 平滑趋势）

residual = 原始值 - 平滑值;

图片
蒙特卡洛模拟（realization_LOWESS_new 函数）

这个函数做了以下事情：

拟合残差的 AR(1) 模型，获得参数 ϕ\phiϕ 和噪声分布

重复 realization_number 次：

用 AR(1) 模型模拟一条残差序列

加到原始 LOWESS 曲线上，形成一条“可能的”趋势曲线

对这条合成曲线再次应用 LOWESS 滤波，获得平滑模拟曲线

最终生成 realization_number 条趋势曲线的集合

这些曲线构成趋势的不确定性估计基础。

这里直接放到函数里了。简单的很：
直接调用函数即可：
[IAP_OHC_0_700m_LOWESS_realise]=realization_LOWESS_new(IAP_OHC_year,IAP_OHC_0_700m,realization_number,window_size);

代码：
.rtcContent { padding: 30px; } .lineNode {font-size: 12pt; font-family: “Times New Roman”, Menlo, Monaco, Consolas, “Courier New”, monospace; font-style: normal; font-weight: normal; }
clear;clc;close all;
% OHC_trend_with_uncertainty_by_MC_LOWESS.m
%
% 通过LOWESS滤波估计海洋热含量(OHC)的趋势曲线，
% 并基于残差序列构建AR(1)模型，使用蒙特卡洛方法生成多条平滑趋势模拟，
% 用于估计趋势的不确定性和置信区间。
%
% 输入：
% - IAP_OHC_year: 年份序列
% - IAP_OHC_0_700m: 0-700m海洋热含量时间序列
% - realization_number: 模拟曲线数量（如1000）
% - window_size: LOWESS滤波窗口宽度（年）
%
% 输出：
% - 多条模拟趋势曲线，可用于不确定性分析和可视化
%%% Number of realization
realization_number=1000;
window_size = 25; % 比如25年窗口，或者改成你想用的数值
%%% Load IAP OHC time seires (annual means)
load IAP_OHC_annual_record.mat IAP_OHC_0_700m IAP_OHC_year
IAP_OHC_0_700m_LOWESS = smooth(IAP_OHC_year,IAP_OHC_0_700m,window_size,‘lowess’,2);
[IAP_OHC_0_700m_LOWESS_realise]=realization_LOWESS_new(IAP_OHC_year,IAP_OHC_0_700m,realization_number,window_size);
figure(‘position’,[10 10 850 550],‘color’,‘w’);
title([‘OHC 0-700m time series, LOWESS smooth and 1000 realization’])
plot(IAP_OHC_year,IAP_OHC_0_700m_LOWESS_realise,‘color’,[0.8,0.8,0.8]);hold on
plot(IAP_OHC_year,IAP_OHC_0_700m,‘b’,‘linewidth’,2);hold on
plot(IAP_OHC_year,IAP_OHC_0_700m_LOWESS,‘r’,‘linewidth’,2);hold on
figure(‘position’,[10 10 850 550],‘color’,‘w’);
title([‘OHC 0-700m time series, LOWESS smooth and 1000 realizations’])
% 灰色：1000 条模拟趋势（取其中一条代表）
h1 = plot(IAP_OHC_year, IAP_OHC_0_700m_LOWESS_realise, ‘color’, [0.8, 0.8, 0.8]); hold on
% 蓝色：原始序列
h2 = plot(IAP_OHC_year, IAP_OHC_0_700m, ‘b’, ‘linewidth’, 2); hold on
% 红色：LOWESS 趋势
h3 = plot(IAP_OHC_year, IAP_OHC_0_700m_LOWESS, ‘r’, ‘linewidth’, 2); hold on
% 添加图例，h1(1) 只用一条灰线代表所有模拟趋势
legend([h1(1), h2, h3], ‘Monte Carlo realizations’, ‘Original OHC’, ‘LOWESS trend’, ‘Location’, ‘NorthWest’);
xlabel(‘Year’);
ylabel(‘OHC (0–700 m)’);
grid on;
box on;

参考：可点击阅读原文：
https://iap.cas.cn/gb/lunbo/202408/t20240805_7240685.html

(https://mp.weixin.qq.com/s/fF2O1vJePnWnYnt33fM2uQ)