选题背景与意义

微分方程是数学中一个重要且广泛应用的领域，涉及到描述变化和相互关系的方程。它是一种包含导数或微分的方程，常用于自然现象建模以及解决科学和工程领域中的问题。微分方程研究是数学中的一个重要分支，涵盖了广泛的领域和应用。通过研究微分方程，我们可以理解和预测自然现象的行为，以及解决科学和工程中的实际问题。同时，微分方程的研究也促进了数学的发展和数学工具在其他学科中的应用。

常微分方程被广泛应用于物理、生物、化学、经济学等各个领域。在许多情况下，观测数据是已知的，而其中隐含的微分方程仍然难以捉摸。因此，数据驱动的常微分方程(或动力系统)发现和反演是一个重要的研究方向。

微分方程反演问题是相对于微分方程的求解而言的：

• 微分方程的求解通常是一个正向问题，即给定初始条件，求解关于未知函数的方程；
• 微分方程反演问题是指根据已知的结果，推导出产生这些结果的微分方程，是从结果向方程的逆向推导和反演。

在数学上，微分方程的反演问题具有一系列的数学挑战和困难：

• 反问题的病态性：微小扰动在反问题的输出中会导致较大的误差，因此需要稳定性分析和正则化方法来解决数值计算中的误差；
• 存在非唯一解情况：需要利用先验信息、约束条件和经验知识等来进行问题的约束和规定，以得到合理的解。

微分方程反演问题研究在数学理论和实际应用中具有重要意义。通过解决微分方程的反演问题，我们可以从有限的观测数据中重建和推断未知的边界条件、初始条件或未知函数，从而深入理解系统的行为和性质，并提供科学、工程和医学等领域的相关应用。

随着计算机技术的发展，数值计算和机器学习在微分方程反演问题研究中的应用也呈现出迅猛发展的趋势。传统的微分方程推断方法依赖于领域专家的知识和经验，而机器学习方法，如神经网络、深度学习和遗传算法等，可以自动从大量数据中学习模式和关系，从而更准确地推断出微分方程模型。如今我们可以获取到大规模的、高维度的数据集，这为从数据中推断微分方程提供了更多的信息和挑战。而传统的模型推断方法在处理大数据集时可能受到维度灾难和计算复杂度的限制。机器学习的强大计算能力和自适应建模能力为解决这些问题提供了新的思路。

利用机器学习从数据中反演微分方程的方法可以分为两个方面：

• 在模型选择方面，机器学习可以通过自动搜索和学习合适的微分方程模型，利用测度函数或结构特征来评估不同模型的性能，并选择最适合数据的微分方程模型；
• 在参数估计方面，机器学习可以利用大数据集中的样本来优化微分方程模型的参数，以最小化模型与实际数据之间的误差。

利用机器学习从数据中反演微分方程是微分方程和机器学习交叉领域的重要研究方向。它利用机器学习的能力和大数据的优势，可以更准确地建模和预测复杂的动力学系统，推进科学研究和实际应用的发展。

针对常微分方程，目前常用的反演方法有如下几种：

• 高斯过程：将高斯过程置于状态函数之上，然后通过最大化边际对数似然性从数据中推断参数，这种方法适用于解决高维问题，但是对系统的形式有限制，并且仅用于估计系统的参数。
• 符号回归：创建并更正与观测数据相对应的符号模型，为控制函数提供更具表现力的函数形式，但缺点是对于大型系统来说计算成本高，并可能容易过拟合。
• 稀疏回归：找到候选基函数的稀疏组合来近似控制函数，其系数由最小二乘法或贝叶斯回归确定。这种方法提供系统的明确形式，不需要太多先验知识，但是依靠适当的候选函数；对于没有简单或稀疏表示的复杂动力系统来说可能是低效的。
• 统计学习：通过最小化经验误差来学习系统在某个假设空间中的相互作用核，避免维度诅咒，可以发现非常高维度的系统，但是仅适用于具有交互作用的核函数的动力系统。
• 物理信息神经网络（PINN）：一种新的深度学习方法，用于解决非线性偏微分方程的正问题和反问题。通过在前馈神经网络中嵌入由偏微分方程描述的物理信息，在不需要标签数据的情况下，将 PINN 训练为偏微分方程的近似解的代理模型。受此启发，进一步发展出了多步神经网络LMNets，这本质上是使用给定相点上的流映射提供的信息来反演动力系统种的未知控制函数f的一种逆过程。

在数值分析中，发展高阶方法是许多应用中的一个重要课题。传统上，在求解动力系统时，高阶离散化技术，如线性多步法和龙格-库塔方法已经得到了发展。近年来，线性多步法也被用于动力系统的发现。随着使用深度学习的发现取得令人满意的进展，对它在动力系统发现的理论理解也在进一步发展。

基于以上所陈述的微分方程反演问题的研究现状，针对无显式方程的非线性系统这一在真实物理世界中更为广泛且常见的情形，我们希望基于高阶离散化技术线性多步法，结合多步神经网络LMNets这一深度学习方法，拟合和反演出数据背后的动力学物理规律。

算法框架与原理

研究对象——常微分方程（组）

首先明确我们的研究对象——常微分方程：
$$
\begin{aligned}&\frac{dx(t)}{dt}=f(x(t),t),\&x(t_{0})=x_{0}\end{aligned}
$$
上述式子是一个最简单的二维常微分方程示例，其中x(t) ∈ R是关于时间t的函数，函数关系由一个微分方程刻画，且给出了一定的初值条件以唯一确定函数x的表达式，从而能够真实地刻画物理过程随时间变化的发展情况。

事实上大多数时候我们需要观测的物理量远远不止x这一个，因此为了后续更容易推广到高维情形（后续无特别说明，均围绕自治常微分方程展开讨论），我们将上述的非自治常微分方程转化为向量化的自治常微分方程：
$$
\begin{aligned}&\frac{d\boldsymbol{u}(t)}{dt}=\hat{f}(\boldsymbol{u}(t)),\&\boldsymbol{u}(t_{0})=\boldsymbol{u}{0}\end{aligned}
$$
其中：
$$
u=(x,y)\in\mathbb{R}^{D+1},\hat{f}=(f,1)^{T},u{0}=(x_{0},t_{0})^{T}
$$
在反演问题中，已知在若干时刻节点t_1、t_2、…、t_n的解x的数据，我们一般不会希望直接解出x(t)的表达式，而是通过确定光滑函数f，从而利用给定的数据反演出微分方程。这样的传统是从简单的反演情形推广并沿用的，因为比起得到一个精确的函数表达式，我们更希望探究这个表达式背后的物理规律，而这往往离不开微分方程，因此保留微分方程的形式具有一定的必要性。

多步神经网络

线性多步法是用于微分方程数值求解的常用方法。传统上，它们用于求解给定常微分方程的解，可以称为微分方程的正问题。但线性多步法也可以用于反演给定解的原微分方程，属于微分方程的反问题，尤其是将线性多步法的经典数值方法与神经网络相结合，例如多步神经网络。

线性多步法

线性多步法是求解常微分方程数值解的一种常见方法。随着计算机技术的发展，线性多步法得到更加广泛的应用。人们逐渐发现，在某些情况下，线性多步法可以提供更高的数值精度和数值稳定性。此外，线性多步法可以与其他数值方法（如龙格-库塔法）结合使用，互补彼此的优点。

对于上述常微分方程（非自治），设其中结点总数为N，时间间隔h为常数:
$$
h = t_{n + 1} - t_n, n = 1, …, N - 1
$$
则第n个结点的步长为M的线性多步法有以下公式：
$$
\sum_{m=0}^{M}[\alpha_{m}x_{n-m}+h\beta_{m}f(x_{n-m},t)]=0,n=M、\ldots、N, \alpha_{_M}\neq0
$$
使用线性多步法求解常微分方程时，必须选择初始的步长M值，以及确定系数α、β的不同方法，然后可以依次计算n≥M的近似值x_n。

事实上，根据确定系数的方法不同，线性多步法也分为多种类型，本项目中主要使用如下的三种常见的线性多步法：

• Adams-Moulton 法：作为一种隐式方法，需要通过迭代或其他数值求解技术来解决每个时间步长的方程组。它在求解非刚性和刚性常微分方程时都表现良好，并且具有较高的数值稳定性，其精度随步长的减小而提高，但响应地可能会产生更高的计算代价。对于高阶常微分方程，Adams-Moulton 法需要结合相应的差分格式将高阶常微分方程转化为一阶方程组进行求解。

$$
\sum_{m=0}^M\alpha_mx_{n-m}=h\sum_{m=0}^M\beta_mf(x_{n-m})
$$

• Adams-Bashforth法：作为一种显式方法，在求解非刚性和刚性常微分方程时都表现良好，计算效率较高，并且容易实现，但可能会收到稳定性条件地限制，其精度随步长的减小而提高。对于高阶常微分方程，Adams-Bashforth法需要结合相应的差分格式将高阶常微分方程转化为一阶方程组进行求解。

$$
\sum_{m=0}^M\alpha_mx_{n-m}=h\sum_{m=0}^M\beta_mf(x(t_{n-m}))
$$

• 后向微分公式法(Backward Differentiation Formula, BDF)：同样是一种隐式方法，但与Adams-Moulton法不同，后向微分公式法是通过解一个非线性方程组来得到当前时间步长的数值解。这个方程组可以用牛顿迭代等数值求解技术来解决。该方法具有较好的数值稳定性，其精度依赖于使用的插值多项式的阶数。一般来说，其高阶形式可以提供更高的精度，但也会增加计算的复杂性。它适用于求解非刚性和刚性常微分方程，并且在稳定性和数值精度方面通常有良好的表现。

$$
\sum_{m=0}^{M}\alpha_{m}x_{n-m}=h\beta f(x(t_{n}))
$$

除此之外，即使采用同一种系数确定方式，选择的步长M不同，得到的系数α、β也有所不同。具体而言，本项目中针对以上三种线性多步法，分别选取步长M=2/3/4的三种情况进行实验测试，下面列出各方法对应不同步长时的系数情况：

线性多步法的优点在于其高阶的精度和较小的计算代价，可以有效地逼近微分方程的数值解。但线性多步法可能会受到稳定性和初始条件的限制，选择适当的步长和求解方法非常重要。

多步神经网络算法

从线性多步法的公式中得到启发，可以用神经网络去反演常微分方程右端的f。通过相等时间间隔h的数值解x的数据训练神经网络，该神经网络的参数可以通过使以下均方误差损失函数MSE最小化来训练得到：
$$
MSE=\frac{1}{N-M+1}\sum_{n=M}^N\parallel y_n\parallel^2
$$
其中：
$$
y_n=\sum_{m=0}^M\left[\alpha_mu_{n-m}+h\beta_mf_{NN}(u_{n-m})\right],n=M,\ldots,N
$$
式中f_NN表示神经网络对函数f的近似，y_n也称为线性多步法（LMM）残差。

在Python中编写函数train()以封装多步神经网络的训练流程：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

# 模型训练 
def train(u_tensor, model, loss_func, h, M): 
    # 参数说明：训练集u_tensor,待训练模型model,损失函数loss_func,时间步长h,步数M,最大训练次数EPOCH
    optimizer = torch.optim.Adam(model.parameters(), lr=LR)
    loss_history = []
    
    for epoch in range(EPOCHS):          
        model.train()        
        batch_u = u_tensor.to(device)        
        train_loss = loss_func(batch_u, model, h, M)   
        optimizer.zero_grad()      
        train_loss.backward()    
        optimizer.step()    
        epoch_avg_loss = train_loss.item()         
        loss_history.append(epoch_avg_loss) 
    
    return loss_history

其中loss_func即为上述定义的均方误差损失函数MSE：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

# 损失函数
def loss_func(u, model, h, M):
	# 参数说明：模型输入数据u,待训练模型model,时间步长h,步数M

    # 根据选定的方法确定系数α、β（示例：Adams-Moulton法，M=4）
    alpha = [1, -1, 0, 0, 0]
    beta = [251/720, 646/720, -264/720, 106/720, -19/720]
        
    loss = 0.0
    for n in range(M, len(u)):
        residual = 0.0
        for m in range(M+1):
            u_nm = u[n - m].unsqueeze(0)
            f_nm = model(u_nm)
            residual += alpha[m] * u[n - m] + h * beta[m] * f_nm.squeeze()
        loss += torch.norm(residual) ** 2
    return loss / (len(u) - M + 1)

基于多步神经网络改进的常微分方程反演算法

改进的多步神经网络算法

多步神经网络做出了比较理想的假设，即给出了所有的数值解数据是无噪声的，以及神经网络可以表示对常微分方程产生零残差的反演。这种假设是理想化的，我们可以在实际情况下使用改进的方法尝试发现未知常微分方程。

对于多步神经网络，除线性多步法残差之外，如果神经网络能够满足一些通用近似性质，尤其是在多步神经网络反演法效果不佳的非自治常微分方程的反演问题中，可以扩展精确和完整数据的收敛结果。考虑到神经网络对函数近似的强大能力（见通用近似定理），这意味着至少在合适的光滑的常微分方程中存在良好的泛化误差。

假设数据集为给定在若干时刻t_1、t_2、…、t_N的方程的数值解u的值，把每一时刻t_n对应的解u(t_n)记为u^(n)，而t_n对应的导数数据一般无从得知。当时间间隔h较小时，可用二阶中心差分法近似求出导数数据，记对于每一时刻t_n求出的导数数据为d^(n)，则：
$$
d^{(n)}=\begin{cases}(-3\boldsymbol{u}^{(n)}+4\boldsymbol{u}^{(n+1)}-\boldsymbol{u}^{(n+2)})/2h,\quad n=1,\(\boldsymbol{u}^{(n+1)}-\boldsymbol{u}^{(n-1)})/2h,\quad1<n<N,\(\boldsymbol{u}^{(n-2)}-4\boldsymbol{u}^{(n-1)}+3\boldsymbol{u}^{(n)})/2h,\quad n=N.&\end{cases}
$$
显然二阶中心差分法求导数d^(n)的误差为o(h^2)。

对于常微分方程的反演，我们改进了多步神经网络的损失函数，其中既包含解数据，又包含导数数据：
$$
L=l_1(u,d)+l_2(u,d,f_{NN})
$$
其中：
$$
l_{1}=\frac{1}{N-M+1}\sum_{n=M}^{N}\left|\sum_{m=0}^{M}[\alpha_{m}u_{n-m}+h\beta_{m}f_{NN}(u_{n-m})]\right|^{2}
$$

$$
l_{2}=\frac{1}{N-M+1}\sum_{n=M}^{N}\left|f_{NN}(u_{n})-d_{n}\right|^{2}
$$

$$
n=M,\ldots,N
$$

式中向量范数采用二范数；u为数值解的数据，d为二阶中心差分法得到的导数数据，则有：

• l_1：线性多步法残差，可以表示常微分方程的近似程度
• l_2：均方误差函数，可以表示网络结构的准确度

把u^(n)作为输入，d^(n)作为期望输出进行训练，得到的f_NN即为对函数f的近似。训练流程与改进后的损失函数Python实现如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

# 模型训练 
def train(u_tensor, model, loss_func, h, M): 
    # 参数说明：训练集u_tensor,待训练模型model,损失函数loss_func,时间步长h,步数M,最大训练次数EPOCH
    optimizer = torch.optim.Adam(model.parameters(), lr=LR)
    loss_history = []
    
    for epoch in range(EPOCHS):
        model.train()
        batch_u = u_tensor.to(device)
        batch_d = d_tensor.to(device)
        train_loss = loss_func(batch_u, batch_d, model, h, M)
        optimizer.zero_grad()
        train_loss.backward()
        optimizer.step()
        epoch_avg_loss = train_loss.item()
        loss_history.append(epoch_avg_loss)
        
    return loss_history

可以看到，与改进前相比，在训练过程中增加了对二阶中心差分法得到的导数的相关计算，进一步利用了训练数据中的信息。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

# 损失函数
def loss_func(u, d, model, h, M):
    # 参数说明：模型输入数据u,输入数据对应的真实微分数据d,待训练模型model,时间步长h,步数M

    # 根据选定的方法确定系数α、β（示例：Adams-Moulton法，M=4）
    alpha = [1, -1, 0, 0, 0]
    beta = [251/720, 646/720, -264/720, 106/720, -19/720]
        
    # 多步法残差 l1
    l1 = 0.0
    for n in range(M, len(u)):
        residual = 0.0
        for m in range(M+1):
            u_nm = u[n - m].unsqueeze(0)
            f_nm = model(u_nm)
            residual += alpha[m] * u[n - m] + h * beta[m] * f_nm.squeeze()
        l1 += torch.norm(residual) ** 2
    l1 = l1 / (len(u) - M + 1)
    
    # 导数拟合误差 l2
    pred_d = model(u)
    l2 = torch.mean(torch.norm(pred_d - d, dim=1) ** 2)
    
    return l1 + l2  # 总损失

值得注意的是，针对非自治常微分方程，由于在将其变换为自治常微分方程时进行了升维操作，这意味着真实的导数向量在y=t这一分量上的值必然是1，因此后续在检验模型训练效果时可利用这一点对模型性能进行先决的把握。

基于误差界理论论证改进算法的优越性

在上述损失函数中，l_1是对目标函数的良好近似，而l_2是神经网络中常用的损失函数。考虑到神经网络对函数近似的强大能力，这意味着即使对于噪声数据，改进的方法也应有良好的泛化能力和鲁棒性。在这里，我们先从理论上来论证改进方法在泛化能力与鲁棒性上的优越性，后续的数值实验部分将采用三类非自治方程的反演来验证改进的多步神经网络的效果。

对于二分类问题（可推广至函数近似的回归问题），当假设空间是有限个函数的集合F={f_1,f_2,…,f_d}时，对任意一个函数f∈F,至少以概率1-δ，以下不等式成立：
$$
R(f)\leq\hat{R}(f)+\varepsilon(d,N,\delta)
$$
其中：
$$
\begin{gathered}R(f)=E[L(Y,f(X))]\\hat{R}(f)=\frac{1}{N}\sum_{1}^{N}L(y_{i},f(x_{i}))\\varepsilon(d,N,\delta)=\sqrt{\frac{1}{2N}(\log d+\log\frac{1}{\delta})}\end{gathered}
$$
式中，R(f)为泛化误差（测试集上的测试风险），\hat{R}(f)为训练集上的经验风险，\hat{R}(f) + ε(d,N,δ)即为泛化误差界。观察上式可知，泛化误差界与样本数N成正比，与假设空间包含的函数数量d成反比。因此，当样本数N越大，泛化误差界越小，当假设空间F包含的函数越多，泛化误差界越大。

根据上述定理，针对改进后的神经网络，有如下不等式成立：
$$
\left|f_{NN}(u^{(n)})-\hat{f}(u^{(n)})\right|\leq\left|f_{NN}(u^{(n)})-d^{(n)}\right|+\left|\hat{f}(u^{(n)})-d^{(n)}\right|\leq w+o(h^{2})
$$
该不等式表明，改进后的损失函数可以拆解成两部分误差：

• 第一部分为通过神经网络得到的导数近似值与通过二阶中心差分法求出的导数之间的误差，设该误差限为w；
• 第二部分为导数真实值与通过二阶中心差分法求出的导数之间的误差，根据二阶中心差分法的基本原理，其误差限为**o(h^2)**。

显然改进后的多步神经网络的泛化误差相较改进前有了进一步的约束，尤其显著减小了噪声对模型的影响，具有更强的泛化能力和鲁棒性。

数值实验与分析

为充分验证改进方法的泛化能力与鲁棒性，我们选取了六个实际的物理问题（对应常微分方程）作为测试用例，对应自治/非自治以及三种不同的线性多步法，在每一个测试用例内选取不同的线性多步法步数以及训练数据所含的高斯噪声方差，对改进前后的两种基于多步神经网络的常微分方程反演方法进行测试。下面我们以“受迫振动方程”这一用例详细展示测试流程，其他的测试用例仅作结果展示。

数值实验流程

数据构建

基于选定的测试用例对应的常微分方程，可采用如下步骤生成原始的真实数据集：

1. 在所给区域内均匀选取间隔为h的N个结点t_n,n=1,2,…,N ,使用Python中的Scipy库求出各结点对应的数值解u^(n)；
2. 将数值解u^(n) 代入二阶中心差分法中得到每个x_n 对应的导数值d^(n) ；
3. 噪声数据构建：在各结点对应的数值解u^(n) 上依次加上期望为0，方差为0、0.01、0.05的高斯噪声。

在Python中，用函数get_data()函数实现了数据构建的代码封装：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

def get_data(t0, t_end, h, y0, noise_std=0.0):     
    t_eval = np.arange(t0, t_end + h, h)     
    
    # 数值解    
    sol = solve_ivp(forced_oscillator, [t0, t_end], y0, t_eval=t_eval)    
    x = sol.y.T  # shape: [N, 2]     
    t = sol.t.reshape(-1, 1)  # shape: [N, 1]     
    
    # 添加高斯噪声     
    if noise_std > 0:         
        x += np.random.normal(0, noise_std, size=x.shape)      
        
    # 构造 u = [x1, x2, t]     
    u = np.hstack([x, t])      
    
    # 二阶中心差分求导数 d ≈ du/dt     
    d = np.zeros_like(u)     
    # 边界点使用前向/后向差分     
    d[0] = (-3 * u[0] + 4 * u[1] - u[2]) / (2 * h)     
    d[-1] = (u[-3] - 4 * u[-2] + 3 * u[-1]) / (2 * h)     
    # 内部点使用中心差分     
    d[1:-1] = (u[2:] - u[:-2]) / (2 * h)      
    
    u_tensor = torch.tensor(u, dtype=torch.float32)     
    d_tensor = torch.tensor(d, dtype=torch.float32)      
    
    return u_tensor, d_tensor, t# 模型训练 

以“受迫振动方程”这一测试框架为例，force_oscillator函数中刻画了该场景下的微分方程规律：

1
2
3
4
5

def forced_oscillator(t, y):
    x1, x2 = y
    dx1dt = x2
    dx2dt = -np.cos(t) * x1 - x2 + t / 50
    return [dx1dt, dx2dt]

生成的无噪声数据如下图所示：

神经网络训练

根据上述算法理论部分列出的训练流程框架，适用数据（u，d）对神经网络进行训练，训练后得到的函数f_NN即为对函数f的近似。

关于神经网络的结构与参数，本项目中所使用的神经网络结构均相同（后续不再赘述），均包含4个隐藏层，每层128个神经元，激活函数为tanh，Python中基于pytorch框架对神经网络架构的实现如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

import torch.nn as nn

class ODEApproximator(nn.Module):
    def __init__(self, input_dim=3, hidden_dim=128, output_dim=3):
        super(ODEApproximator, self).__init__()
        self.net = nn.Sequential(
            nn.Linear(input_dim, hidden_dim),
            nn.Tanh(),
            nn.Linear(hidden_dim, hidden_dim),
            nn.Tanh(),
            nn.Linear(hidden_dim, hidden_dim),
            nn.Tanh(),
            nn.Linear(hidden_dim, hidden_dim),
            nn.Tanh(),
            nn.Linear(hidden_dim, output_dim)
        )

    def forward(self, x):
        return self.net(x)

训练过程中实时记录了每一个epoch的损失函数值，下图分别绘制了改进前后的多步神经网络在无噪声数据集上训练过程中损失函数值的变化情况（总Epoch=100），可以看到损失函数均呈收敛态势，且改进后的多步神经网络收敛更快：

效果检验

为定量地比较改进前后的多步神经网络的泛化性能与鲁棒性，采取如下的实验步骤进行训练后模型的效果检验：

1. 用训练得到的函数f_NN解方程得到数值解u_NN^(n) ,观察近似效果，计算均方误差mse_x 和平均绝对误差mae_x并与多步神经网络作比较；
2. 用数值解u^(n)代入训练得到的函数f_NN得到的导数的近似值d_NN^(n)，将u^(n)代入给定的常微分方程中的 f 得到导数的真实值d^(n),计算均方误差mse_f和平均绝对误差mae_f并与多步神经网络作比较。

其中用于比较性能的指标均方误差mse和平均绝对误差mae的定义式如下所示：
$$
mse_{x}=\frac{1}{ND}\sum_{n=1}^{N}\left|u_{NN}^{(n)}-u^{(n)}\right|_{2}^{2}
$$

$$
mse_{f}=\frac{1}{ND}\sum_{n=1}^{N}\left|d_{NN}^{(n)}-d^{(n)}\right|_{2}^{2}
$$

$$
mae_{x}=\frac{1}{ND}\sum_{n=1}^{N}\left|u_{NN}^{(n)}-u^{(n)}\right|_{1}
$$

$$
mae_{f}=\frac{1}{ND}\sum_{n=1}^{N}\left|d_{NN}^{(n)}-d^{(n)}\right|_{1}
$$

Python中将利用训练好的模型进行预测的过程封装成函数predict_u，指标计算直接调用scikit-learn函数库中的相关函数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

def rhs(t, u_np):      
    """用于 solve_ivp：输入 u = [x1, x2, t] 输出 du/dt"""      
    u_tensor = torch.tensor(u_np, dtype=torch.float32).unsqueeze(0).to(device)      
    with torch.no_grad():          
        du_dt = model(u_tensor).squeeze().cpu().numpy()      
    return du_dt    

def predict_u(u_tensor):      
    # 初始条件与时间点一致     
    u0 = u_tensor[0].cpu().numpy()  # [x1, x2, t]    
    t_span = (0, 100)  
    t_eval = np.arange(0, 100 + 0.05, 0.05)    
    
    sol_nn = solve_ivp(rhs, t_span, u0, t_eval=t_eval)   
    return sol_nn.y.T  # shape: [N, 3]

在“受迫振动方程”这一测试案例中，基于训练轮数Epoch=50的模型，运行以上检验流程，得到如下检验指标数据：

可以看到，随着噪声强度增加，两种方法的误差均增大，但改进方法对数值解的预测误差在各种情况下均小于多步神经网络。

可以看到，当数据无噪声时，本文提出的改进方法对于受迫振动方程的近似效果强于多步神经网络;随着噪声强度的增加，多步神经网络的误差大幅增加，而改进的方法的误差一直保持较小且始终强于多步神经网络。

除此之外，针对在无噪声数据集上训练轮数分别为Epoch=100和Epoch=500时的模型，将其作为微分方程中的f用于预测数据的生成，效果分别如下图所示：

可以看到，随着训练轮数的增加，模型对于实际微分方程的拟合效果越来越好，生成的预测数据也越来越接近用于训练的真实数据；事实上，从理论上讲，大概在训练轮数Epoch达到10^4数量级时才会有较为精准的拟合效果，本项目由于用于训练的计算资源有限，仅进行了最大训练轮数为1000的测试。在后续的结果分析部分，将直接采用高训练轮数下的精确结果进行展示。

数值实验结果分析

非自治常微分方程

受迫振动方程——Adams-Moulton法

受迫振动方程(Forced Oscillator)是描述一个振动系统在外界力的作用下进行振动的方程。这个方程可以用来研究各种物理系统中的振动现象，例如弹簧振子、摆锤、电路振荡等。受迫振动方程的核心思想是将外界力引入方程中，以描述振动系统在外界激励下的行为。受迫振动方程的研究对于科学、工程和技术领域具有重要意义。它不仅有助于我们深入理解振动现象的本质和规律，还为我们设计和优化振动控制系统、减振装置等提供了理论基础。同时，它在电子设备、结构工程、交通运输等方面也有广泛的应用。

实验使用的受迫振动方程如下：
$$
\begin{aligned}&\frac{dx_1}{dt}=x_2,\&\frac{dx_2}{dt}=-cosy\cdot x_1-x_2+\frac{y}{50},\&\frac{dy}{dt}=1\end{aligned}
$$
使用[0,1,0]T作为初值，生成从t=0到t=100，间隔h=0.05的数据作为训练集，最终得到的多步神经网络对常微分方程的反演效果如下图所示，其中红色曲线是原方程数值解的曲线，蓝色点是用无噪声数据训练的函数f_NN代替f得到的数值解:

线性标量方程——Adams-Bashforth法

线性标量方程(Linear Scalar Equation)是一种描述某个未知函数与其导数之间的关系的方程，其中未知函数的导数的最高次数为1。线性标量方程是微分方程中的一类重要问题，它在数学物理等领域中有广泛应用。线性标量方程的研究对于数学、物理学和工程学等领域具有重要意义。它可以用来描述动力学系统的行为、传热和传质过程、量子力学中的波函数演化等。通过分析和求解线性标量方程，可以深入理解系统的特征、稳定性、渐近行为等，为问题的模拟、控制和优化提供理论基础。在现代科学和工程中，线性标量方程的研究得到了进一步推广和延伸。例如，非线性标量方程、偏微分方程等是线性标量方程的扩展和推广，它们更好地描述了一些复杂的物理现象和现实问题。

实验使用的线性标量方程如下：
$$
\frac{dx_1}{dt}=-x_1(sin4y+1)+cos\frac{y^2}{1000},
$$

$$
\frac{dy}{dt}=1
$$

使用[2,0]T作为初值，生成从t=0到t=20，间隔h=0.05的数据作为训练集，最终得到的多步神经网络对常微分方程的反演效果如下图所示，其中红色曲线是原方程数值解的曲线，蓝色点是用无噪声数据训练的函数f_NN代替f得到的数值解:

食饵-捕食者模型方程——后向微分公式法

食饵-捕食者模型(Predator-prey Model)是生态学中研究食物链和生物群落动力学的重要模型之一，描述了食物链中食饵(被捕食者)和捕食者之间的相互作用关系。这种模型的研究对于我们理解生态系统的稳定性、物种相互作用以及生态系统中物种数量的变化具有重要意义。通过食饵-捕食者模型，我们可以研究食饵与捕食者之间的数量关系以及它们之间的相互作用。一般来说，食饵的数量被认为是捕食者的食物来源，并且捕食者的数量受到食饵的供给和捕食行为的影响。该模型通常描述了食饵数量随时间的变化，以及捕食者数量随时间的变化。食饵-捕食者模型的研究对于生态学、环境保护和资源管理等领域具有重要意义。通过该模型，我们可以深入了解生态系统中物种数量的变化规律，从而预测物种灭绝和生态系统崩溃的风险，以及寻找保护和管理生物多样性的方法。

实验使用的食饵-捕食者模型方程如下：
$$
\begin{aligned}&\frac{dx_1}{dt}=x_1-x_1\cdot x_2+sin\frac{y}{2}+cosy+2,\&\frac{dx_2}{dt}=x_1\cdot x_2-x_2,\&\frac{dy}{dt}=1\end{aligned}
$$
使用[3,3,0]T作为初值，生成从t=0到t=50，间隔h=0.05的数据作为训练集，最终得到的多步神经网络对常微分方程的反演效果如下图所示，其中红色曲线是原方程数值解的曲线，蓝色点是用无噪声数据训练的函数f_NN代替f得到的数值解:

自治常微分方程

线性常微分方程——Adams-Moulton法

线性常微分方程(Linear ODEs)是微分方程中的一类重要问题，它描述了未知函数及其导数之间的线性关系。线性常微分方程的研究对于数学、物理学和工程学等领域具有重要意义。它们出现在物理学中许多基础问题的数学描述中，例如弹簧振动、电路分析、传热过程等。在工程学中，线性常微分方程也广泛应用于控制系统的建模和分析。随着科学技术的发展，研究者提出了各种各样的数值方法和近似方法，以求解复杂的线性常微分方程。这些方法包括欧拉方法、龙格-库塔方法、有限差分方法、变分法等，它们为实际问题的求解提供了有效的数值工具。随着对非线性动力学系统的研究和认识的深入，研究者们开始将非线性常微分方程引入到线性常微分方程中，以描述更为复杂的现象和现实问题。

实验使用的线性常微分方程如下：
$$
\frac{dx_1}{dt}=x_1-4x_2,
$$

$$
\frac{dx_2}{dt}=4x_1-7x_2
$$

使用[0,-1]T作为初值，生成从t=0到t=10，间隔h=0.01的数据作为训练集，最终得到的多步神经网络对常微分方程的反演效果如下图所示，其中红色曲线是原方程数值解的曲线，蓝色点是用无噪声数据训练的函数f_NN代替f得到的数值解:

阻尼三次振子方程——Adams-Bashforth法

阻尼三次振子(Damped Cubic Oscillator)是一种物理系统，它描述了受到阻尼力和弹力作用的振动系统。阻尼指的是系统中存在能量损耗的因素，如摩擦力，它会导致振动的逐渐减弱；三次振子表示系统的势能函数是一个三次函数，它具有非线性的特性。阻尼三次振子最早出现在振动力学和非线性动力学的研究中。对于振动系统的研究是物理学、工程学和应用数学等领域的重要课题之一。阻尼三次振子的研究对于理解和预测复杂动力学系统的行为具有重要意义。通过分析阻尼三次振子的运动规律和稳定性，我们可以深入了解非线性振动系统的动力学特性，例如振动的周期性、混沌行为等。这些研究也为控制系统、电子电路、力学系统等领域的工程应用提供了参考和启示。

实验使用的阻尼三次振子方程如下：
$$
\frac{dx_1}{dt}=-0.1x_1^3+2x_2^3,
$$

$$
\frac{dx_2}{dt}=-2x_1^3-0.1x_2^3
$$

使用[1,1]T作为初值，生成从t=0到t=25，间隔h=0.01的数据作为训练集，最终得到的多步神经网络对常微分方程的反演效果如下图所示，其中红色曲线是原方程数值解的曲线，蓝色点是用无噪声数据训练的函数f_NN代替f得到的数值解:

阻尼简谐摆方程——后向微分公式法

阻尼简谐摆(Damped Pendulum)是一个经典力学中的物理系统，它由一个具有质量的物体通过一根轻质绳或杆与一个固定支点相连接组成。阻尼简谐摆在受到重力作用下，沿着一条弧线进行周期性振动。阻尼简谐摆的研究对于理解和预测振动系统的行为具有重要意义。通过分析阻尼简谐摆的运动规律和稳定性，我们可以深入了解振动系统在存在阻尼时的响应特性，例如振动的振幅、频率等。这些研究不仅在理论物理学中具有重要意义，也为工程学中的控制系统、机械振动和结构动力学等领域的应用提供了基础。

实验使用的阻尼简谐摆方程如下：
$$
\begin{aligned}&\frac{dx_1}{dt}=x_2,\&\frac{dx_2}{dt}=-\alpha x_2-\beta\sin x_1,\&\alpha=0.2,\beta=8.91\end{aligned}
$$
使用[-1.193,-3.876]T作为初值，生成从t=0到t=20，间隔h=0.01的数据作为训练集，最终得到的多步神经网络对常微分方程的反演效果如下图所示，其中红色曲线是原方程数值解的曲线，蓝色点是用无噪声数据训练的函数f_NN代替f得到的数值解:

研究结论与创新

研究结论

本项目针对常微分方程的反演问题提出了一种基于多步神经网络改进的常微分方程反演算法，旨在从数据中提取内含的常微分方程。通过将导数数据融入训练数据，我们改进了多步神经网络的数据集，并优化了神经网络的损失函数，以提高模型的准确性和泛化能力。经过训练，改进的神经网络成功拟合原方程中的函数，实现了对非自治常微分方程的反演。我们还推导并讨论了改进算法的误差界。我们进一步将改进的算法应用于自治常微分方程反演，通过训练网络拟合原方程函数，拓展了算法在自治常微分方程反演中的适用性。

本项目通过实验以非自治常微分方程中的受迫振动方程、线性标量方程和食饵-捕食者模型方程，以及自治常微分方程中的线性常微分方程、阻尼三次振子方程和阻尼简谐摆方程为例，展示了算法的有效性。在实验中，我们添加了不同强度的噪声，并利用改进算法对常微分方程进行反演。实验结果表明，改进方法在大多数情况下优于多步神经网络，在处理有噪声数据时表现更为出色。这些结果充分验证了本文提出的方法在反演问题中的可行性和优越性。

研究方法创新性

传统方法通常利用多步神经网络来处理常微分方程的反演问题，主要集中在处理自治常微分方程，对于非自治常微分方程的反演效果不尽理想，特别是在使用带有噪声数据进行训练时表现不佳。本研究首次将研究重点聚焦于非自治常微分方程的反演，通过将非自治方程转化为自治方程的形式，提出了一种新的研究思路。针对带有噪声数据的方程反演问题，采用了改进损失函数的方法以提升算法的鲁棒性。此外，将这种新方法扩展应用于自治常微分方程的反演问题，结果显示其反演效果明显优于多步神经网络。这一新思路对解决常微分方程反演问题带来了一种新的方法。

参考文献

[1] 王开荣,杨大地编著.应用数值分析[M].高等教育出版社,2010.

[2] 付长铠.常微分方程反演的机器学习方法研究[D].长春工业大学,2024.

[3] Raissi M, Perdikaris P, Karniadakis G E. Multistep neural networks for data-driven discovery of nonlinear dynamical systems[J]. arXiv preprint arXiv:1801.01236, 2018.

[4] 李航. 统计学习方法[M]. 第二版. 北京：清华大学出版社, 2019.

[5] 陈新海, 刘杰, 万仟, 等. 一种改进的基于深度神经网络的偏微分方程求解方法[J]. 计算机工程与科学, 2022, 44(11): 1932-1940.

研究背景及其意义

图像放大与旋转是数字图像处理中最基础的几何变换操作，其核心在于如何通过插值算法重建原始图像中不存在的像素信息。当对图像进行放大操作时，输出图像的像素网格会超出原始图像的采样范围，需要通过插值来填补这些新增像素点的颜色值；而在旋转操作中，即使保持图像尺寸不变，原始像素的整数坐标经过旋转变换后也会落在新图像的非整数位置，同样需要通过插值来重新确定每个输出像素的颜色值。

图像插值是利用原图像中的颜色值通过一定的方法计算出待插入像素点的颜色值的过程。对图像进行插值一般有两个步骤：首先定义一个图像插值公式，然后利用该插值公式计算待插入点的颜色值。常见的图像插值算法有双线性法、最近邻法、非均匀法、双三次卷积插值法、双立方法、Lagrange法、 样条插值法、 克里金（Krijing） 插值法等。这些插值方法通常定义一个插值数据点的隐式函数，再提取该函数的等值面作为图像插值方法，常用的插值核包括线性插值核、样条插值核等。

• 最近邻插值作为最简单的算法，直接将距离待插值点最近的已知像素值作为结果，虽然计算效率极高（时间复杂度O(1)），但会产生明显的块状伪影（“马赛克”）和锯齿形边缘；
• 双线性插值通过考虑2×2邻域内四个像素的加权平均，在计算成本（O(n)）和视觉效果之间取得平衡，但仍会导致高频信息丢失和边缘模糊；
• 更高阶的双三次插值（使用4×4邻域）和样条插值虽然能提供更平滑的结果，但计算复杂度显著增加（O(n²)），且可能引入不必要的振铃效应。

现有算法的根本局限在于采用统一的插值核函数处理整幅图像，忽视了图像不同区域的特征差异。例如，在平坦区域使用复杂插值会造成计算资源浪费，而在纹理丰富区域使用简单插值又会导致细节损失。基于此，我们希望通过改良的四平面插值算法对图像的放大与旋转效果进行优化，根据图像局部特征自适应地选择不同的插值策略，以规避用同一个插值公式对所有像素进行插值存在的不足。

常用图像插值算法

课本在6.5节中提到，在插值节点数量较多时，为避免Runge振荡现象的发生，并不提倡用高次多项式进行插值，而宁可用低次多项式作分段插值。在图像处理这一特定的应用场景中，需要处理的图像尺寸规模往往较大，且同一行（列）的所有像素颜色值显然并不具有可以用一个多项式函数显式表达的规律，但相邻的像素点颜色值之间又存在一定的关联性，因此分段插值仅考虑局部特征的特性在这里能够良好地契合所需性能。根据对于待插入像素点周围已有的像素点信息的利用情况，这里列举了几种常见的图像插值算法：

• 最近邻法：仅利用待插值像素点转换至原图像坐标后距离其最近的一个像素点的颜色值，将其直接作为待插值像素点的颜色值
• 双线性法：利用待插值像素点转换至原图像坐标后距离其最近的四个像素点的颜色值，加权平均后作为待插值像素点的颜色值
• 双立方法：利用待插值像素点转换至原图像坐标后距离其最近的十六个像素点的颜色值，加权平均后作为待插值像素点的颜色值

最近邻法

如上图所示，在一维最近邻插值中，坐标轴上各点 xi-1，xi，xi+1 … 两两对半等分间隔 (红色虚线划分)，从而非边界的各坐标点都有一个等宽的邻域，并根据每个坐标点的值构成一个类似分段函数的函数约束，从而使各插值坐标点的值等同于所在邻域原坐标点的值。例如，插值点 x 坐落于 坐标点 xi 的邻域，那么其值 f(x) 就等于 f(xi)。

在二维的图像插值场景中，可以对上述一维最近邻插值进行推广，如下图所示：

可以看到，(x0, y0)、(x0, y1)、(x1, y0)、(x1, y1) 都是原图像上的坐标点，颜色值分别对应为 Q11、Q12、Q21、Q22。而颜色值未知的插值点 (x, y)（需转换至原图像坐标），根据最近邻插值方法的约束，其与坐标点 (x0, y0) 位置最接近 (即位于 (x0, y0) 的邻域内)，故插值点 (x, y) 的颜色值 P = Q11。

总而言之，最近邻法的基本思想即：将待插入点的坐标进行四舍五入，再以该行列坐标都是整数点的颜色值（灰度值）替代待插入点(x, y)处的颜色值。事实上，这也正是机器学习中KNN（K-Nearest Neighbor）算法在K=1时的情形。

基于以上算法思想，编写python函数代码实现图像放缩与旋转过程中的最近邻法插值：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

# 最近邻法插值实现图像放缩
def nearest_neighbor_interpolation(image, scale_factor):
    h, w, channel = image.shape
    new_h, new_w = int(h * scale_factor), int(w * scale_factor)
    resized_image = np.zeros((new_h, new_w, int(channel)), dtype=image.dtype)
    
    for i in range(new_h):
        for j in range(new_w):
            src_i = int(round((i + 1) / scale_factor, 0))
            src_j = int(round((j + 1) / scale_factor, 0))
            resized_image[i, j] = image[src_i - 1, src_j - 1]
    
    return resized_image

# 最近邻法插值实现图像旋转
def nearest_neighbor_rotation(image, angle):
    h, w, channel = image.shape
    angle_rad = math.radians(angle)
    
    # 计算旋转后的图像尺寸
    cos_theta = abs(math.cos(angle_rad))
    sin_theta = abs(math.sin(angle_rad))
    new_w = int(h * sin_theta + w * cos_theta)
    new_h = int(h * cos_theta + w * sin_theta)
    
    # 旋转中心
    cx, cy = w / 2, h / 2
    new_cx, new_cy = new_w / 2, new_h / 2
    
    rotated_image = np.zeros((new_h, new_w, channel), dtype=image.dtype)
    
    for i in range(new_h):
        for j in range(new_w):
            # 将新图像坐标转换回原图像坐标
            x = (j - new_cx) * math.cos(angle_rad) + (i - new_cy) * math.sin(angle_rad) + cx
            y = -(j - new_cx) * math.sin(angle_rad) + (i - new_cy) * math.cos(angle_rad) + cy
            
            # 最近邻插值
            if 0 <= x < w and 0 <= y < h:
                src_x = int(round(x))
                src_y = int(round(y))
                rotated_image[i, j] = image[src_y - 1, src_x - 1]
    
    return rotated_image

双线性法

如上图所示，在一维的线性插值中，坐标轴上各点 xi-1，xi，xi+1 … 的值“两两直接相连”为线段，从而构成了一条连续的约束函数。而插值坐标点例如 x，根据约束函数其值应为 f(x)。因为每两个坐标点之间的约束函数曲线是一次线性的线段，对插值结果而言是“线性” 的，所以该方法称为线性插值。基于线性函数的特性，可以便捷地求取原图像上的两个像素点间任一待插值点的颜色值：

可以看到，图中 x0 和 x1 都是原有的坐标点，颜色值分别对应为 y0 和 y1，此时根据线性插值法约束，在 (x0, y0) 和 (x1, y1) 构成的一次函数上，颜色值未知的插值点 x的颜色值 y 即为：
$$
y=y_0+(x-x_0)\frac{y_1-y_0}{x_1-x_0}=y_0+\frac{(x-x_0)y_1-(x-x_0)y_0}{x_1-x_0}
$$
实际上，即便 x 不在 x0 与 x1 之间，该公式也成立（此时为线性外插），但图像处理中不需涉及此情形。

从一维的线性插值出发，很容易拓展到二维图像的双线性插值，通过三次一阶线性插值（本质为加权求和）获得最终结果，下图便展示了该过程的定性斜视与定量俯视示意图：

其中，(x0, y0)、(x0, y1)、(x1, y0)、(x1, y1) 均为原图像上的像素坐标点，颜色值分别对应为 f(x0, y0)、f(x0, y1)、f(x1, y0)、f(x1, y1)。而颜色值未知的插值点 (x, y)，根据双线性插值法的约束，可以先由像素坐标点 (x0, y0) 和 (x0, y1) 在 y 轴向作一维线性插值得到 f(x0, y)、由像素坐标点 (x1, y0) 和 (x1, y1) 在 y 轴向作一维线性插值得到 f(x1, y)，然后再由 (x0, y) 和 (x1, y) 在 x 轴向作一维线性插值得到插值点 (x, y) 的灰度值 f(x, y)。

事实上，一维线性插值先作 x 轴向再作 y 轴向，得到的结果完全相同，仅为顺序先后的区别。这里不妨先由像素坐标点 (x0, y0) 和 (x1, y0) 在 x 轴向作一维线性插值得到 f(x, y0)、由像素坐标点 (x0, y1) 和 (x1, y1) 在 x 轴向作一维线性插值得到 f(x, y1)：
$$
f(x,y_0)=\frac{x_1-x}{x_1-x_0}f(x_0,y_0)+\frac{x-x_0}{x_1-x_0}f(x_1,y_0)
$$

$$
f(x,y_1)=\frac{x_1-x}{x_1-x_0}f(x_0,y_1)+\frac{x-x_0}{x_1-x_0}f(x_1,y_1)
$$

然后再由 (x, y0) 和 (x, y1) 在 y 轴向作一维线性插值得到插值点 (x, y) 的灰度值 f(x, y)：
$$
f(x,y)=\frac{y_1-y}{y_1-y_0}f(x,y_0)+\frac{y-y_0}{y_1-y_0}f(x,y_1)
$$
合并上述式子，得到最终的双线性插值结果：
$$
f(x,y)=\frac{(y_1-y)(x_1-x)}{(y_1-y_0)(x_1-x_0)}f(x_0,y_0)+\frac{(y_1-y)(x-x_0)}{(y_1-y_0)(x_1-x_0)}f(x_1,y_0)+\frac{(y-y_0)(x_1-x)}{(y_1-y_0)(x_1-x_0)}f(x_0,y_1)+\frac{(y-y_0)(x-x_0)}{(y_1-y_0)(x_1-x_0)}
$$
值得注意的是，在实际的图像插值处理过程中，为尽量保证插值效果的准确性，往往仅采用距离待插值点（转换至原图像坐标）最近的四个点，即:（[]符号表示待插值点转换至原图像坐标后向下取整）
$$
x_0=[x]，y_0=[y]
$$

$$
x_1=x_0+1，y_1=y_0+1
$$

从加权求和的角度理解，可以进一步地将双线性插值结果改写为如下形式：
$$
p=x-[x], q=y-[y]
$$

$$
\begin{array}{rcl}f(x,y)=(1-q){(1-p)f([x][y])+pf([x]+1,[y])}+q{(1-p)f([x],[y]+1)+pf([x]+1,[y]+1)}\end{array}
$$

基于以上算法思想，编写python函数代码实现图像放缩与旋转过程中的双线性法插值：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72

# 双线性法插值实现图像放缩
def bilinear_interpolation(image, scale_factor):
    h, w, channel = image.shape
    new_h, new_w = int(h * scale_factor), int(w * scale_factor)
    resized_image = np.zeros((new_h, new_w, int(channel)), dtype=image.dtype)
    
    for i in range(new_h):
        for j in range(new_w):
            x = (j + 1) / scale_factor
            y = (i + 1) / scale_factor
            x1 = int(x)
            y1 = int(y)
            x2 = x1 + 1
            y2 = y1 + 1
            p = x - x1
            q = y - y1

            # 边界问题处理
            if x2 == w + 1:
                x2 = x1
            if y2 == h + 1:
                y2 = y1
                
            resized_image[i ,j] = (1 - q) * ((1 - p) * image[y1 - 1, x1 - 1] + p * image[y1 - 1, x2 - 1]) + q * ((1 - p) * image[y2 - 1, x1 - 1] + p * image[y2 - 1, x2 - 1])
    
    return resized_image

# 双线性法插值实现图像旋转
def bilinear_rotation(image, angle):
    h, w, channel = image.shape
    angle_rad = math.radians(angle)
    
    # 计算旋转后的图像尺寸
    cos_theta = abs(math.cos(angle_rad))
    sin_theta = abs(math.sin(angle_rad))
    new_w = int(h * sin_theta + w * cos_theta)
    new_h = int(h * cos_theta + w * sin_theta)
    
    # 旋转中心
    cx, cy = w / 2, h / 2
    new_cx, new_cy = new_w / 2, new_h / 2
    
    rotated_image = np.zeros((new_h, new_w, channel), dtype=image.dtype)
    
    for i in range(new_h):
        for j in range(new_w):
            # 将新图像坐标转换回原图像坐标
            x = (j - new_cx) * math.cos(angle_rad) + (i - new_cy) * math.sin(angle_rad) + cx
            y = -(j - new_cx) * math.sin(angle_rad) + (i - new_cy) * math.cos(angle_rad) + cy
            
            # 双线性插值
            if 0 <= x < w-1 and 0 <= y < h-1:
                x1, y1 = int(x), int(y)
                x2, y2 = min(x1 + 1, w - 1), min(y1 + 1, h - 1)
                
                # 计算权重
                a = x - x1
                b = y - y1
                
                # 边界处理
                if x2 >= w:
                    x2 = x1
                if y2 >= h:
                    y2 = y1
                
                # 插值计算
                rotated_image[i, j] = (1 - a) * (1 - b) * image[y1, x1] + \
                                     a * (1 - b) * image[y1, x2] + \
                                     (1 - a) * b * image[y2, x1] + \
                                     a * b * image[y2, x2]
    
    return rotated_image    

双立方法

双立方法插值又称立方卷积插值/双三次插值，这也是数值分析中最常用的二维插值方法。在这种方法中，插值点 (x, y) 的像素颜色值 f(x, y) 通过矩形网格中最近的十六个采样点的加权平均得到，而各采样点的权重由该点到待求插值点的距离确定，此距离包括水平和竖直两个方向上的距离。相比之下，双线性插值仅由周围的四个采样点加权得到。

如上图所示，设（转换至原图像中）待求插值点坐标为 (i+u, j+v)【i、j为整数部分，u、v为小数部分】，已知其周围的 16 个像素坐标点 (网格) 的颜色值，还需要计算 16 个点各自的权重。以像素坐标点 (i, j) 为例，因为该点在 y 轴和 x 轴方向上与待求插值点 (i+u, j+v) 的距离分别为 u 和 v，所以其权重为 w(u) × w(v)，其中 w(·) 是插值权重核 (可以理解为定义的权重函数)。同理可得其余 15 个像素坐标点各自的权重。那么，待求插值点 (i+u, j+v) 的颜色值 f(i+u, j+v) 将通过如下计算得到：
$$
f(i+u,j+v)=A\times B\times C
$$
其中各项由向量或矩阵表示为：
$$
\mathrm{A}=[w(1+u)w(u)w(1-u)w(2-u)]
$$

$$
\mathrm{B}=\begin{bmatrix}f(i-1,j-1)&f(i-1,j+0)&f(i-1,j+1)&f(i-1,j+2)\f(i+0,j-1)&f(i+0,j+0)&f(i+0,j+1)&f(i+0,j+2)\f(i+1,j-1)&f(i+1,j+0)&f(i+1,j+1)&f(i+1,j+2)\f(i+2,j-1)&f(i+2,j+0)&f(i+2,j+1)&f(i+2,j+2)\end{bmatrix}
$$

$$
\mathbb{C}=[w(1+v)w(v)w(1-v)w(2-v)]^T
$$

插值权重核 w(·) 为：
$$
w(x)=\begin{cases}1-2|x|^2+|x|^3&,|x|<1\4-8|x|+5|x|^2-|x|^3&,1\leq|x|<2\0&,|x|\geq2&\end{cases}
$$
插值权重核 w(·) 的函数图像：

为方便后续算法实现，将以上加权求和过程各步骤展开，合并后化简得到待插入点的颜色值计算公式：
$$
f(i+u,j+v)=\sum_{m=0}^{3}\sum_{n=0}^{3}a_{mn}u^{m}v^{n}
$$
其中多项式的系数a_{mn}计算公式如下：(式中p {qr}与上述矩阵B中元素一一对应，如p 00=f(i-1,j-1))
$$
\begin{aligned}
&a{00}=p{11}\&a_{01}=-\frac{1}{2}p_{10}+\frac{1}{2}p_{12}\&a_{02}=p_{10}-\frac{5}{2}p_{11}+2p_{12}-\frac{1}{2}p_{13}\&a_{03}=-\frac{1}{2}p_{10}+\frac{3}{2}p_{11}-\frac{3}{2}p_{12}+\frac{1}{2}p_{13}\&a_{10}=-\frac{1}{2}p_{01}+\frac{1}{2}p_{21}\&a_{11}=\frac{1}{4}p_{00}-\frac{1}{4}p_{02}-\frac{1}{4}p_{20}+\frac{1}{4}p_{22}\&a_{12}=-\frac{1}{2}p_{00}+\frac{1}{4}p_{01}-p_{02}+\frac{1}{4}p_{03}+\frac{1}{2}p_{20}-\frac{5}{4}p_{21}+p_{22}-\frac{1}{4}p_{23}\&a_{13}=\frac{1}{4}p_{00}-\frac{3}{4}p_{01}+\frac{3}{4}p_{02}-\frac{1}{4}p_{03}-\frac{1}{4}p_{20}+\frac{3}{4}p_{21}-\frac{3}{4}p_{22}+\frac{1}{4}p_{23}\
&a_{20}=p_{01}-\frac{5}{2}p_{11}+2p_{21}-\frac{1}{2}p_{31}\
&a_{21}=-\frac{1}{2}p_{00}+\frac{1}{2}p_{02}+\frac{5}{4}p_{10}-\frac{5}{4}p_{12}-p_{20}+p_{22}+\frac{1}{4}p_{30}-\frac{1}{4}p_{32}\&a_{22}=p_{00}-\frac{5}{2}p_{01}+2p_{02}-\frac{1}{2}p_{03}-\frac{5}{2}p_{10}+\frac{25}{4}p_{11}-5p_{12}+\frac{5}{4}p_{13}+2p_{20}-5p_{21}+4p_{22}-p_{23}-\frac{1}{2}p_{30}+\frac{5}{4}p_{31}-p_{32}+\frac{1}{4}p_{33}\
&a_{23}=-\frac{1}{2}p_{00}+\frac{3}{2}p_{01}-\frac{3}{2}p_{02}+\frac{1}{2}p_{03}+\frac{5}{4}p_{10}-\frac{15}{4}p_{11}+\frac{15}{4}p_{12}-\frac{5}{4}p_{13}-p_{20}+3p_{21}-3p_{22}+p_{23}+\frac{1}{4}p_{30}-\frac{3}{4}p_{31}+\frac{3}{4}p_{32}-\frac{1}{4}p_{33}\
&a_{30}=-\frac{1}{2}p_{01}+\frac{3}{2}p_{11}-\frac{3}{2}p_{21}+\frac{1}{2}p_{31}\
&a_{31}=\frac{1}{4}p_{00}-\frac{1}{4}p_{02}-\frac{3}{4}p_{10}+\frac{3}{4}p_{12}+\frac{3}{4}p_{20}-\frac{3}{4}p_{22}-\frac{1}{4}p_{30}+\frac{1}{4}p_{32}\&a_{32}=-\frac{1}{2}p_{00}+\frac{5}{4}p_{01}-p_{02}+\frac{1}{4}p_{03}+\frac{3}{2}p_{10}-\frac{15}{4}p_{11}+3p_{12}-\frac{3}{4}p_{13}-\frac{3}{2}p_{20}+\frac{15}{4}p_{21}-3p_{22}+\frac{3}{4}p_{23}+\frac{1}{2}p_{30}-\frac{5}{4}p_{31}+p_{32}-\frac{1}{4}p_{33}\&a_{33}=\frac{1}{4}p_{00}-\frac{3}{4}p_{01}+\frac{3}{4}p_{02}-\frac{1}{4}p_{03}-\frac{3}{4}p_{10}+\frac{9}{4}p_{11}-\frac{9}{4}p_{12}+\frac{3}{4}p_{13}+\frac{3}{4}p_{20}-\frac{9}{4}p_{21}+\frac{9}{4}p_{22}-\frac{3}{4}p_{23}-\frac{1}{4}p_{30}+\frac{3}{4}p_{31}-\frac{3}{4}p_{32}+\frac{1}{4}p_{33}
\end{aligned}
$$
基于以上算法思想，编写python函数代码实现图像放缩与旋转过程中的双立方法插值：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148

# 双立方法插值实现图像放缩
def bicubic_interpolation(image, scale_factor):
    h, w, channel = image.shape
    new_h, new_w = int(h * scale_factor), int(w * scale_factor)
    resized_image = np.zeros((new_h, new_w, channel))
    
    for i in range(new_h):
        for j in range(new_w):
            x = i / scale_factor
            y = j / scale_factor
            
            # 确定16个邻域像素的坐标
            x0 = max(0, int(np.floor(x)) - 1)
            x1 = x0 + 1
            x2 = x0 + 2
            x3 = min(w-1, x0 + 3)
            
            y0 = max(0, int(np.floor(y)) - 1)
            y1 = y0 + 1
            y2 = y0 + 2
            y3 = min(h-1, y0 + 3)
            
            # 获取16个邻域像素的值
            p = np.zeros((4, 4, 3))
            for m in range(4):
                for n in range(4):
                    xi = x0 + n
                    yi = y0 + m
                    xi = min(max(xi, 0), w-1)  # 边界处理
                    yi = min(max(yi, 0), h-1)
                    p[m, n] = image[yi, xi]
            
            # 计算相对位置
            dx = x - x1
            dy = y - y1

            # 系数
            a = np.zeros((4, 4, channel))
            a[0, 0] = p[1, 1]
            a[0, 1] = -0.5*p[1, 0] + 0.5*p[1, 2]
            a[0, 2] = p[1, 0] - 2.5*p[1, 1] + 2*p[1, 2] - 0.5*p[1, 3]
            a[0, 3] = -0.5*p[1, 0] + 1.5*p[1, 1] - 1.5*p[1, 2] + 0.5*p[1, 3]
            
            a[1, 0] = -0.5*p[0, 1] + 0.5*p[2, 1]
            a[1, 1] = 0.25*p[0, 0] - 0.25*p[0, 2] - 0.25*p[2, 0] + 0.25*p[2, 2]
            a[1, 2] = -0.5*p[0, 0] + 1.25*p[0, 1] - p[0, 2] + 0.25*p[0, 3] + 0.5*p[2, 0] - 1.25*p[2, 1] + p[2, 2] - 0.25*p[2, 3]
            a[1, 3] = 0.25*p[0, 0] - 0.75*p[0, 1] + 0.75*p[0, 2] - 0.25*p[0, 3] - 0.25*p[2, 0] + 0.75*p[2, 1] - 0.75*p[2, 2] + 0.25*p[2, 3]
            
            a[2, 0] = p[0, 1] - 2.5*p[1, 1] + 2*p[2, 1] - 0.5*p[3, 1]
            a[2, 1] = -0.5*p[0, 0] + 0.5*p[0, 2] + 1.25*p[1, 0] - 1.25*p[1, 2] - p[2, 0] + p[2, 2] + 0.25*p[3, 0] - 0.25*p[3, 2]
            a[2, 2] = p[0, 0] - 2.5*p[0, 1] + 2*p[0, 2] - 0.5*p[0, 3] - 2.5*p[1, 0] + 6.25*p[1, 1] - 5*p[1, 2] + 1.25*p[1, 3] + 2*p[2, 0] - 5*p[2, 1] + 4*p[2, 2] - p[2, 3] - 0.5*p[3, 0] + 1.25*p[3, 1] - p[3, 2] + 0.25*p[3, 3]
            a[2, 3] = -0.5*p[0, 0] + 1.5*p[0, 1] - 1.5*p[0, 2] + 0.5*p[0, 3] + 1.25*p[1, 0] - 3.75*p[1, 1] + 3.75*p[1, 2] - 1.25*p[1, 3] - p[2, 0] + 3*p[2, 1] - 3*p[2, 2] + p[2, 3] + 0.25*p[3, 0] - 0.75*p[3, 1] + 0.75*p[3, 2] - 0.25*p[3, 3]
            
            a[3, 0] = -0.5*p[0, 1] + 1.5*p[1, 1] - 1.5*p[2, 1] + 0.5*p[3, 1]
            a[3, 1] = 0.25*p[0, 0] - 0.25*p[0, 2] - 0.75*p[1, 0] + 0.75*p[1, 2] + 0.75*p[2, 0] - 0.75*p[2, 2] - 0.25*p[3, 0] + 0.25*p[3, 2]
            a[3, 2] = -0.5*p[0, 0] + 1.25*p[0, 1] - p[0, 2] + 0.25*p[0, 3] + 1.5*p[1, 0] - 3.75*p[1, 1] + 3*p[1, 2] - 0.75*p[1, 3] - 1.5*p[2, 0] + 3.75*p[2, 1] - 3*p[2, 2] + 0.75*p[2, 3] + 0.5*p[3, 0] - 1.25*p[3, 1] + p[3, 2] - 0.25*p[3, 3]
            a[3, 3] = 0.25*p[0, 0] - 0.75*p[0, 1] + 0.75*p[0, 2] - 0.25*p[0, 3] - 0.75*p[1, 0] + 2.25*p[1, 1] - 2.25*p[1, 2] + 0.75*p[1, 3] + 0.75*p[2, 0] - 2.25*p[2, 1] + 2.25*p[2, 2] - 0.75*p[2, 3] - 0.25*p[3, 0] + 0.75*p[3, 1] - 0.75*p[3, 2] + 0.25*p[3, 3]
            
            # 计算插值结果
            value = np.zeros(channel)
            for m in range(4):
                for n in range(4):
                    value += a[m, n] * (dx**n) * (dy**m)
            
            resized_image[i, j] = np.clip(value, 0, 255)
    
    return resized_image.astype(np.uint8)

# 双立方法插值实现图像旋转
def bicubic_rotation(image, angle):
    h, w, channel = image.shape
    angle_rad = math.radians(angle)
    
    # 计算旋转后的图像尺寸
    cos_theta = abs(math.cos(angle_rad))
    sin_theta = abs(math.sin(angle_rad))
    new_w = int(h * sin_theta + w * cos_theta)
    new_h = int(h * cos_theta + w * sin_theta)
    
    # 旋转中心
    cx, cy = w / 2, h / 2
    new_cx, new_cy = new_w / 2, new_h / 2
    
    rotated_image = np.zeros((new_h, new_w, channel))
    
    for i in range(new_h):
        for j in range(new_w):
            # 将新图像坐标转换回原图像坐标
            x = (j - new_cx) * math.cos(angle_rad) + (i - new_cy) * math.sin(angle_rad) + cx
            y = -(j - new_cx) * math.sin(angle_rad) + (i - new_cy) * math.cos(angle_rad) + cy
            
            if 0 <= x < w and 0 <= y < h:
                # 确定16个邻域像素的坐标
                x0 = max(0, int(np.floor(x)) - 1)
                x1 = x0 + 1
                x2 = x0 + 2
                x3 = min(w-1, x0 + 3)
                
                y0 = max(0, int(np.floor(y)) - 1)
                y1 = y0 + 1
                y2 = y0 + 2
                y3 = min(h-1, y0 + 3)
                
                # 获取16个邻域像素的值
                p = np.zeros((4, 4, channel))
                for m in range(4):
                    for n in range(4):
                        xi = x0 + n
                        yi = y0 + m
                        xi = min(max(xi, 0), w-1)  # 边界处理
                        yi = min(max(yi, 0), h-1)
                        p[m, n] = image[yi, xi]
                
                # 计算相对位置
                dx = x - x1
                dy = y - y1

                # 系数
                a = np.zeros((4, 4, channel))
                a[0, 0] = p[1, 1]
                a[0, 1] = -0.5*p[1, 0] + 0.5*p[1, 2]
                a[0, 2] = p[1, 0] - 2.5*p[1, 1] + 2*p[1, 2] - 0.5*p[1, 3]
                a[0, 3] = -0.5*p[1, 0] + 1.5*p[1, 1] - 1.5*p[1, 2] + 0.5*p[1, 3]
                
                a[1, 0] = -0.5*p[0, 1] + 0.5*p[2, 1]
                a[1, 1] = 0.25*p[0, 0] - 0.25*p[0, 2] - 0.25*p[2, 0] + 0.25*p[2, 2]
                a[1, 2] = -0.5*p[0, 0] + 1.25*p[0, 1] - p[0, 2] + 0.25*p[0, 3] + 0.5*p[2, 0] - 1.25*p[2, 1] + p[2, 2] - 0.25*p[2, 3]
                a[1, 3] = 0.25*p[0, 0] - 0.75*p[0, 1] + 0.75*p[0, 2] - 0.25*p[0, 3] - 0.25*p[2, 0] + 0.75*p[2, 1] - 0.75*p[2, 2] + 0.25*p[2, 3]
                
                a[2, 0] = p[0, 1] - 2.5*p[1, 1] + 2*p[2, 1] - 0.5*p[3, 1]
                a[2, 1] = -0.5*p[0, 0] + 0.5*p[0, 2] + 1.25*p[1, 0] - 1.25*p[1, 2] - p[2, 0] + p[2, 2] + 0.25*p[3, 0] - 0.25*p[3, 2]
                a[2, 2] = p[0, 0] - 2.5*p[0, 1] + 2*p[0, 2] - 0.5*p[0, 3] - 2.5*p[1, 0] + 6.25*p[1, 1] - 5*p[1, 2] + 1.25*p[1, 3] + 2*p[2, 0] - 5*p[2, 1] + 4*p[2, 2] - p[2, 3] - 0.5*p[3, 0] + 1.25*p[3, 1] - p[3, 2] + 0.25*p[3, 3]
                a[2, 3] = -0.5*p[0, 0] + 1.5*p[0, 1] - 1.5*p[0, 2] + 0.5*p[0, 3] + 1.25*p[1, 0] - 3.75*p[1, 1] + 3.75*p[1, 2] - 1.25*p[1, 3] - p[2, 0] + 3*p[2, 1] - 3*p[2, 2] + p[2, 3] + 0.25*p[3, 0] - 0.75*p[3, 1] + 0.75*p[3, 2] - 0.25*p[3, 3]
                
                a[3, 0] = -0.5*p[0, 1] + 1.5*p[1, 1] - 1.5*p[2, 1] + 0.5*p[3, 1]
                a[3, 1] = 0.25*p[0, 0] - 0.25*p[0, 2] - 0.75*p[1, 0] + 0.75*p[1, 2] + 0.75*p[2, 0] - 0.75*p[2, 2] - 0.25*p[3, 0] + 0.25*p[3, 2]
                a[3, 2] = -0.5*p[0, 0] + 1.25*p[0, 1] - p[0, 2] + 0.25*p[0, 3] + 1.5*p[1, 0] - 3.75*p[1, 1] + 3*p[1, 2] - 0.75*p[1, 3] - 1.5*p[2, 0] + 3.75*p[2, 1] - 3*p[2, 2] + 0.75*p[2, 3] + 0.5*p[3, 0] - 1.25*p[3, 1] + p[3, 2] - 0.25*p[3, 3]
                a[3, 3] = 0.25*p[0, 0] - 0.75*p[0, 1] + 0.75*p[0, 2] - 0.25*p[0, 3] - 0.75*p[1, 0] + 2.25*p[1, 1] - 2.25*p[1, 2] + 0.75*p[1, 3] + 0.75*p[2, 0] - 2.25*p[2, 1] + 2.25*p[2, 2] - 0.75*p[2, 3] - 0.25*p[3, 0] + 0.75*p[3, 1] - 0.75*p[3, 2] + 0.25*p[3, 3]
                
                # 计算插值结果
                value = np.zeros(channel)
                for m in range(4):
                    for n in range(4):
                        value += a[m, n] * (dx**n) * (dy**m)
                
                rotated_image[i, j] = np.clip(value, 0, 255)
    
    return rotated_image.astype(np.uint8)

图像插值算法优化：基于四平面

在上述的多种基于分段插值的图像插值算法中，均采用f(i, j)来表示图像的像素点坐标处的颜色值，其中ｉ表示行坐标，ｊ表示列坐标。为进一步地体现图像的局部特征差异并将其用于插值过程，我们引入“平面”的概念，并对图像数据进行升维处理，用三维空间点(i, j, f(i, j))来表示一个像素，并将其对应至空间坐标系中的一个点(x, y, z)。

对一个待插入点而言，可以通过坐标平移将其周围4 个像素点转换为：（注意：此处z0z3为像素坐标点s0s3的颜色值，下同）
$$
s_0(0,0,z_0),s_1(0,1,z_1),s_2(1,0,z_2),s_3(1,1,z_3)
$$
从上述４个点的坐标可以看出它们任意３个点一定不在同一条直线上， 不在同一直线上的３个点可以确定一个平面， 下面讨论具体的插值方法：

1. 先求出这４个点可能的４个平面方程

   已知空间平面的一般方程为：
   $$
   Ax+By+Cz+D=0
   $$
   将s0、s1、s2分别带入上式可得：
   $$
   \begin{cases}Cz_0+D=0\B+Cz_1+D=0\A+Cz_2+D=0&\end{cases}
   $$
   则有：
   $$
   D=-Cz_0,B=C(z_0-z_1),A=C(z_0-z_2)
   $$
   再将其带回空间平面方程，整理后用f(x, y)代替z得到插值公式：
   $$
   f(x,y)=(z_{2}-z_{0})x+(z_{1}-z_{2})y+z_{0}
   $$
   同理，将s0、s1、s3带入空间平面方程可得插值公式：
   $$
   f(x,y)=(z_{3}-z_{1})x+(z_{1}-z_{0})y+z_{0}
   $$
   将s0、s2、s3带入空间平面方程可得插值公式：
   $$
   f(x,y)=(z_{2}-z_{0})x+(z_{3}-z_{2})y+z_{0}
   $$
   将s1、s2、s3带入空间平面方程可得插值公式：
   $$
   f(x,y)=(z_{3}-z_{1})x+(z_{3}-z_{2})y+(z_{2}+z_{1}-z_{3})
   $$

2. 如果s0、s1、s2、s3这４ 个点在同一平面上， 则使用上述任意一个插值公式进行插值均可。 【平面法】

   判断这４个点是否在同一平面上， 只需要比较z1+z2 与 z0+z3是否相等：

   线段s0s3中点坐标为
   $$
   (\frac12,\frac12,\frac{z_0+z_3}2)
   $$
   线段s1s2中点坐标为
   $$
   (\frac12,\frac12,\frac{z_1+z_2}2)
   $$
   如果它们的中点坐标相同，则说明两条线段相交，相交的两条直线可以决定一个平面，即如果待插人点周围的四个点满足：
   $$
   z_1+z_2=z_0+z_3
   $$
   则这它们就是同一平面上的 4 个点，否则就不是同一平面上的 4 个点。

3. 从４个可能的平面中选择一个平面进行插值【四平面法】

   如果它们不是同一平面上的４个点， 情况比较复杂， 需认真讨论，s0、s1、s2、s3４个点的位置关系如下图所示：

   

   在插值的过程中如果一半的区域选择由s0、s1、s2 所确定的平面进行插值， 则另一半必须选择由s1、s2、s3所确定的平面进行插值， 以保证对角线的每一边都是在同一个平面上， 避免出现 “锯齿形” 边缘，为了便于描述， 称s0、s1、s2所确定的平面为 “左下平面”，s1、s2、s3 所确定的平面为 “右上平面”，s0、s1、s3 所确定的平面 “左上平面”，s0、s2、s3所确定的平面 “右下平面”。 为此， 需要参考周围其他点的情况以决定选择哪个平面进行插值。 具体情况如下图所示（黑点是待插入点周围的４个点，白点是参考点）：

   

   1. 对于“左下平面”， 只能参考s0、s1、s2三点左面和下面的点， 即s0、s1、s2三点与s4、s5、s6、s8四个点中的任意一点在同一平面上即可。
   2. 对于“右上平面”，只能参考s1、s2、s3三点右面和上面的点， 即s1、s2、s3三点与s7、s9、s10、s11四个点中的任意一点在同一平面上即可。
   3. 对于 “左上平面”，只能参考s0、s1、s3三点左面和上面的点， 即s0、s1、s3三点与s6、s8、s10、s11四个点中的任意一点在同一平面上即可。
   4. 对于 “右下平面”，只能参考s0、s2、s3三点右面和下面的点， 即s0、s2、s3三点与s4、s5、s7、s9四 个点中的任意一点在同一平面上即可。

   针对1、2两种情况， 当y = 1 + x时，用 “左下平面” 进行插值， 否则用 “右上平面” 进行插值；针对3、4两种情况， 当y = x ^ 3时，用 “左上平面” 进行插值， 否则用 “右下平面” 进行插值。

   判断４个点在同一平面上的方法：（以情况1为例）

   • 对于判断s0、s2、s1、s8 ４ 点是否在同一平面上， 只需要判断z0 + z1与z2 + z8是否相等即可；

   • 对于s0、s1、s2、s5 ４点， 只需要判断z0 + z2与z1 + z5是否相等即可；

   • 对于s0、s1、s2、s6 ４点：如果s0、s2、s6 ３点在同一直线上， 则直线外一点s1与该直线就可以确定一个平面，而要判断这三点是否在同一直线上，只需判断z2 + z6与2 * z0是否相等即可【线段s2(1, 0, z2) s6(-1, 0, z6) 的中点坐标为(0, 0, z2 + z6)，若z2 + z6 = 2 * z0， 则点s0(0, 0, z0)就是它们的中点坐标，当然这３点就在同一条直线上】；

   • 对于s0、s1、s2、s4 ４点， 与s0、s1、s2、s6 ４ 点的情况相同。

4. 如果2和3两点中的情形均不满足， 说明待插入点周围的情况太复杂（不符合平面插值）， 此时采用双线性法进行插值。

基于以上算法思想，编写python函数代码实现图像放缩与旋转过程中的四平面法插值：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104

# 四平面法插值实现图像放缩
def four_plane_interpolation(img, scale):
    H, W, C = img.shape
    new_H, new_W = int(H * scale), int(W * scale)
    output = np.zeros((new_H, new_W, C), dtype=img.dtype)
    
    for y in range(new_H):
        for x in range(new_W):
            # 计算原图对应坐标（浮点数）
            src_x = x / scale
            src_y = y / scale
            
            # 获取周围4个整数坐标点
            x0, y0 = int(np.floor(src_x)), int(np.floor(src_y))
            x1, y1 = min(x0 + 1, W - 1), min(y0 + 1, H - 1)
            
            # 获取4个点的颜色值（z坐标）
            s0 = img[y0, x0]
            s1 = img[y0, x1]
            s2 = img[y1, x0]
            s3 = img[y1, x1]
            
            # 计算相对位置（归一化到[0,1]）
            dx = src_x - x0
            dy = src_y - y0
            
            # 判断四点是否共面（z1 + z2 ≈ z0 + z3）
            if np.allclose(s1 + s2, s0 + s3, atol=1e-6):
                # 共面时，选择任意平面（此处用左下平面）
                interpolated = s0 + (s2 - s0) * dx + (s1 - s0) * dy
            else:
                # 不共面时，动态选择平面
                # 获取周围12个参考点（简化实现，仅取最近邻）
                # 注：论文中需判断参考点是否共面，此处简化逻辑
                if dy > 1 - dx:  # 对角线 y = 1 - x 上方
                    # 选择右上平面
                    interpolated = (s3 - s1) * dx + (s3 - s2) * dy + (s2 + s1 - s3)
                else:
                    # 选择左下平面
                    interpolated = s0 + (s2 - s0) * dx + (s1 - s0) * dy
            
            # 边界检查
            interpolated = np.clip(interpolated, 0, 255 if img.dtype == np.uint8 else 1.0)
            output[y, x] = interpolated
    
    return output

# 四平面法插值实现图像旋转
def four_plane_rotation(image, angle):
    h, w, channel = image.shape
    angle_rad = math.radians(angle)
    
    # 计算旋转后的图像尺寸
    cos_theta = abs(math.cos(angle_rad))
    sin_theta = abs(math.sin(angle_rad))
    new_w = int(h * sin_theta + w * cos_theta)
    new_h = int(h * cos_theta + w * sin_theta)
    
    # 旋转中心
    cx, cy = w / 2, h / 2
    new_cx, new_cy = new_w / 2, new_h / 2
    
    rotated_image = np.zeros((new_h, new_w, channel), dtype=image.dtype)
    
    for i in range(new_h):
        for j in range(new_w):
            # 将新图像坐标转换回原图像坐标
            x = (j - new_cx) * math.cos(angle_rad) + (i - new_cy) * math.sin(angle_rad) + cx
            y = -(j - new_cx) * math.sin(angle_rad) + (i - new_cy) * math.cos(angle_rad) + cy
            
            # 边界检查
            if 0 <= x < w and 0 <= y < h:
                # 获取周围4个整数坐标点
                x0, y0 = int(np.floor(x)), int(np.floor(y))
                x1, y1 = min(x0 + 1, w - 1), min(y0 + 1, h - 1)
                
                # 获取4个点的颜色值
                s0 = image[y0, x0]
                s1 = image[y0, x1]
                s2 = image[y1, x0]
                s3 = image[y1, x1]
                
                # 计算相对位置
                dx = x - x0
                dy = y - y0
                
                # 判断四点是否共面
                if np.allclose(s1 + s2, s0 + s3, atol=1e-6):
                    # 共面时，选择任意平面（此处用左下平面）
                    interpolated = s0 + (s2 - s0) * dy + (s1 - s0) * dx
                else:
                    # 不共面时，动态选择平面
                    if dy > 1 - dx:  # 对角线 y = 1 - x 上方
                        # 选择右上平面
                        interpolated = (s3 - s1) * dx + (s3 - s2) * dy + (s2 + s1 - s3)
                    else:
                        # 选择左下平面
                        interpolated = s0 + (s2 - s0) * dy + (s1 - s0) * dx
                
                # 边界检查
                interpolated = np.clip(interpolated, 0, 255 if image.dtype == np.uint8 else 1.0)
                rotated_image[i, j] = interpolated
    
    return rotated_image

实验测试结果分析

一个理想的插值算法对一幅图像逆时针旋转若干度，再顺时针旋转若干度，应该与原图像相同；同理，对一幅图像放大若干倍，再缩小若干倍，也应该与原图像相同。 基于此，将下面的4幅图像分别用4种算法先逆时针旋转45°，再顺时针旋转45°；先放大４倍，再缩小４倍，然后分别用峰值信噪比（PSNR）验证各算法的优劣。

从定性实验的效果角度，上述四幅图像通过常用的三种分段插值算法完成上述的放大与旋转任务后得到的结果如下图所示：

从实验结果上来看，最近邻算法的边缘颜色“最醒目”，且出现了较为严重的“锯齿形”边缘现象；双线性算法的边缘颜色“最暗淡”；双线性算法和双三次算法也有“锯齿形”边缘现象， 但视觉效果相比最近邻算法而言并不明显。

通过改进的四平面插值算法，对上述四幅图像完成上述的放大与旋转任务后得到的结果如下图所示：

可以看到，四平面插值算法处理后的图像斜线边缘部分是 “光滑连续” 的， 视觉效果比较好，同时有效避免了“锯齿形”边缘现象和“马赛克”现象。

从定量实验的数据角度，我们对于各图像用不同算法完成上述旋转与放缩任务后得到的图像峰值信噪比与算法运行时间进行了计算与统计，结果如下表所示：

峰值信噪比(PSNR)用于表示信号的最大可能功率与影响其表示的保真度的破坏噪声的功率之间的比率。PSNR在图像处理上主要用于量化受有损压缩影响的图像和视频的重建质量。

PSNR 通过均方误差( MSE ) 定义。

给定一个无噪声的m×n单色图像I及其噪声近似值K，MSE定义为：
$$
MSE=\frac{1}{mn}\sum_{i=0}^{m-1}\sum_{j=0}^{n-1}[I(i,j)-K(i,j)]^2.
$$
故PSNR定义为：
$$
\begin{aligned}\mathrm{PSNR}&=10\cdot\log_{10}\left(\frac{MAX_I^2}{MSE}\right)\&=20\cdot\log_{10}\left(\frac{MAX_I}{\sqrt{MSE}}\right)\&=20\cdot\log_{10}(MAX_I)-10\cdot\log_{10}(MSE).\end{aligned}
$$
一般而言，通过PSNR来判断处理后图像的失真情况有如下通用结论：

• PSNR > 30 dB：图像质量较好，失真不明显。
• PSNR 20~30 dB：中等质量，存在可察觉失真。
• PSNR < 20 dB：质量较差，失真显著。

实际计算时，采用opencv自带的PSNR方法cv2.PSNR(img, output)对原始图像与处理后图像的PSNR进行比较计算。

通过对比上述定量实验结果可以发现，在传统的三种分段插值算法中，随着运算阶数（采样待插值点周围的原图像像素点颜色值信息）的增加，图像经过放缩与旋转处理后的失真程度有明显降低，但仍大致处于存在可察觉失真的区间，且算法运行用时也逐渐增加（事实上双立方法的实现可以在编程层面实现优化，这里只是为更直观地展现O（n^2）时间复杂度在图像大小达到一定规模时的显著影响）；而引入的四平面算法不仅在失真程度上较传统的插值算法均有显著改善，算法运行用时也明显优于传统算法中效果最好的双立方法。

综合以上的定性与定量实验结果及分析，本文提出的基于四平面的图像插值算法在图像处理效果（失真）与运行效率上均较传统算法有明显提升，这充分证明了该算法的有效性。

将上文提到的全部四种算法及旋转与放缩两种功能集成到基于python的gui可视化系统中，并打包成exe可执行文件，制作了一个基于插值的图像处理系统，基本功能演示如下图所示：

参考文献

[1] 王开荣,杨大地编著.应用数值分析[M].高等教育出版社,2010.

[2] 毛伟伟,于素萍,石念峰.一种基于四平面的图像插值算法[J].洛阳理工学院学报(自然科学版),2024,34(01):76-81.

[3] 刘显德,李笑.任意大小图像的量子描述及双线性插值方法[J].计算机工程与设计,2024,45(08):2423-2432.

[4] 张喜民,詹海生.基于双三次插值的Canny-Devernay亚像素图像边缘检测算法[J].现代制造工程,2025,(03):107-114.

[5] 陈玲玲,周宁,殷永,等.插值方法在光声图像重建中的应用[J].计算机与数字工程,2013,41(10):1676-1677+1694.

1 项目背景与研究意义

1.1 项目背景

肺炎作为一种常见的呼吸系统疾病，对人类健康构成了长期威胁，特别是随着COVID-19新冠肺炎疫情的全球爆发，其公共健康影响更为显著。COVID-19肺炎具有传播速度快、感染范围广、诊断难度大的特点，对全球医疗系统和社会经济产生了深远影响。特别是在疫情高峰期，医疗资源的短缺和诊断效率的瓶颈，进一步突显了快速、准确诊断工具的重要性。

传统肺炎诊断方法主要依赖医生对胸部X光片或CT图像的人工分析，既耗时又容易受到经验和疲劳的影响，尤其在COVID-19疫情期间，大量影像数据的涌现使得人工诊断难以满足需求。与此同时，COVID-19的影像表现与其他类型肺炎的重叠性增加了诊断的复杂性，这进一步加剧了对智能化诊断系统的需求。

随着人工智能技术的快速发展，深度学习为医疗影像分析带来了全新的解决方案。卷积神经网络（Convolutional Neural Network, CNN）凭借其强大的图像特征提取能力，在自动化诊断中展现了巨大潜力。同时，由于不同医疗机构的CT扫描设备性能差异显著，特别是在医疗资源较为匮乏的地区，CT影像常常受到设备老化、分辨率低或操作不规范等因素的影响，图像质量参差不齐，这不仅增加了诊断的复杂性，还对自动化系统的鲁棒性提出了更高要求；而自编码器（Autoencoder）作为一种有效的降噪工具，为医疗影像数据预处理提供了重要支持。通过自编码器的引入，可以有效消除图像中的噪声干扰，减少不同设备间的成像差异，为后续的分类和识别模型提供高质量的输入数据。因此，本项目提出结合自编码器和卷积神经网络的深度学习框架，开发一套针对肺炎（包括COVID-19）CT影像识别的智能诊断系统。

1.2 研究意义

本项目以新型冠状病毒肺炎为切入点，面向未来医学智能化需求，开发的诊断系统不仅能够应对当下疫情挑战，还具有推广至其他医学影像诊断场景的潜力，从而为全球公共健康事业的发展提供有力支持。

面向新冠肺炎疫情期间大规模肺部影像数据的快速诊断需求，本系统依托深度学习技术，有效提升了诊断效率，为疫情防控和患者管理提供重要的技术支持。通过先进的模型算法，系统能够精准识别肺部CT影像中的病变特征，减少人为误差，确保诊断结果的准确性和一致性，大幅降低不同医疗机构和医生之间的诊断差异，避免误诊和漏诊风险，从而更好地保障患者安全。

针对基层医院或偏远地区医疗资源匮乏的现状，本系统可作为一种可靠的辅助诊断工具，为医生提供科学的决策支持，帮助缓解诊断能力不足带来的压力。其高效的处理能力不仅提高了基层医疗服务水平，也为疫情防控的全面推进提供了技术保障，为应对紧急医疗需求的地区解决实际困难。

此外，本系统在高效处理和分析海量肺部CT影像数据的基础上，还为研究新冠肺炎的病理特征及流行规律提供了宝贵的数据支持。这些分析结果可进一步应用于疫情传播趋势预测和公共卫生政策制定，为疫情防控策略的科学性和有效性奠定了坚实基础。

同时，本系统也为患者病情的动态管理提供了重要帮助。通过智能分析新冠肺炎影像特征的变化趋势，系统能够为临床医生提供精确的病情评估建议，有助于及时调整治疗方案。这种基于影像数据的技术支持，不仅提高了患者管理的科学性和有效性，还为医疗资源的合理分配提供了重要依据，进一步推动了疫情防控工作的高效开展。

2 数据集获取与预处理

2.1 数据集介绍

本项目使用的肺部X-光片数据集从Kaggle网站（链接：https://www.kaggle.com/datasets/alsaniipe/chest-x-ray-image）获取，共分为三类标签：新冠肺炎COVID19、正常NORMAL和普通肺炎PNEUMONIA，该数据集已经事先划分好了训练集与测试集。此外，为测试模型对于不同质量CT影像的识别精度，还对测试集中的部分图像使用高斯噪声进行扰动，以模拟实际CT扫描的成像质量差异。

数据分布如下表所示：

表2.1：数据集数据分布情况

2.2 编程环境搭建

本项目中所有的代码编写与运行均是在配备NVIDIA GeForce RTX 3060显卡、16GB运行内存、12th Gen Inter(R) Core(TM) i7-12700@2.10GHz处理器与Microsoft Windows 11操作系统的工作站上使用Python编程语言完成的。

软件环境方面，采用Conda进行环境管理，在控制台中通过命令”conda create -n covid python=3.12”创建虚拟环境，并在激活环境后使用pip install命令依次安装所需的各种依赖库；全部安装并测试完成后，通过命令”pip freeze >
requirements.txt”将虚拟环境中安装的所有依赖库及对应版本写入文件requirements.txt，后续移植时可在新的运行环境中运行命令”pip install -r requirements.txt”完成环境的一键配置。

值得注意的是，在PyTorch（包括torch、torchvision与torchaudio库）安装时需要根据自己电脑使用的CUDA版本（使用CPU则直接在命令行中使用pip install安装即可）在PyTorch官网中找到对应的安装命令进行安装。我使用的CUDA版本为12.6，安装命令的选取如下图所示：

下面对于项目中使用到的主要依赖库进行简要介绍：

其中，plt用于绘图，nn中包含了用于构建神经网络的隐藏层（全连接层、卷积层等），F中包含了各种激活函数（ReLu、Sigmoid等），DataLoader用于在训练时加载数据，datasets和transforms用于读取和处理数据集，tqdm用于进度条可视化，torchmetrics用于模型精度的测试，torchviz和torchsummary用于以图形与文字的方式描述模型架构概况。

2.3 图像数据读取

torch对于一些常用的数据集做了封装，可以直接调用，例如datasets.MNIST()。但此处我们使用的是本地的图片数据，可以使用ImageFolder将一个文件夹下的图片读取成数据集并完成数据增强工作。在读取完数据集后，还需要定义DataLoader用于加载数据为可分批次（batch）读取的迭代器以供后续使用。为使得代码更加简洁，将上述的数据读取与加载过程为封装在getDataLoader函数中，并在主函数中通过指定不同的目录加载训练集、测试集或是含噪声测试集。

可以看到，其中构建了数据增强器transform，在读取数据时进行相应处理：

• Grayscale: 指以灰度图的形式读取。

• Resize: 由于图像尺寸各不相同，在训练前需将它们重塑成相同尺寸256*256。

• ToTensor: 将图片格式转换成张量形式，torch的计算以张量的形式进行。

除此之外，在构建数据加载器时需要指定一个批次（batch）中的图片数据数量batch_size，在模型训练时训练批次大小TRAIN_BATCH_SIZE也是会影响最终模型性能的重要超参数之一。在训练过程中，设定TRAIN_BATCH_SIZE为32，而在测试过程中，为提高测试效率，将TEST_BATCH_SIZE设置为66并对函数进行对应修改。

2.4 叠加噪声函数

不论是构建噪声测试集，还是在利用无噪声的训练集进行训练时，都需要手动添加噪声，故编写add_noise函数，默认的噪声强度为0.5，并在添加噪声后进行归一化以确保图像值位于[0,1]范围内。

3 模型构建与网络训练

3.1 整体模型框架

整体模型框架由两个核心部分组成，分别是用于去噪的数据预处理模块和负责分类的卷积神经网络（CNN）。去噪模块采用自编码器（Autoencoder）的架构，专注于从输入数据中去除噪声，以提升后续分类的准确性；分类模块基于卷积神经网络，其强大的特征提取和模式识别能力使其成为分类任务的理想选择。

两个模块相辅相成，通过有效的数据处理和特征提取，确保模型能够在噪声干扰较大的环境中实现高精度分类。噪声数据首先经过自编码器处理，生成质量优化的特征表示，然后被CNN接收并完成分类任务。这一整体框架设计非常适合肺炎图像识别任务，通过结合去噪和分类两大模块的优势，模型不仅能够有效提高数据质量，还能充分挖掘数据中的有用特征，从而能够在复杂的医学影像处理中表现出卓越的鲁棒性和准确性，满足肺炎诊断的实际需求。

3.2 自编码器

自编码器模型用于处理输入数据中的噪声问题，提升后续分类的准确性。其核心思想是通过编码器将输入数据压缩至低维潜在表示（latent representation），再由解码器将其还原至去噪后的重构数据，从而实现降噪效果。

3.2.1 网络结构设计

自编码器网络结构由编码器encoder与解码器decoder组成：

• 编码器由两层卷积（Conv2d）和两次池化（MaxPool2d）操作组成，用于提取特征；

• 解码器通过两次反卷积（ConvTranspose2d）和两次上采样（UpsamplingNearest2d）逐步恢复图像尺寸到原始大小；

• 最后使用Sigmoid激活函数将输出值限制在[0,1]区间。

模型定义代码如下：

模型继承自nn.Module类，在__init__()函数中定义模型的结构，在forward()函数中定义模型的前向传播过程。

通过调用torchviz和torchsummary库，可以输出该模型结构的基本信息：

3.2.2 模型训练

基本的训练流程集成在函数train_autoencoder_process中，如下图所示：

其中指定优化器optimizer为Adam，损失函数为均方误差MSE，并使用超参数：训练轮数Epochs=50、学习率lr=0.001。每轮（Epoch）训练中均需要以多个batch的形式遍历训练集中的所有数据，并在每个batch后对模型进行更新，具体而言每次更新均需执行如下操作：

• 从加载器中获取输入数据

• 使用add_noise函数对干净图像加噪

• 将加噪后图像输入自编码器模型并计算模型输出

• 根据模型输出和标签计算损失Loss

• 清空梯度

• 反向传播

• 更新模型

值得注意的是，由于用于训练的图像数据没有噪声，因此训练时首先需要对输入的图像进行加噪处理，再输入自编码器模型进行训练。

训练过程中还利用tqdm进度条函数对训练进程进行可视化，并在每轮训练完成后打印出当轮训练过程中模型的平均损失：

在训练过程中，将每轮训练的平均损失存储在列表中，并在训练结束后将平均损失的变化过程以图像形式呈现：

可以看到，经过多轮训练，模型的损失函数值在不断减小且逐渐趋近于0，这意味着该自编码器的模型训练过程是收敛的，模型具有较稳定的工作性能。

3.3 卷积神经网络

卷积神经网络负责从图像中提取多层次的空间特征，通过逐步减少图像尺寸和增加特征通道来捕捉关键信息，从而实现去噪后肺部CT图像的分类功能。CNN以其强大的特征提取能力，能够有效处理图像的局部依赖性和空间不变性，高效处理结构化数据（如图像、时序数据）。模型简单且高效，具有较强的泛化能力，适合处理小规模数据集的图像分类问题。

3.3.1 网络结构设计

卷积神经网络结构（如上图，通过NN-SVG工具绘制）由两层卷积层（Conv2d）和池化层（MaxPool2d）组成，激活函数均选用ReLU，逐步提取特征并将输入图像的尺寸从原始大小减小到64×64。卷积后的特征图展平后通过三个全连接层（Linear），分别将特征维度从32×64×64降至128，再降至32，最后输出3个类别（Covid19、Normal、Pneumonia）的预测结果。

模型定义代码如下：

模型继承自nn.Module类，在__init__()函数中定义模型的结构，在forward()函数中定义模型的前向传播过程。

通过调用torchviz和torchsummary库，可以输出该模型结构的基本信息：

3.3.2 模型训练

基本的训练流程集成在函数train_cnn_process中，如下图所示：

其中指定优化器optimizer为Adam，损失函数为交叉熵损失CrossEntropy，并使用超参数：训练轮数Epochs=50、学习率lr=0.001。每轮（Epoch）训练中均需要以多个batch的形式遍历训练集中的所有数据，并在每个batch后对模型进行更新，具体而言每次更新均需执行如下操作：

• 从加载器中获取输入数据

• 使用add_noise函数对干净图像加噪

• 将加噪后图像输入训练好的自编码器模型trained_autoencoder_model

• 将经过自编码器去噪后的图像输入CNN模型并计算模型输出

• 根据模型输出和标签计算损失Loss

• 清空梯度

• 反向传播

• 更新模型

值得注意的是，由于用于训练的图像数据没有噪声，为与实际的输入情况一致，首先需要对输入的图像进行加噪处理，再利用训练好的自编码器模型进行降噪（为了不在更新CNN的同时更新自编码器，这一步不需要产生梯度），才能输入CNN分类模型进行训练。

训练过程中还利用tqdm进度条函数对训练进程进行可视化，并在每轮训练完成后打印出当轮训练过程中模型的平均损失与在训练集上的测试精度：

在训练过程中，将每轮训练的平均损失与模型在训练集上的测试精度存储在列表中，并在训练结束后将两者的变化过程以图像形式呈现：

可以看到，经过多轮训练，模型的损失函数值在不断减小且逐渐趋近于0，这意味着该自编码器的模型训练过程是收敛的，模型具有较稳定的工作性能；同时随着训练轮数增加，模型在训练集上的精度也逐渐增高（波动上升），在模型训练完成时，卷积神经网络在训练集上的分类精度已经可以达到99.59%（一度达到99.90%），接近百分之百，说明模型的分类能力较好。

4 模型测试及应用

4.1 自编码器降噪效果

在自编码器模型的训练过程中，每隔10轮对模型参数进行了一次存档；在测试过程中，分别使用训练轮数为10、20、30、40、50的自编码器模型对于加噪后的模型进行降噪处理，效果如下图所示：

通过对比不同训练轮数的自编码器模型降噪效果可以发现，随着训练轮数的增加，自编码器模型的降噪效果在逐渐提升，但在Epoch到达30之后，训练带来的降噪效果提升就不如先前显著了。尽管由于较大的噪声强度（0.5）导致降噪后的图像仍然比较模糊，但通过肉眼还是能粗略观察处肺部骨骼的轮廓等特征，后续实验也证明了卷积神经网络确实可以从这样清晰度的图像中提取相应的特征来进行分类，该自编码器模型的设计有效。

4.2 卷积神经网络分类精度

在卷积神经网络的分类精度上，训练过程中已经实时对于每一轮训练后的模型在训练集上进行了精度测试（3.2.2节中已有提及），而在测试集上，可以编写与训练过程类似的代码利用torchmetrics库对模型分类精度进行测试，只是不会更新模型，代码如下：

可以看到，由于我们的测试集分为含噪声和不含噪声两类，因此编写了不同的函数对模型分类精度进行测试。两个函数的主要差别就在于，由于含噪声测试集是已经加噪的图片（噪声与手动通过add_noise函数添加的不同），因此在含噪声测试集的测试代码中不必再次手动添加噪声，而是直接将图像输入自编码器降噪后再输入CNN分类模型中进行分类；而对于不含噪声的测试集而言，为模拟与训练集同样的处理流程，会先进行手动加噪再通过自编码器降噪之后才输入CNN分类模型中进行分类。

运行测试代码后，得到模型在含噪测试集上的分类精度为96.97%，在不含噪声的测试集上的分类精度为94.57%，在两个测试集上的分类精度水平均较高，说明该模型具有良好的分类效果。

4.3 模型应用：基于CT影像的肺炎诊断Web服务

通过对比多组超参数的模型降噪与分类效果，最终选定如下的超参数：

• 训练轮数Epochs=50；

• 学习率LR=0.001；

• 训练批次大小Train_Batch_Size=32。

选定参数后，将整体代码抽离为model.py（包含模型定义类代码），run.py（服务端代码）和train.py（训练函数），并将模型部署到实际应用中，使用Flask作为服务端，以Web形式用户提供操作接口以上传图片进行诊断。由于主要功能是提供接口，故网页只做了很简易的一个index.html，给用户提供上传图片的按钮，并在用户上传有噪声的CT影像后返回诊断结果及去噪后的图像。除此之外，还将挂载在本地端口上的Web通过内网穿透映射到公网，以供实时访问。

网页初始界面如下图所示：

接下来分别测试当输入COVID19、NORMAL和PNEUMONIA三个组别的图片，模型能否正确判断：

可以发现，模型在大多数情况下可以正确识别图像来源，但也会出现错误识别的情况，这和Test 集上的Accuracy相符合；此外，在测试时还注意到，模型识别结果偶尔会出现不稳定的现象，即输入同一张图像有时识别为某一类别，有时又会识别为另一类别，这是由模型内部部分随机参数导致的，这也反映了模型在一些模棱两可的情况下（两类别概率接近）做出判断时的不稳定性。在实际应用中，为尽可能减少误诊对于患者带来的各方面影响，还需要采取更多优化措施提升模型性能，并对模型在模棱两可的情况下做出的判断进行合理的限制。

5 总结与展望

本项目全部代码（不包含数据集）已上传至Github仓库，仓库URL地址：https://github.com/Asgard-Tim/Pneumonia-Image-Recognition

5.1 项目总结

本项目基于深度学习技术，结合自编码器和卷积神经网络，开发了一套智能诊断系统，用于快速、高效地识别肺部的CT影像并判断该患者是否患有肺炎（包括COVID-19）。自编码器模块有效去除了噪声，提升了图像质量，而卷积神经网络以其强大的特征提取能力，实现了高精度的分类。本项目在数据预处理、模型设计、网络训练及测试等环节中均采用了创新性的技术方案，最终实现了在含噪声测试集上96.97%和在无噪声测试集上94.57%的分类精度，表现出了较高的鲁棒性和实用价值。同时，系统已通过Flask框架部署为Web服务，能够实时接收CT影像并给出诊断结果，为疫情期间大规模影像数据的快速诊断及基层医疗资源匮乏地区的医疗支持提供了重要的技术保障。

5.2 课程收获与反思

本次选修《智能图像处理》这门课程确实让我学到了很多东西，其实自己之前也自己看过一些机器学习方面的内容，有一定的知识基础与环境搭建经验，但由于各方面原因总是没有系统性的去学习计算机视觉的相关知识，也缺乏足够的实战代码与项目经验。通过这门课程的学习，很大程度上锻炼了我Python的代码能力，也在Coding的过程中不断熟悉OpenCV、Pytorch等库的使用，更在实践的过程中不断加深对于各种算法模型（AlexNet、ResNet、YOLO等）的理解。

本次项目让我完整地经历了从数据集获取、论文调研及算法代码实现，再到代码调试与模型训练测试，最终将模型应用到实际系统中的全过程，在项目实现的过程中收获了很多课程教学与实验中涉及不到的东西，包括数据集的收集、模型的选择以及作为一个完整项目的代码实现等等多个方面，这也是我第一次使用GPU资源去进行。虽然由于时间等条件的限制，在模型选择上并没有进行深入的调研与充分的对比试验，只是基于自己已知的一些知识对于架构较为简单的自编码器模型与卷积神经网络进行了复现与设计，最终模型的分类精度还有一定的提升空间，但是这也为我后续的自主学习打下了一个良好的基础，希望未来我能在计算机视觉方面有更加深入的学习与探索，也感谢老师的耐心指导与悉心教学。

参考文献

[1] Nosa-Omoruyi M, Oghenekaro L U. AutoEncoder Convolutional Neural Network for Pneumonia Detection[J]. arXiv preprint arXiv:2409.02142, 2024.

[2] Ratiphaphongthon W, Panup W, Wangkeeree R. An improved technique for pneumonia infected patients image recognition based on combination algorithm of smooth generalized pinball SVM and variational autoencoders[J]. IEEE Access, 2022, 10: 107431-107445.

[3] Gayathri J L, Abraham B, Sujarani M S, et al. A computer-aided diagnosis system for the classification of COVID-19 and non-COVID-19 pneumonia on chest X-ray images by integrating CNN with sparse autoencoder and feed forward neural network[J]. Computers in biology and medicine, 2022, 141: 105134.

[4] García-Ordás M T, Benítez-Andrades J A, García-Rodríguez I, et al. Detecting respiratory pathologies using convolutional neural networks and variational autoencoders for unbalancing data[J]. Sensors, 2020,20(4): 1214.

[5] Xia Y. Enhanced Pneumonia Detection in Chest X-Rays Based on Integrated Denoising Autoencoders and Convolutional Neural Networks[J].

[6] El-Shafai W, El-Nabi S A, El-Rabaie E S M, et al. Efficient Deep-Learning-Based Autoencoder Denoising Approach for Medical Image Diagnosis[J]. Computers, Materials & Continua, 2022, 70(3).

[7] Rana N, Marwaha H. Auto encoder-guided Feature Extraction for Pneumonia Identification from Chest X-ray Images[C]//E3S Web of Conferences. EDP Sciences, 2024, 556: 01011.

[8] Ankayarkanni B, Sangeetha P. An Autoencoder-BiLSTM framework for classifying multiple types of lung diseases from CXR images[J]. Multimedia Tools and Applications, 2024: 1-30.

[9] 孙敬,丁嘉伟,冯光辉.一种基于自编码器降维的神经卷积网络入侵检测模型[J/OL].电信科学,1-7[2025-01-05].

[10] 张淙越,杨晓玲.基于卷积神经网络的新冠肺炎CT图像识别系统[J].电脑与信息技术,2022,30(03):12-14+40.