基于Python卷积的N体模拟数据合成光谱高效生成方法求助

背景与现有方案痛点 #

处理N体模拟数据时，需要生成单个像素或Voronoi bin区域内粒子的合成光谱，每个粒子都对应一条视线（LOS）速度。当前采用的方案是：对每个粒子的恒星模板光谱做Doppler位移，再插值回原始波长点后求和通量。但面对800万粒子的规模时，该方法运行速度极慢，亟需优化。

尝试路径与卡点 #

为提升效率，尝试改用卷积方法：已通过直方图拟合得到视线速度分布（LOSVD）的Gauss-Hermite函数，但在将其与恒星模板结合实现卷积时遇到核心问题：

• 单位匹配障碍：已知对数波长均匀采样对应速度均匀采样，但不知如何在代码中实际完成这一单位转换
• 代码实现报错：尝试的代码因数组形状不匹配无法运行

需要针对以下需求提供帮助：

• 完整的卷积实现代码（含单位转换逻辑）
• 卷积原理中单位转换的详细说明

现有方法代码（示例） #

# 现有逐粒子处理的示例代码
import numpy as np
from scipy.interpolate import interp1d

def doppler_shift(wave, flux, v_los):
    c = 299792.458  # 光速，单位km/s
    shift_factor = np.sqrt((1 + v_los/c)/(1 - v_los/c))
    shifted_wave = wave * shift_factor
    # 插值回原波长网格
    interpolator = interp1d(shifted_wave, flux, bounds_error=False, fill_value=0)
    return interpolator(wave)

# 假设wave是原始波长网格，template是恒星模板光谱，v_los_list是所有粒子的LOS速度
wave = np.logspace(np.log10(4000), np.log10(7000), 1000)
template = np.random.randn(1000)  # 示例模板
v_los_list = np.random.randn(8_000_000) * 200  # 示例LOS速度

# 逐粒子处理（慢）
total_flux = np.zeros_like(wave)
for v in v_los_list:
    shifted_flux = doppler_shift(wave, template, v)
    total_flux += shifted_flux

尝试的卷积代码（报错版本） #

import numpy as np
from scipy.signal import convolve

# 假设已拟合得到Gauss-Hermite参数：v_mean, sigma, h3, h4
v_mean = 50  # km/s
sigma = 100  # km/s
h3 = 0.1
h4 = 0.05

# 生成Gauss-Hermite LOSVD
def gauss_hermite(v, v_mean, sigma, h3, h4):
    x = (v - v_mean)/sigma
    gauss = np.exp(-0.5*x**2) / (sigma * np.sqrt(2*np.pi))
    hermite = 1 + h3*np.polyval([1,0,-3], x)/np.sqrt(6) + h4*np.polyval([1,0,-6,0,3], x)/np.sqrt(24)
    return gauss * hermite

# 问题：速度网格与波长网格的单位转换
v_grid = np.linspace(-500, 500, 1000)
losvd = gauss_hermite(v_grid, v_mean, sigma, h3, h4)

# 尝试卷积（形状不匹配报错）
convolved_flux = convolve(template, losvd, mode='same')

相关公式 #

1. Doppler位移公式：
   $$\lambda' = \lambda \sqrt{\frac{1 + v_{\text{los}}/c}{1 - v_{\text{los}}/c}} \approx \lambda (1 + v_{\text{los}}/c) \quad (v_{\text{los}} \ll c)$$

2. 对数波长与速度的关系：
   $$\ln\lambda' = \ln\lambda + \ln\left(\sqrt{\frac{1 + v/c}{1 - v/c}}\right) \approx \ln\lambda + v/c$$
   当波长网格为均匀对数采样时，$\Delta\ln\lambda$为常数，对应$\Delta v = c \cdot \Delta\ln\lambda$，即速度网格为均匀采样。

3. Gauss-Hermite LOSVD公式：
   $$f(v) = \frac{1}{\sigma\sqrt{2\pi}} e^{-(v - v_0)^2/(2\sigma2)} \left[ 1 + h_3 H_3\left(\frac{v - v_0}{\sigma}\right) + h_4 H_4\left(\frac{v - v_0}{\sigma}\right) \right]$$
   其中$H_3(x) = x^3 - 3x$，$H_4(x) = x^4 - 6x^2 + 3$为厄米多项式。

内容的提问来源于stack exchange，提问作者mamark