pylops

用途

PyLops 是一个线性算子库，专为需要无矩阵（matrix-free）表示的大规模反问题设计。它不显式构造矩阵，从而在内存受限时仍能处理大规模问题。适合在代理需要定义线性算子、执行正演/伴随运算、组合算子或求解反问题时使用。

主要能力

基本算子

可创建并应用各类基本算子，支持正向与伴随运算：

import numpy as np
import pylops

# 创建算子并执行正演/伴随
A = pylops.FirstDerivative(n=100, dtype='float64')
y = A @ x        # 正演：y = A @ x
x_adj = A.H @ y  # 伴随：x_adj = A.H @ y
x_est = A / y    # 求解反问题

对角算子：

D = pylops.Diagonal(np.array([1., 2., 3.]))
y = D @ x; x_adj = D.H @ y

导数与梯度：

D1 = pylops.FirstDerivative(n, dtype='float64')
D2 = pylops.SecondDerivative(n, dtype='float64')
G = pylops.Gradient(dims=(64, 64), dtype='float64')

卷积：

wavelet = np.sin(np.linspace(0, 2*np.pi, 21)) * np.hanning(21)
C = pylops.signalprocessing.Convolve1D(n, h=wavelet, offset=10)
y = C @ x        # 卷积
x_adj = C.H @ y  # 相关（伴随）

算子组合与堆叠

链式组合：

# y = C @ B @ A @ x
composed = pylops.Smoothing1D(5, n) @ pylops.FirstDerivative(n) @ pylops.Identity(n)

算子堆叠：

V = pylops.VStack([A, B])     # 垂直堆叠：(2n, n)
H = pylops.HStack([A, B])     # 水平堆叠：(n, 2n)
BD = pylops.BlockDiag([A, B]) # 块对角：(2n, 2n)

关键类包括：LinearOperator（所有算子的基类）、VStack/HStack（垂直/水平算子堆叠）、BlockDiag（块对角算子组合）。

求解反问题

简单最小二乘：

x_est = A / y

正规方程：

x_est = pylops.optimization.leastsquares.NormalEquationsInversion(A, None, y)

带平滑度的正则化反演：

Reg = pylops.SecondDerivative(n)
x_est = pylops.optimization.leastsquares.RegularizedInversion(
    A, [Reg], y, epsRs=[0.1]
)

迭代求解器：

x_lsqr = pylops.optimization.solver.lsqr(A, y, iter_lim=100)[0]
x_cgls = pylops.optimization.solver.cgls(A, y, niter=100)[0]

促进稀疏性的反演：

x_l1 = pylops.optimization.sparsity.fista(A, y, niter=100, eps=0.1)[0]

验证伴随（点积测试）

A = pylops.FirstDerivative(100)
pylops.utils.dottest(A, 100, 100, verb=True)  # Dot test passed!

常见模式

地震反卷积

C = pylops.signalprocessing.Convolve1D(n, h=wavelet, offset=len(wavelet)//2)
seismic = C @ reflectivity

# 反卷积（正则化）
Reg = pylops.SecondDerivative(n)
reflectivity_est = pylops.optimization.leastsquares.RegularizedInversion(
    C, [Reg], seismic, epsRs=[0.01]
)

使用 TV 进行图像去噪

ny, nx = image.shape
G = pylops.Gradient(dims=(ny, nx))
x_tv = pylops.optimization.sparsity.splitbregman(
    pylops.Identity(ny*nx), G, image.ravel(),
    niter_inner=5, niter_outer=10, mu=1.0, epsRL1s=[0.1]
)[0].reshape(ny, nx)

使用方式与工作流程

正则化地震反卷积典型流程：

定义子波数组并创建 Convolve1D 算子。
生成或加载地震道数据。
运行点积测试验证算子伴随：pylops.utils.dottest()。
设置正则化算子（如 SecondDerivative 用于平滑）。
运行 RegularizedInversion(C, [Reg], data, epsRs=[eps])。
如有真实值，比较估计反射率与真实值。
调节 epsRs 参数：值越高越平滑，越低越锐利。
若需稀疏解，改用 pylops.optimization.sparsity.fista()。

与替代方案的比较

无矩阵线性算子用于大规模反问题：PyLops——专用、内存高效。
已知结构的稀疏矩阵运算：scipy.sparse——标准、文档完善。
简单卷积/反卷积：PyLops——Convolve1D 提供简洁 API。
小规模问题的自定义算子：NumPy/SciPy——无额外依赖。
GPU 加速线性代数：PyLops——传入 CuPy 数组即可自动使用 GPU。
地震反卷积或成像算子：PyLops——丰富的信号处理算子库。

选择 PyLops 的场景：需要无矩阵线性算子以扩展到大规模问题而不显式构造矩阵；其算子代数（@、VStack、BlockDiag）与内置求解器（LSQR、FISTA、Split Bregman）使反问题工作流简洁。

避免 PyLops 的场景：问题规模小到可用显式矩阵（使用 NumPy/SciPy）；或需要非线性算子（PyLops 严格为线性）。

使用提示

不要构造完整矩阵——使用 .matvec() 和 .rmatvec() 以节省内存。
检查形状——算子有类似矩阵的 .shape 属性。
验证伴随——对自定义算子始终运行 pylops.utils.dottest()。
从简单开始——先在小问题上测试再扩展。
使用 GPU——传入 CuPy 数组即可自动 GPU 加速。

适用场景与限制

适用场景包括：地震反卷积、成像算子等地球物理应用，以及需要无矩阵线性算子的大规模反问题，尤其是涉及卷积、导数、梯度与正则化反演的工作流。

限制：PyLops 严格为线性算子库，不适用于非线性算子需求；对于小规模问题，直接使用 NumPy/SciPy 显式矩阵更合适。

脚本示例

scripts/deconvolution.py —— 地震反卷积示例。

参考

Operators Reference —— 可用算子的完整列表。
Solvers Reference —— 求解器方法及配置选项。