0 前言

按照自己实现论文代码的思路，去研究作者的代码，找到自己的知识盲区和不足，提升编码技能。

本模块主要介绍代码实现思路。细节分析详见后续博客。

1、合成地震数据

利用波动方程？？？合成数据，也就是标签（Ground truth）。

代码实现在 generateHyperbolic.m中，

使用的函数为：D = hyperbolic_events(dt, f0, tmax, offset, tau, v, amp, snr, L);
% dt：采样间隔（秒）
%f0：中心频率（Hz）
%tmax：模拟的最长时间（秒）
%h：偏移矢量，单位为米
%tau，v，amp：截距矢量，均方根速度和每个线性事件的振幅，v以m/s为单位，tau以秒为单位）
%snr：信噪比（清洁的最大振幅信号/噪声的最大幅度）
%L：随机噪声是L个样本的平均值

论文中参数设置如下：

dt = 2./1000;
tmax = 2.;
n = 100;
offset = (-n:n)*10;
tau = [.5, .8, 1., 1.4];
v = [1700, 1800, 2000, 2300];
amp = [.4, .4, .6, .5];
f0 = 30;
snr = Inf;
L = 20;

生成的地震数据如下图：

但参数的设置是如何影响图像中曲线形态的，以及和实际地震数据采集过程中的，地震道之间的距离、采样间隔等实际指标是如何对应的，不明确？

2、生成含噪声的下采样的地震数据

输入：合成地震数据

输出：含噪声的下采样地震数据

数据预处理包括四个步骤：归一化、加噪声、下采样、预插值。

2.1 归一化

1）方法：

Min-Max Normalization
x' = (x - X_min) / (X_max - X_min)

2）作用：

加速模型收敛，节约训练时间

2.2 加噪声

此处为随机噪声。噪声水平[0,255]。

2.3 下采样

1）生成下采样模板

这里有两种方式可选择，按照指定模板下采样和自定义下采样。这里介绍自定义下采样。

这里的下采样可以随机去行、随机去列、有规律去列、无规律去列。

这里采用的是无规律去列。

如：mask = projMask(D, Ratio, sampleType)

参数 D 为输入图像，Ratio 为保留原图像的比率，sampleType是采样类型，这里选择的是iregc（无规律去除列数据）。

实现代码如下：

index = randperm(n); % n为原图像列的数目，将n个整数随机排列
sample = index(1:fix(n*r)); % fix是向下取整的意思，r是保留的比例 
mask(:, sample) = 1;

最终得到只含有0 和 1 的mask，0表示原图像对应位置的值去除，1表示保留。

2）得到下采样数据

将原始数据和下采样模板点乘可得到下采样数据，注意，这里还做了一个小处理，用原图像的均值将缺失值补齐，最起码缺失部分看起来没有那么突兀。

% Down-sampling input.
input = nlabel.*mask;   % nlabel为加噪声的地震数据
input(mask==0) = mean(nlabel(:));  % 将缺失值用均值补齐

2.4 预插值

作用：预插值可以提高效率和准确性，就像神经网络初始化权重时不是随机初始化权重，而是使用某种方法进行初始化。

这里用的是Shepard插值方法，又称为与距离成反比的加权法，其基本实现是将插值函数定义为各数据点函数值的加权平均，权函数定义为与距离成反比。

什么类型的插值算法不重要，重要的是得有这一步，好的初始化很重要！

3、实现IRCNN-FPOCS的迭代过程

3.1 论文算法的整体描述如下：

%  论文算法步骤5   
maxv = 30;
epsilon = 10;
LambdaS = maxv * exp(((0:totalIter-1) * (log(epsilon) - log(maxv))) / (totalIter-1));  

for itern = 1 : totalIter
    %%%%%%%%%%%%%%%
    
    d_old = output;
    
    % 对应论文算法的步骤3、4、7，这里认为 alpha=1，tt为论文算法步骤3里面的st.
    % 步骤7作用有点类似于网络训练中的 optimize.zero_grad() 放在前面后面都可以
    tt = next_t(t);  
    beta = (t-1)/tt;
    output = mask.*input + (1 - mask).*(output + beta * (output - d_old));
   
    d_old = output;
    %%%%%%%%%%%%%%%%%%%%%
    
    if ns(itern+1) ~= ns(itern)
        [net] = loadmodel(LambdaS(itern), CNNdenoiser);    
        net = vl_simplenn_tidy(net); %修复不完整或过时的网络
        if useGPU
            net = vl_simplenn_move(net, 'gpu');
        end
    end
    
    for k = 1 : inIter
        res    = vl_simplenn(net,output,[],[],'conserveMemory',true,'mode','test');
        output = output - res(end).x;
    end
end

论文步骤5的目标是求噪声方差 $\sigma _{t}$ ，其最大值设置为30，最小值设置为10，代码中的LambdaS是采用指数方式迭代30次，产生30个 30和10之间的值，作为迭代过程中IRCNN 去噪模型的参数。