AGN光变曲线Structure Function计算参数相关异常问题排查

项目背景 #

正在开展活动星系核（AGNs）相关研究，通过模拟光变曲线并应用Structure Function（参考Macleod 2012、Kozlowski 2017）进行分析。模型采用带泊松等待时间的爆发函数，每个爆发对应自定义脉冲形状（如指数型、幂律型），生成模拟光变曲线后计算Structure Function。

模型参数设置 #

• 核心控制参数：指数脉冲衰减时间td、爆发平均间隔时间ts，当前将二者设为等值（如ts=20、td=20）
• 时间步长：取td/ts的1%，示例中为0.2
• Structure Function时滞设置：
  • 用np.linspace生成时滞序列k，范围为kmin=timestep*2到kmax=td*30，示例为[0.4,0.6,...,600]
  • 将k转换为相对时间步长差值：lag = k/timestep，示例中为[2,3,...,3000]

核心问题 #

当td/ts取值在8-20区间时，Structure Function的双对数图会出现偏离平滑曲线的异常点，通常对应时滞序列的第5和第10个取值；提高模拟分辨率后，异常现象反而更严重。已通过截取中间1/3数据修正边缘效应，但问题仍存在。

代码实现 #

爆发生成函数（当前increments=100） #

from numpy import random
import numpy as np

lambd = 1/ts 
mu = lambd/increments
p1 = 1 - np.exp(-mu)

def burst_function(timemax, p1, increments):
    bursts = np.zeros(len(timemax)) # 创建爆发计数数组
    for j in range(len(timemax)):
        # 按泊松分布生成每个时间步的爆发次数
        unittime = random.choice([1,0], size=increments, p=[p1,1-p1])
        bursts[j] = np.sum(unittime)
    return bursts

指数脉冲光变曲线生成 #

先生成目标时长3倍的光变曲线，再截取中间1/3用于Structure Function计算：

from math import exp

s = 0
result = []
for i in range(len(timemax)):
    if bursts[i] > 0: # 有爆发发生时更新亮度
        s = bursts[i] + (s * exp(-1/td))
        result.append(s)
    else: # 无爆发时按衰减常数更新
        s = s * exp(-1/td)
        result.append(s)

Structure Function计算代码 #

from tqdm import tqdm
import numpy as np

def structure_function(timelength, lag, noisedata):
    SF = []
    noisedata = np.array(noisedata, dtype=np.float64)
    for m in tqdm(range(len(lag))):
        v_values = []
        k = int(lag[m])
        maximum = int(len(timelength)-k)
        for j in range(maximum):
            init = (noisedata[j+k] - noisedata[j])**2
            v_values.append(init)
        SF.append(np.sqrt(2 * np.mean(v_values)))
    return SF

异常现象对比 #

• 异常案例：当td=ts=20时，双对数图中第5、10个时滞点明显偏离平滑曲线
• 正常案例：当td=ts=25时，双对数曲线平滑无异常

已尝试的排查操作 #

1. 修正边缘效应：截取光变曲线中间1/3数据用于计算，排除边界影响
2. 提高模拟分辨率：增大increments或缩小时间步长，结果异常更严重

求助需求 #

排查异常产生的原因，并提供可行的解决方法。

内容的提问来源于stack exchange，提问作者agn-aud