Astropy版本：5.0.1

astropy.modeling的文档：Models and Fitting (astropy.modeling) — Astropy v5.0.1

astropy提供了astropy.modeling模块用于对数据进行特定模型下的拟合，笔者结合自身的使用经验，在本文中以二维高斯模型与moffat模型拟合为例子介绍该模块的各项功能。
此外, astropy.modeling还提供了线性模型、三角函数模型、反三角函数模型等常见函数模型。

1. 基本功能介绍——以二维高斯拟合为例

已知二维高斯函数的定义为：
$$
f(x, y) = A e^{-a\left(x - x_{0}\right)^{2} -b\left(x - x_{0}\right) \left(y - y_{0}\right) -c\left(y - y_{0}\right)^{2}}
$$
其中$a,b,c$为与$x$轴、$y$轴方向上的方差及函数半长轴指向有关的常数。
简单的拟合程序示例如下：

1
2
3
4
5
6
7
8
9
import numpy as np
from astropy.modeling import models,fitting

#data为拟合数据
data = data
model_1 = models.Gaussian2D(amplitude=60,x_mean=170,y_mean=250,x_stddev=1, y_stddev=1)
fit_1 = fitting.LevMarLSQFitter()
x, y=np.mgrid[:400, :400]
data1 = fit_1(model_1, x, y, data)#模型输出的数据

models.Gaussian2D()内的各项参数为模型的初始值，对应含义如下：
amplitude：高斯模型的初始峰值，即函数中的$A$。
x_mean与y_mean：高斯模型中心点的初始位置，即函数中的$x_0$与$y_0$。
x_stddev与y_stddev：高斯模型的初始方差，即函数中的$\sigma_x$与$\sigma_y$。

fitting.LevMarLSQFitter()为选用的拟合器（Fitter），注意不同的模型需要使用对应的拟合器。

限制参数

考虑到该模型有六个自由参数，在拟合上有一定难度，可通过限制参数减少自由度：

1
2
#使x_stddev与y_stddev保持初始值不变。
gaussian = models.Gaussian2D(amplitude= 80 , x_mean= 8.5, y_mean=9, x_stddev=1, y_stddev=1，fixed={'x_stddev': True,'y_stddev':True})

或者限制参数的变化范围：

1
2
3
4
5
def tie_stddev(model):
	return model.y_stddev_0
#使x_stddev与y_stddev的变化范围限制在0.1到1.5之间
gaussian = models.Gaussian2D(amplitude= 180 , x_mean= 8.5, y_mean=9, x_stddev=1.3, y_stddev=1.3,bounds={'x_stddev': (0.1,1.5),'y_stddev':(0.1,1.5)})
gaussian.x_stddev.tied = tie_stddev

多个模型叠加

在实际处理过程中，有时需要多个模型叠加对数据进行拟合。如用二维高斯模型拟合一个天体源时，考虑到背景的影响，应在此基础上再加一个常数模型拟合。astropy中模型叠加的操作很简单：

1
2
3
c = models.Const2D(amplitude=10)
gaussian = models.Gaussian2D(amplitude= 180 , x_mean= 8.5, y_mean=9, x_stddev=1.3, y_stddev=1.3,)
model = gaussian + c

获得拟合后的模型参数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
#获取高斯函数参数
>>>model_1.amplitude
Parameter('amplitude', value=200.0)

>>>model_1.x_mean
Parameter('x_mean', value=8.5)

>>>model_1.x_mean.value
8.5

#获取高斯函数的fwhm
>>>model_1.x_fwhm
2.3548200450309493

2. PSF拟合——Moffat函数拟合

在进行测光时，我们需要对星点进行拟合以获得其半高全宽，一般可采用二维高斯函数、二维moffat函数或一维moffat函数进行拟合。
若使用二维moffat函数，代码示例如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
#x_cen与y_cen为星点在数据中的对应位置
#peak为星点在数据中的流量最大值

#选定拟合的范围
 analysis_radius = 10

#计算出该圆形对应的各坐标点
 xmin = x_cen - analysis_radius
 xmax = x_cen + analysis_radius
 ymin = y_cen - analysis_radius
 ymax = y_cen + analysis_radius

 x, y=np.mgrid[:20, :20]

model = models.Moffat2D(amplitude=1.0,x_0=0,y_0=0,gamma=2,alpha=3.5)
fitter = fitting.SimplexLSQFitter()
fitted_model = fitter(model, x,y,data[ymin:ymax,xmin:xmax] / peak)

考虑到有时星点的成像质量不佳，此时采用一维moffat函数拟合可取得更好的效果，其具体思路也是以目标星点为中心，划定一定大小的圆，按像素点与圆心的距离，将二维数据重新组合为以距离为变量的一维函数：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
#x_cen与y_cen为星点在数据中的对应位置
#peak为星点在数据中的流量最大值

flux_counts = []
pixel_distance = []

#选定拟合的范围
analysis_radius = 10

xmin = x_cen - analysis_radius
xmax = x_cen + analysis_radius
ymin = y_cen - analysis_radius
ymax = y_cen + analysis_radius

x, y=np.mgrid[:20, :20]

for x in range(x_cen - analysis_radius, x_cen + analysis_radius):
	for y in range(y_cen - analysis_radius, y_cen + analysis_radius):
		flux_counts.append(((data[y][x]) / peak))
		pixel_distance.append(np.linalg.norm((x_cen - x, y_cen - y)))

model = models.Moffat1D(amplitude=1.0, x_0=0, gamma=2., alpha=3.5)
fitter = fitting.SimplexLSQFitter()
fitted_model = fitter(model, pixel_distance, flux_counts)