如何消除Matlab搭建的DTS仿真系统曲线的锯齿失真？

我在Matlab中开发了一套分布式光纤温度传感（DTS）系统的数值仿真程序，整体结果已接近预期，但生成的曲线图中出现了锯齿状失真（蓝色圈注区域），而非平滑曲线，现寻求解决该问题的方案。

clc;
clear;
close all;

% 光纤参数
z = linspace(0, 10, 50000); % 10 km，共50000个采样点
dz = z(2) - z(1);           % 空间步长（km）

% 物理常数
h = 6.62607015e-34;         % 普朗克常数（J·s）
c = 299792458;              % 光速（m/s）
k = 1.380649e-23;           % 玻尔兹曼常数（J/K）
delta_sigma = 440 * 100;    % 拉曼频移，从440 cm⁻¹转换为m⁻¹
S_IE = h * c * delta_sigma / k; % 温度系数（K）

% 波长参数
lambda_p = 1550e-9;                      % 泵浦波长（m）
lambda_s = 1/(1/lambda_p - delta_sigma);   % 斯托克斯波长（m）
lambda_as = 1/(1/lambda_p + delta_sigma);  % 反斯托克斯波长（m）

% 温度分布（单位：°C）
T_C = 20 * ones(size(z));   
T_C(z >= 0.800 & z <= 1.050) = 50;      % 0.8-1.05km区域温度设为50°C
T_K = T_C + 273.15;                     % 转换为开尔文温度

% 拉曼背向散射效率
eta_as = 0.000015; % 反斯托克斯散射效率
eta_s = 0.00015;   % 斯托克斯散射效率

% 衰减系数（dB/km）
alpha_p_dB = -0.2;     % 泵浦光衰减
alpha_s_dB = -0.21;    % 斯托克斯光衰减
alpha_as_dB = -0.25;   % 反斯托克斯光衰减

% 基于波长的效率修正
eta_s_prime = eta_s * (lambda_p/lambda_s)^4;
eta_as_prime = eta_as * (lambda_p/lambda_as)^4;

% 平均次数与噪声参数
num_avg = 7500;          % 平均次数
noise_level = 0.05;      % 噪声水平

% 衰减系数差值（用于计算增量误差）
delta_alpha_per_km = abs(alpha_as_dB - alpha_s_dB);

% 初始化计算温度的累加变量
T_calculated_K_sum = zeros(size(z));

for n = 1:num_avg
    % 泵浦光（正向）、斯托克斯/反斯托克斯光（反向）的衰减
    att_p = 10.^(alpha_p_dB * z / 10);
    att_s = 10.^(alpha_s_dB * z / 10);
    att_as = 10.^(alpha_as_dB * z / 10);
    
    % 无噪声情况下的理论背向散射功率
    potencia_ps = 1 * eta_s_prime .* (att_p .* att_s);
    potencia_pas = 1 * eta_as_prime .* (att_p .* att_as) .* exp(-S_IE./T_K);
    
    % 添加与初始功率成正比的高斯噪声
    noise_ps = noise_level * potencia_ps(1) * randn(size(z));
    noise_pas = noise_level * potencia_pas(1) * randn(size(z));
    
    % 添加噪声后的功率（确保为正值）
    potencia_ps_noisy = abs(potencia_ps + noise_ps);
    potencia_pas_noisy = abs(potencia_pas + noise_pas);
    
    % 实测拉曼比值
    R_measured = potencia_pas_noisy ./ potencia_ps_noisy;
    
    % 校准常数
    C0 = eta_as_prime / eta_s_prime;
    
    % 矢量计算所有点的温度分布
    term_dB = (alpha_as_dB - alpha_s_dB) * z; % 累计衰减差值
    term_linear = 10.^(term_dB / 10);
    
    argument = C0 ./ (R_measured .* term_linear);
    argument(argument <= 0) = NaN; % 避免对非正数取对数
    
    T_calculated_K = zeros(size(z));
    T_calculated_K(:) = NaN; % 初始化
    
    valid_idx = ~isnan(argument);
    T_calculated_K(valid_idx) = S_IE ./ log(argument(valid_idx));
    
    % ---- 计算增量误差（原逻辑）----
    delta_alpha_acumulado = delta_alpha_per_km * z; % 累计衰减差值（dB）
    erro_temperatura_incremental = floor(delta_alpha_acumulado / 0.03); % 每0.03 dB对应1K误差
    
    % ------- 新增：利用50°C区域修正误差 -------
    idx_regiao_aquecida = (z >= 0.800 & z <= 1.050); % 高温区域索引
    
    erro_local = T_K(idx_regiao_aquecida) - T_calculated_K(idx_regiao_aquecida);
    % 计算局部误差的平均值（忽略NaN）
    erro_medio_local = mean(erro_local,'omitnan');
    
    % 将局部平均误差应用到整个温度分布
    T_calculated_K_corr = T_calculated_K + erro_medio_local;
    
    % 叠加修正后的分布与离散增量误差
    T_calculated_K_corr = T_calculated_K_corr + erro_temperatura_incremental;
    
    % 累加用于最终平均
    T_calculated_K_sum = T_calculated_K_sum + T_calculated_K_corr;
    
end

% 计算最终平均温度并转换为摄氏度
T_calculated_K_avg = T_calculated_K_sum / num_avg;
T_calculated_C_avg = T_calculated_K_avg - 273.15;

% 参考区域索引（50°C线圈）
idx_ref = (z >= 0.800 & z <= 1.050);

% 参考区域计算温度的统计值
media_calc_ref = mean(T_calculated_C_avg(idx_ref), 'omitnan');
desvio_calc_ref = std(T_calculated_C_avg(idx_ref), 'omitnan');

% 计算温度分布的Z-score归一化（以参考区域为基准）
T_calculated_C_norm = (T_calculated_C_avg - media_calc_ref) / desvio_calc_ref;

% 真实温度分布的Z-score归一化（使用计算值的统计量，避免真实温度方差为0导致的除零）
T_C_norm = (T_C - media_calc_ref) / desvio_calc_ref;

% 绘制归一化温度分布
figure;
plot(z, T_calculated_C_norm, 'r-', 'LineWidth', 1.5, 'DisplayName', '归一化计算温度（Z-score）');
hold on;
plot(z, T_C_norm, 'k', 'LineWidth', 2, 'DisplayName', '归一化真实温度（Z-score）');
xlabel('距离（km）');
ylabel('归一化温度（Z-score）');
title('基于参考线圈的归一化温度分布（Z-score）');
legend;
grid on;
xlim([0.5 2.0]);

内容的提问来源于stack exchange，提问作者LEO101