光谱数据处理：核回归平滑

做光谱建模时，原始反射率曲线几乎都会带有随机噪声。噪声来源很多，比如传感器误差、环境干扰、采集过程不稳定等。

如果直接拿这类数据建模，叶绿素 a 浓度反演模型通常会受到影响。一个常见且有效的做法，就是先做平滑，尽量压低噪声，同时保留真实的光谱变化趋势。

核回归平滑（Kernel Regression Smoothing）之所以常用于这类任务，核心原因是：

• 它是非参数方法，不依赖固定函数形式。
• 对非线性光谱特征的适应性更好。
• 相比简单移动平均，更不容易把有效细节“抹平”。

在水体叶绿素 a 遥感反演中，核回归平滑通常能提升光谱一致性，从而提高模型稳定性和跨数据集泛化能力^[1]。

1. 核回归函数（MATLAB）

首先先准备一个一维核回归函数 ksr（Nadaraya-Watson + 高斯核）：

▶

核回归函数ksr.m

function r=ksr(x,y,h,N)
% KSR   Kernel smoothing regression
%
% r=ksr(x,y) returns the Gaussian kernel regression in structure r such that
%   r.f(r.x) = y(x) + e
% The bandwidth and number of samples are also stored in r.h and r.n
% respectively.
%
% r=ksr(x,y,h) performs the regression using the specified bandwidth, h.
%
% r=ksr(x,y,h,n) calculates the regression in n points (default n=100).
%
% Without output, ksr(x,y) or ksr(x,y,h) will display the regression plot.
%
% Algorithm
% The kernel regression is a non-parametric approach to estimate the
% conditional expectation of a random variable:
%
% E(Y|X) = f(X)
%
% where f is a non-parametric function. Based on the kernel density
% estimation, this code implements the Nadaraya-Watson kernel regression
% using the Gaussian kernel as follows:
%
% f(x) = sum(kerf((x-X)/h).*Y)/sum(kerf((x-X)/h))
%
% See also gkde, ksdensity

% Example 1: smooth curve with noise
%{
x = 1:100;
y = sin(x/10)+(x/50).^2;
yn = y + 0.2*randn(1,100);
r=ksr(x,yn);
plot(x,y,'b-',x,yn,'co',r.x,r.f,'r--','linewidth',2)
legend('true','data','regression','location','northwest');
title('Gaussian kernel regression')
%}
% Example 2: with missing data
%{
x = sort(rand(1,100)*99)+1;
y = sin(x/10)+(x/50).^2;
y(round(rand(1,20)*100)) = NaN;
yn = y + 0.2*randn(1,100);
r=ksr(x,yn);
plot(x,y,'b-',x,yn,'co',r.x,r.f,'r--','linewidth',2)
legend('true','data','regression','location','northwest');
title('Gaussian kernel regression with 20% missing data')
%}

% By Yi Cao at Cranfield University on 12 March 2008.
%

% Check input and output
error(nargchk(2,4,nargin));
error(nargoutchk(0,1,nargout));
if numel(x)~=numel(y)
    error('x and y are in different sizes.');
end

x=x(:);
y=y(:);
% clean missing or invalid data points
inv=(x~=x)|(y~=y);
x(inv)=[];
y(inv)=[];

% Default parameters
if nargin<4
    N=100;
elseif ~isscalar(N)
    error('N must be a scalar.')
end
r.n=length(x);
if nargin<3
    % optimal bandwidth suggested by Bowman and Azzalini (1997) p.31
    hx=median(abs(x-median(x)))/0.6745*(4/3/r.n)^0.2;
    hy=median(abs(y-median(y)))/0.6745*(4/3/r.n)^0.2;
    h=sqrt(hy*hx);
    if h<sqrt(eps)*N
        error('There is no enough variation in the data. Regression is meaningless.')
    end
elseif ~isscalar(h)
    error('h must be a scalar.')
end
r.h=h;

% Gaussian kernel function
kerf=@(z)exp(-z.*z/2)/sqrt(2*pi);

r.x=linspace(min(x),max(x),N);
r.f=zeros(1,N);
for k=1:N
    z=kerf((r.x(k)-x)/h);
    r.f(k)=sum(z.*y)/sum(z);
end

% Plot
if ~nargout
    plot(r.x,r.f,'r',x,y,'bo')
    ylabel('f(x)')
    xlabel('x')
    title('Kernel Smoothing Regression');
end

这个函数适合把带噪声曲线平滑成趋势线，尤其适用于光谱曲线预处理。

2. ksr 的常见调用方式

最常用的是以下 4 种：

1. ksr(x, y)

• 自动估计带宽 h
• 默认在 100 个点上输出平滑结果
• 如果不接收输出，会直接画图

2. r = ksr(x, y)

• 返回结果结构体 r，不自动画图
• r 里有：
  • r.x：平滑后横坐标（等间距）
  • r.f：平滑后纵坐标
  • r.h：带宽
  • r.n：用于回归的数据点数（剔除 NaN 后）

3. r = ksr(x, y, h)

• 手动指定带宽 h（平滑强度）
• h 越大越平滑，细节越少
• h 越小越贴近原始数据，可能保留噪声

4. r = ksr(x, y, h, N)

• 手动指定带宽和输出点数 N
• N 是结果曲线采样点数，不是原始数据点数

3. 光谱数据中的典型用法

波长 + 反射率的简单示例：

x = wavelength;      % 例如 400:2500
y = reflectance;     % 与 x 等长
r = ksr(x, y, 15, numel(x));
plot(x, y, 'k.', r.x, r.f, 'r-', 'LineWidth', 1.5);

带宽选择可以按这个思路来：

1. 先跑自动带宽：r = ksr(x, y)
2. 再围绕自动结果微调 h（比如 0.5 x、1 x、2 x）
3. 如果希望输出点与原始波长一一对应，设 N = numel(x)

使用时还要注意：

1. x 和 y 必须等长
2. 这个函数会自动删除 NaN 数据点
3. h 和 N 必须是标量
4. 如果数据几乎没有变化，函数会报“回归无意义”错误（带宽过小/数据无变异）

4. 处理示例

先准备一段示例数据（350-357 nm）：

样本	350	351	352	353	354	355	356	357
1	0.006236933	0.006293467	0.006349733	0.006406133	0.0064808	0.0065554	0.006630067	0.006620733
2	0.0058478	0.005906133	0.005964333	0.006022467	0.006095667	0.006168733	0.006241933	0.006241667
3	0.005712733	0.0057696	0.005826667	0.005883667	0.0059576	0.0060316	0.006105467	0.0061084
4	0.005123643	0.005191357	0.005259357	0.005326857	0.005413929	0.005501071	0.005588143	0.005599071

下面对每组反射率进行核回归平滑，并导出为新的 Excel 文件：

▶

处理脚本示例

[num, ~, ~] = xlsread('data.xlsx');
wavelength = num(1, :);
reflectance = num(2:end, :);

% 创建一个结构体数组，用于保存每组数据的回归结果
results = struct('x', {}, 'f', {});

% 创建一个矩阵，用于保存整理后的回归结果
regression_results = zeros(size(reflectance));

figure; % 初始化图形
hold on; % 设置保持图形的状态

for i = 1:size(reflectance, 1)
    y = reflectance(i, :);
    
    % 进行内核平滑回归分析
    r = ksr(wavelength, y);
    
    % 保存回归结果到结构体数组和矩阵
    results(i).x = r.x;
    results(i).f = r.f;
    regression_results(i, :) = interp1(r.x, r.f, wavelength);
    
    % 绘制当前组的回归曲线和原始数据点
    plot(r.x, r.f, 'r')
    scatter(wavelength, y, 'b', 'filled')
end

hold off; % 完成所有绘图，释放图形状态

ylabel('反射率')
xlabel('波长')
title('反射率的内核平滑回归')

regression_results = [wavelength; regression_results];

% 导出回归结果到 xlsx 文件
output_file = 'data_hph.xlsx';
writematrix(regression_results, output_file); % 或者使用 writecell 函数

这段脚本完成了三件事：

1. 对每组反射率逐行平滑
2. 将平滑结果插值回原始波长坐标
3. 输出到 data_hph.xlsx，同时绘制原始点与平滑曲线做可视化对比

5. 小结

核回归平滑是光谱预处理中很实用的一步：

• 能在降噪和保留细节之间取得较好的平衡
• 参数可解释性强（尤其是带宽 h）
• 与后续建模流程衔接简单

如果你在做水体、土壤或植被等光谱反演任务，这个方法通常值得先跑一版基线结果。

对于光谱数据还可以考虑使用OpenSA。针对光谱分析过程中涉及的常用训练数据集分割、光谱预处理、波长选择和校准模型算法，建立了的完整的算法库，命名为 opensa（开放光谱分析）。

参考文献

1. 程春梅,韦玉春,王国祥,等.使用光谱平滑提高浑浊水体叶绿素a浓度估算模型的应用精度[J].遥感技术与应用,2013,28(06):941-948. ↩︎