星系测光：理论基础与实操

十万行出头天

已于 2023-07-22 11:39:50 修改

阅读量1.3k

点赞数 3

分类专栏： 2023年CSST星系科学研究生暑期学校-学习记录文章标签：天文星系测光 Python photutils

于 2023-07-02 21:01:58 首次发布

本文链接：https://blog.csdn.net/qq_41959720/article/details/131505270

版权

2023年CSST星系科学研究生暑期学校-学习记录专栏收录该内容

2 篇文章 4 订阅

订阅专栏

在望远镜相机拍摄到星系之后，我们便可以对其开展大量的科学研究。做为星系科学研究的一个重要手段，星系测光便显得尤为重要

星系测光的内容与目的

在星系测光中，我们试图从图像里获得各星系的：

本征参数：星系的物理性质。
表征参量：星系的几何形态或其他性质。在测量表征参量时，我们即便不知道红移等信息也可以进行。此次学习主要围绕的便是表征测量。

如何从图像中探测到星系

通常，在望远镜拍摄的一幅图像中，会出现不止一个星系。尤其是对于现在的大视场望远镜而言，在每一幅图像中所出现的星系数量，都是我们肉眼所看不过来的。这个时候，就需要用程序来帮助我们去探测这些星系了。

星系探测的大致算法如下：设置一个阈值（例如面亮度、天光背景起伏），在图像中超出阈值的连续像素的数目较多时，则认为这些连续的像素组成一个星系。通常在选取阈值的时候，会让这个值比较地小，以确保连续性。

常用星系参数及测量

星系位置

星系位置的确定方法：

将最亮像素的所在位置，作为星系的位置：不适用于形状不规则的星系。
求出图像的质心（一阶矩），即对流量取加权平均值。现在一般用的便是这个方法。

面亮度

星系的面亮度 $I$ ，表示单位立体角接收到的流量，常用单位为 $L_\odot/pc^2$ 、 $mag/arcsec^2$ ，这个面亮度与观测的距离无关。

星等

孔径星等

用圆形或椭圆形的孔径，对星系进行流量测量。

自动孔径星等

由二阶矩的测量结果，来决定用于测光的椭率和方位角。

总星等

在孔径测光的孔径足够大时，测得星等趋向的值。

标准星等

定义B波段面亮度为25 $mag/arcsec^2$ （有时候也用26.5）的等照度线是星系的边缘，用等照度线之内的像素流量积分，计算出来的便是标准星等。

模型星等

用模型拟合星系，从而计算得到的星等。

Petrosian星等

由定义的Petrosian半径进行计算。

形状

对图像求二阶矩，可以得到长轴、短轴以及方位角等信息。

形状参数的应用：

根据椭率和方位角，研究星系棒。
与物理性质的相关性。

如果求出图像的三阶矩，则可以得到有关图像扭曲的信息。

半光半径

流量达到一半时，所对应的孔径测光半径。

聚度指数

定义为：

$CI = m_{aper1} - m_{aper2}$

即用小的孔径星等减去大的孔径星等。

利用mag(total) v.s. 聚度指数（例如，mag(3 pixels) - mag(8 pixels)），则可以区分点源和展源。

颜色

星系的颜色为这个星系在不同波段上星等的差。在测量星系颜色时，要求在不同波段用相同孔径进行测光。

形态

星系的形态可以分为早型和晚型：

早型（椭圆星系）：外观平滑，没有结构特征，颜色偏红，金属丰度高。
晚型（漩涡/棒旋星系）：颜色较蓝。
从早型到晚型的规律：恒星越来越年轻。

非模型分类

中心聚集度

按百分位流量划分的大圈与小圈半径比值，值越大表示中心聚集度越大。早型星系的聚集度更高。

非对称性

原始图像减去旋转180度后的图像。

基尼系数

表征各个像素流量分布的均匀性，早型星系的基尼指数更大，晚型更小。

表面亮度轮廓

对于椭圆星系而言：遵循 $R^{1/4}$ 律，更准确地，Sersic law： $R^{1/n}$ 。在观测上，HST高分辨率图像发现并不能用简单的Sersic law描述椭圆星系成分。相应的模型改进：core+sersic.

对于旋涡星系而言：核球成分（与椭圆星系类似）+盘成分。在观测上，HST+2MASS+DSS H band的观测拟合结果，揭示出了中心有star cluster的成分。拟合发现盘星系存在不同的类型。

等照度线

表面亮度相同的等值线，通常可以用椭圆来描述。

由等照度线，可以定义星系大小：

B波段25 / 26.5 $mag\ arcsec^{-2}$
Petrosian $R_{50}$
模型拟合得到的半光度/半质量半径。

星系测光的主要步骤

扣除天光背景

天光背景来源于大气辉光、黄道光、热辐射、月光等一系列现象。

要扣除天光背景，首先就要进行天光背景的计算：

一般使用众数（可以通过均值和中位值来估算）。
将图像切成网格，计算网格的天光背景，做平滑处理。
网格取太小时，会出现过亮的情况，应该用周围网格的中位值替代。
网格大小的经验选取方法：10 PSF FWHM.
利用星系的外围区域计算（外围区域不受星系影响），取高斯分布的均值。天光背景有起伏。
利用天体附近空白的区域，计算平均强度。
不能使用常数表示天光背景，应该有一个二维的天光背景图像。
- 划分格子，用中值滤波进行平滑。
- 将空白的区域拟合成一个二维斜面。

背景测量分为全局背景测量和局域背景测量：

全局背景测量：同前，盒子的选择，求均值并取平滑。
局域背景测量：用于目标源的测光和拟合。

如果大背景没有减好，则还需要对局域背景进行修正。

评估天光背景测量的质量：

空白区域强度值应该分布在0附近，并且没有大尺度的空间变化。

扣除周围恒星影响

星系周围的恒星，会污染我们所测量到的流量，因此需要在测光之前就将这些恒星的影响给去除掉。去除的方式为使用mask进行掩膜。

为了扣除明亮的污染源，可以使用中值滤波的方式来进行：对图像做一个中值滤波，构建平滑图像，以去除掉其中的噪声，从而得到一个大背景。

等照度线拟合

第一次等照度线拟合，目的：获取中心的位置。

第二次等照度线拟合，固定中心位置，目的：得到星系盘的性质（倾角、椭率）。

最后，固定所有参数，提取不同半长轴处的亮度曲线，从而计算得到半光半径、聚度等信息。

据此，可以进行棒结构的辨认：椭率和PA会在棒的半径附近有一个转折。但值得注意的是，仅强棒有很明显的特征。

对图像中的星系进行批量测光

前面提到的方法，主要是针对单个星系的情形。通常我们对望远镜拍摄到的图像，是会直接用一个pipeline，来对从里面探测到的星系进行测光的。这样的pipeline，通常会对背景进行一个测量，去除掉明亮污染源的影响，然后设置阈值，并由此在图像中进行源的抓取。

一般情况下阈值设置得会比较小，因此这样一来选出的源会有很多是混在一起的，这个时候就需要进行一个名为“deblending”的操作了。deblending的关键参数DEBLEND_MINCONT，对结果的影响非常大。为了确保漩涡星系不被拆分，在对漩涡星系deblending时，用比较大的值；对密集星系团场而言，用比较小的值，确保小的源能被区分开。

在完成源探测与deblending后，便是对这些星系（或点源）的测光操作。测光的方式有很多种（例如ISO、AUTO、ISOCORR、APERTURE等），用不同的方式测出的结果都有所差异，具体选择哪一种来做为结果，需要看自己的需求。

测光完成后，结合图像的星等零点，便可以得到星系的源表，里面包含了大量的相关星系参数。这样一来，我们便从一幅图像中提取出了大量星系并对其进行了批量测光。

常用星系测光工具：photutils

photutils是Python的一个用于星系测光的模块，可以通过pip install的方式安装。

值得注意的是，photutils在版本更新时，会有部分函数的用法被更改，因此不同版本的用法可能会有所不同。这里使用的photutils的版本为1.8.0，Python版本为3.11.3.

在此案例中，背景已被事先扣除，因此将会从mask的制作开始做起。下图为完成背景扣除后的UGC 9476（i波段）图像：

mask制作——扣除亮源污染

我们可以看到在这幅图像中，除了中心的星系以外，还有非常多的点源散布在图像各处。这些点源的存在会对后续的测光带来一定的影响，因此，我们事先要做的，就是将这些点源给“遮住”，即前面说的mask.

首先打开fits图像文件：

from astropy.io import fits

# 打开图像
path_img = 'UGC9476_i.fits'
img = fits.getdata(path_img)
img = img.byteswap().newbyteorder()

随后，设置一个阈值，在图像中进行源探测，并将探测到的源转换为mask，最后计算得到亮源的标准差：

from astropy.stats import sigma_clipped_stats, SigmaClip
from photutils.segmentation import detect_threshold, detect_sources

# 设置阈值
sigma_clip = SigmaClip(sigma=3.0, maxiters=10)
threshold = detect_threshold(img, nsigma=2.0, sigma_clip=sigma_clip)

# 根据阈值，在图像中进行源探测
segment_img = detect_sources(img, threshold, npixels=3)

# 将探测到的源作成mask
premask = segment_img.make_source_mask(size=5)

# 统计mask的均值、中值与标准差
mean, median, std = sigma_clipped_stats(img, sigma=3.0, mask=premask)
print(mean, median, std)

根据刚才得到的标准差的值，再次对图像进行正式的源探测，并从得到的mask中将星系所对应的那一部分移除掉：

import sep
import copy

# 从图像中找出各源，seg_map_all为各源对应的mask
objects, seg_map_all = sep.extract(img, thresh=2.0, err=std, minarea = 5, 
                           deblend_nthresh=300, deblend_cont=0.1, segmentation_map=True)

# 星系的中心坐标
x0 = 255
y0 = 259

# 将星系所对应的那一块从mask中去除掉
def seg_remove_obj(seg, x, y):
    seg_copy = copy.deepcopy(seg)   
    seg_copy[seg == seg[int(y), int(x)]] = 0  # 移除所有具备与(x,y)一样的值的pixel
    return seg_copy

seg_map = seg_remove_obj(seg_map_all, x0, y0)

# 将mask转化为bool类型
seg_map_True = seg_map[:,:] != 0

最后，对得到的mask进行卷积（应该是为了起到一个类似于PSF的作用），并筛选出大于一定阈值的像素，得到最终的mask过后的图像：

from astropy.convolution import Gaussian2DKernel,convolve
import numpy as np

# 对mask进行卷积，并筛选出大于一定阈值的像素
conv = convolve(seg_map_True.astype(int), Gaussian2DKernel(2)) 
seg_mask = (conv >= 0.1)
mmask = seg_mask.astype(int)

# 最终得到mask过后的图像
img_masked = np.ma.array(img, mask=mmask)

将mask前后的图像展示出来：

from astropy.visualization.mpl_normalize import ImageNormalize
import matplotlib.pyplot as plt
from astropy.visualization import LogStretch

# mask前的图像
plt.figure(figsize=(10, 5))
plt.subplot(121)
norm = ImageNormalize(stretch=LogStretch(a=2000), vmin=0, vmax=6)
plt.imshow(img, origin='lower', cmap='Greys_r', norm=norm)
plt.xticks([])
plt.yticks([])

# mask后的图像
plt.subplot(122)
norm = ImageNormalize(stretch=LogStretch(a=2000), vmin=0, vmax=6)
plt.imshow(img_masked, origin='lower', cmap='Greys_r', norm=norm)
plt.xticks([])
plt.yticks([]);

并将mask保存为fits图像文件：

# 保存mask后的图像
fits.writeto("UGC9476_i_mask.fits", mmask, overwrite=True)

从生成的mask图像中，便可以发现星系周围的亮源基本上都被遮挡了。

第一次等照度线拟合——确定星系中心位置

准备工作

读入科学图像文件：

from astropy.io import fits

path_img = 'UGC9476_i.fits'

# 大致估计的星系中心坐标
ra = 220.383
dec = 44.513

# 读取科学图像文件
hdu = fits.open(path_img)
img = hdu[0].data
head = hdu[0].header

将估算的星系中心坐标，由赤经、赤纬转化为像素坐标：

from astropy.wcs import WCS

# 将pixelSize的单位由degree转为arcsec
pixelSize = abs(head['CD2_2']) * 3600

# 将估计的星系中心坐标，根据header中的wcs信息，转化为像素坐标
wcs = WCS(head)                          
objXY = wcs.all_world2pix([ra], [dec], 0)

objXp = objXY[0][0]  
objYp = objXY[1][0] 
posi=[objXp, objYp]
print(posi)

读取mask图像，并对科学图像进行mask操作：

import numpy as np

# 读取mask图像，并对科学图像进行mask操作
path_mask = 'UGC9476_i_mask.fits'
mask = fits.getdata(path_mask)

img_masked = np.ma.array(img, mask=mask)

开始拟合

在拟合前，先设置将要用于拟合的必要初始参数：

# 准备开始第一次等照度线拟合，先进行初始参数的设置
params = {}

# 椭率
params['e0'] = 0.1

# 椭圆的方位角（photutils中的方位角定义是从x轴开始算起的）
params['pa0'] = np.pi / 4

# 椭圆的半长轴尺寸（单位：像素），最好设置成略小于星系的尺寸
params['sma0'] = 90

# 估算的星系中心位置
params['xcen0'] = posi[0]
params['ycen0'] = posi[1]

# 像角元尺寸，即一个pixel对应多大的张角
params['pixscl'] = pixelSize

随后开始进行等照度线的拟合：

import time
from photutils.isophote import EllipseGeometry
from photutils.isophote import Ellipse

time_start = time.time()

geometry = EllipseGeometry(params['xcen0'], params['ycen0'], params['sma0'], params['e0'], params['pa0'])

ellipse = Ellipse(img_masked, geometry)

# 进行等照度线拟合
# 参数含义：不固定中心，不固定方位角，不固定椭率。拟合时，半长轴从1逐步增加到180像素，步长为0.15
iso_free = ellipse.fit_image(fix_center=False, fix_pa=False, fix_eps=False, 
    minsma=1, maxsma=180, step=0.15, maxgerr=1)

time_end = time.time()
print('totally cost: ',int(time_end-time_start))

在等照度线拟合期间，会使用不同大小的椭圆去拟合亮度轮廓，从而得到各椭圆参数随半长轴大小的变化。将这样的关系画出来：

import matplotlib.pyplot as plt

# 将拟合失败的点去除掉
ig = iso_free.stop_code != 4

# 创建总图像
plt.figure(figsize=(10, 7.5))
plt.subplots_adjust(top=0.92, bottom=0.08, left=0.10, right=0.95, hspace=0.4, wspace=0.35)

# 椭率随半长轴的变化
plt.subplot(321)
plt.errorbar(iso_free.sma[ig]*params['pixscl'], iso_free.eps[ig], yerr=iso_free.ellip_err[ig], 
             fmt='o', markersize=3, color='black')
plt.xlabel('SMA (arcsec)')
plt.ylabel('Ellipticity')
plt.ylim(0, 0.4)

# 方位角随半长轴的变化
plt.subplot(322)
plt.errorbar(iso_free.sma[ig]*params['pixscl'], iso_free.pa[ig]/np.pi*180.-90, yerr=iso_free.pa_err[ig]/np.pi*180., 
             fmt='o', markersize=3, color='black')
plt.xlabel('SMA (arcsec)')
plt.ylabel('PA (deg)')

# 中心x坐标随半长轴的变化
plt.subplot(323)
plt.errorbar(iso_free.sma[ig]*params['pixscl'], iso_free.x0[ig], yerr=iso_free.x0_err[ig], fmt='o', 
             markersize=3, color='black')
plt.xlabel('SMA (arcsec)')
plt.ylabel('xcen')

# 中心y坐标随半长轴的变化
plt.subplot(324)
plt.errorbar(iso_free.sma[ig]*params['pixscl'], iso_free.y0[ig], yerr=iso_free.y0_err[ig], fmt='o', 
             markersize=4, color='black')
plt.xlabel('SMA (arcsec)')
plt.ylabel('ycen')

# 面亮度轮廓
plt.subplot(325)
kk = iso_free.intens > 0
mu = -2.5*np.log10(iso_free.intens[kk]/params['pixscl']**2) + 22.5
mu_err = 2.5*iso_free.int_err[kk]/iso_free.intens[kk]/np.log(10)
plt.errorbar(iso_free.sma[kk]*params['pixscl'], mu, yerr=mu_err, color='black', fmt='o', markersize=3)
plt.ylim(30, 16)
plt.ylabel('$\mu$ mag/arcsec$^2$')
plt.xlabel('SMA (arcsec)')

# 生长曲线
plt.subplot(326)
plt.scatter(iso_free.sma*params['pixscl'], iso_free.tflux_e/np.max(iso_free.tflux_e), color='black', marker='.')
plt.ylabel('fraction')
plt.xlabel('SMA (arcsec)')
plt.ylim(0, 1.1)

# 保存图片
plt.savefig('Example/iso_step1.png', bbox_inches='tight')

可以看到下述图像：

左上：椭率
右上：方位角
左中：x轴中心位置
右中：y轴中心位置
左下：面亮度
右下：生长曲线（该椭圆中包含了占比多少的流量）

可以发现，x、y坐标的中心位置，在半长轴较小的时候并没有发生太大的变化，而随着半长轴的增大，这两个值开始发生了较为剧烈的变化。由于星系在中心的结构相对地还比较对称，因此，对于这两个值，我们可以取中心区域的平均值，作为这个星系的中心。注意到这两个值在SMA < 5的时候还算比较稳定，因此，我们对SMA < 5的区域取平均，从而得到星系中心的位置：

# 找出SMA < 5的index
kk = iso_free.sma*params['pixscl'] < 5

# 对SMA < 5的区域取平均，得到星系中心的像素位置并相应地更新参数
xcen = np.mean(iso_free.x0[kk])
ycen = np.mean(iso_free.y0[kk])
params.update({'xcen':xcen, 'ycen':ycen})
print('The galaxy center: ', xcen, ycen)

画出第一次等照度线拟合的示意图：

from astropy.visualization.mpl_normalize import ImageNormalize
from astropy.visualization import PercentileInterval,AsinhStretch
import matplotlib.patches as patches

# 新建画布
plt.figure(figsize=(6, 6))
ax1=plt.gca()

# 设置图像归一化与对比度
norm1 = ImageNormalize(img, interval=PercentileInterval(99.7), stretch=AsinhStretch())

# 显示星系图像
plt.imshow(img, origin='lower', cmap='Greys_r', norm=norm1)
plt.xticks([])
plt.yticks([])
plt.title('UGC 9476; i band image; free parameters')

# 画出在各半长轴长度下拟合的等照度线
Niso = len(iso_free.sma)
for k in range(0, Niso):
    ellipse = patches.Ellipse((iso_free.x0[k], iso_free.y0[k]), 2*iso_free.sma[k], 
                      (1-iso_free.eps[k])*2*iso_free.sma[k], 
                      angle=iso_free.pa[k]*180/np.pi, alpha=0.5, fill=False, color='r')
    ax1.add_artist(ellipse)

# 保存图片
plt.savefig('Example/img_iso_step1.png', bbox_inches='tight')

十万行出头天

关注

3
点赞
踩
15

收藏

觉得还不错? 一键收藏
0
评论
星系测光：理论基础与实操

星系探测的大致算法如下：设置一个阈值（例如面亮度、天光背景起伏），在图像中超出阈值的连续像素的数目较多时，则认为这些连续的像素组成一个星系。尤其是对于现在的大视场望远镜而言，在每一幅图像中所出现的星系数量，都是我们肉眼所看不过来的。星系周围的恒星，会污染我们所测量到的流量，因此需要在测光之前就将这些恒星的影响给去除掉。（有时候也用26.5）的等照度线是星系的边缘，用等照度线之内的像素流量积分，计算出来的便是标准星等。早型星系的聚集度更高。表征各个像素流量分布的均匀性，早型星系的基尼指数更大，晚型更小。
复制链接

扫一扫