pylops

PyLops - 线性算子库

快速参考

import numpy as np
import pylops

# Create operator and apply forward/adjoint
A = pylops.FirstDerivative(n=100, dtype='float64')
y = A @ x        # Forward: y = A @ x
x_adj = A.H @ y  # Adjoint: x = A.H @ y
x_est = A / y    # Solve inverse problem

关键类

类	用途
LinearOperator	所有算子的基类
VStack/HStack	垂直／水平堆叠算子
BlockDiag	块对角算子组合

基本操作

基本算子

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

# Derivatives
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      # Convolve
x_adj = C.H @ y  # Correlation (adjoint)

组合和堆叠算子

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

# Stack operators
V = pylops.VStack([A, B])     # Vertical: (2n, n)
H = pylops.HStack([A, B])     # Horizontal: (n, 2n)
BD = pylops.BlockDiag([A, B]) # Block diagonal: (2n, 2n)

求解逆问题

# Simple least squares
x_est = A / y

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

# Regularized inversion with smoothness
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

# Deconvolve (regularized)
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)

与替代方案的适用场景对比

场景	建议
面向大型逆问题的无矩阵线性算子	PyLops - 专为此用途构建，内存效率高
具有已知结构的稀疏矩阵运算	scipy.sparse - 标准且文档完善
简单的卷积／反卷积	PyLops - 通过 Convolve1D 提供简洁 API
小型问题的自定义算子	Custom NumPy/SciPy - 无需额外依赖
GPU 加速线性代数	PyLops - 传入 CuPy 数组即可自动使用 GPU
地震反卷积或成像算子	PyLops - 丰富的信号处理算子库

在以下情况下选择 PyLops：需要可扩展到大型问题的无矩阵线性算子， 无需构造显式矩阵。其算子代数（@、VStack、BlockDiag）和 内置求解器（LSQR、FISTA、Split Bregman）可简化逆问题工作流。

在以下情况下避免使用 PyLops：问题规模足够小，可以使用显式矩阵（请使用 NumPy/SciPy）， 或者需要非线性算子（PyLops 严格来说仅支持线性算子）。

常见工作流

正则化地震反卷积

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

提示

1. 切勿构造完整矩阵 - 使用 .matvec() 和 .rmatvec() 提高内存效率
2. 检查形状 - 算子具有类似矩阵的 .shape 属性
3. 验证伴随算子 - 自定义算子务必运行 pylops.utils.dottest()
4. 从简单开始 - 扩大规模前，先在小型问题上测试
5. 使用 GPU - 传入 CuPy 数组即可自动启用 GPU 加速

参考资料

• 算子参考 - 可用算子的完整列表
• 求解器参考 - 求解方法和配置选项

脚本

• scripts/deconvolution.py - 地震反卷积示例