大数据相加_FAST数据处理流程[凑数]

FAST相关的内容在这个专辑里，再发布一次凑数

https://mp.weixin.qq.com/mp/appmsgalbum?action=getalbum&album_id=1479243467993202689&__biz=MzI1NjQyNzI0OQ==#wechat_redirect

我又看了一遍感觉差不多就这么些个过程，简单的说就是把谱全都加起来，可能rfi的部分可以用 https://github.com/DuFanXin/RFI-Net 这里的程序标记一下

我在这个文章

https://www.aanda.org/articles/aa/full_html/2020/06/aa38483-20/aa38483-20.html

附录里也仔细描写了数据处理过程，等我得到2020B数据后可能会写一个tracking模式的数据处理流程，估计是差不多的内容，纯属凑数

20200819

今天这个笔记用来解释FAST数据处理，考虑到FAST现在还在风险共担的调试阶段，数据的格式和发布方式还可能会变，所以我这里没有给出来大段大段的程序，而是说说算法，方便自己以后还能回忆起一些基本的思路，希望能说明白：

FAST的数据是个catalog fits，如果是一秒一次读出，五分钟曝光后，每个beam的数据大约313M， 数据在 extenction 1，主要内容有：

OBSNUM          LONG64       1

SCAN            LONG64       2

OBSTYPE         STRING       'ON              '

QUALITY         STRING       'T'

UTOBS           DOUBLE        58620.412

DATE_OBS        STRING   '2019-05-17T09:52:54.000Z'

EXPOSURE        DOUBLE        1.0066330

NCHAN           LONG64        65536

FREQ            DOUBLE        1000.0036

CHAN_BW         DOUBLE        0.0076293945

DATA            FLOAT     Array[4, 65536]

这些，还有一些描述观测时的状态的一些关键词，比如风速等等都会陆陆续续加进来，这些值可以帮助我们去除一些噪音可能比较大的数据，对于一个简单的数据处理来说，目标就是看看大概观测到了什么，此时需要从文件里提取频率，流量信息和做叠加

提取频率会用到的参数有：

FREQ ：起始频率

NCHAN：频率的通道个数

CHAN_BW：每个频率通道的宽度

所以可以生成一个自然数列：0, 1, 2, ..., NCHAN-1，乘以CHAN_BW 加上FREQ，即是FAST扫描的频率范围数列，也就是：

freq = findgen(tab[0].NCHAN) * tab[0].CHAN_BW + tab[0].FREQ

DATA里存的是数据，包括四种偏振，每种偏振有65536个流量，对应于FAST扫描的各个频率，现在这个设置是1秒读出一个谱，噪声管的周期是两秒

这个catalog文件里包括299个这样的DATA文件，是5分钟的曝光得到的数据，用idl读入这个catalog后的指标结构是：

tab[五分钟得到的299条谱].data[4种偏振,65536个频率处的流量]

tab[i].data[j,k]表示第i次曝光的时候收到的第j偏振的谱里的第k个数，tab[i].data[j,*]表示第i次曝光的时候收到的第j偏振的谱的总共的65536个数，对应前面定义的freq里65536个频率bin，可以用

plot, freq, tab[0].data[1,*]

查看一下谱形

这里多解释一点点偏振，在19波束接收机的每个接收喇叭里有个十字排列的导线，相当于用接收机喇叭汇聚信号后，接收的天线，十字形状天线会收集到偏振沿着两个天线方向的信号，因此0号，1号偏振加在一起代表的是接受到的总的信号，根据可靠情报，对于我这种仅仅想测个星系HI流量的科技工作者来说，另外两个偏振可以无视，因此把0号，1号偏振方向收到的数据按照相同频率叠加起来，就是总共收集到的来自指向方向的流量了

于是这里有个鉴定信号和干扰的机会，星系的21cm线信号是没有偏振的，如果是干扰的电磁波可能会有偏振，也就是在一个方向的强度会大于另一个方向，当这两个偏振方向接受到的流量很不一样的时候，可能是个有偏振的信号的干扰

接下来需要把偏振0，偏振1的这299根谱叠加起来，如果就这么按照每个频率的值全部相加的话，几乎就看到最终接受到的 on the target 的流量了，之所以说几乎，是因为这直接相加的是不对的，数据还需要定标，也就是确定一下纵坐标的刻度和温度的关系，FAST直接收到的流量数据值有时候是10的15次方，需要找到这个数据值和温度的换算关系，重新调整前面画的那个谱的纵坐标刻度

定标的方法是用噪音管在接收的流量上间歇的加一个已知温度的噪音，我们先看一下噪音是怎么加进去的

假设现在噪音是1秒打开，1秒关闭，那么添加噪音的周期就是两秒，噪音对于FAST来说也是收集到的数据，并且FAST接收到的谱每隔一秒读出一次，那么如果我们把整个曝光时间内的某个频率处的流量画出来，会看到这个流量的值整体有个2秒的周期，流量的振幅就是我们加进去的噪音了，假设这个振幅高的地方读数是A，低的地方读数是B，我们加进去的噪音10K对应的刻度是A-B，因此单位纵坐标对应 10/(A-B) K

具体操作的时候可以这么着，先找到目标谱线附近的一段儿频率范围:

index_freq_noise = where(freq gt center_freq - 40 and freq lt center_freq - 20)

把这段儿频率范围的值简单平均一下，得到这段儿频率附近的值

noise_check = fltarr(299)

for i = 0, 298 do noise_check[i]= median(tab[i].data[0,index_freq_noise])

这就在noise_check里存了299个数，对应从第一秒曝光到最后一秒曝光的某个频率附近的测量值，如果画出来：

plot, noise_check, psy = -symcat(16), /ynoz

看到的是一个横坐标单位是曝光时间的时间序列，纵坐标是每一秒收到的特征流量，此时明显会看到有两秒的周期，然后可以这么区分一下高低：

ind = where(noise_check gt median(noise_check))

noise_on = alog10(median(noise_check[ind]))

这个noise_on就是前面说的A的对数了，第二秒读数的时候没噪音：

ind = where(noise_check lt median(noise_check))

noise_off = alog10(median(noise_check[ind]))

noise_off 相当于前面说的B的对数，我这里取对数再用10次方算回来是为了保持数字精度，这么噪音引起的FAST接收流量值的改变是 10^noise_on - 10^noise_off ，目前对于FAST的噪音管 high 模式来说是10K：

on_0_K = 10./(10^noise_on - 10^noise_off)

这是 on the target 的5分钟曝光文件的 第0偏振 的谱的纵坐标刻度值到温度值的换算系数，也就是前面说的10/(A-B)

所以把第0偏振的每条谱的纵坐标乘以 10/(A-B) 就把每秒读出来的谱纵坐标变成温度了

我们还要记得，得到的光谱里一半儿是加了 10K 噪音的，一半儿没加这个噪音，因此要记得把那些加过 10K 噪音的谱的纵坐标扣掉手工增加的 10K

这个用噪音校准的过程要分别对0号，1号偏振来做，各自叠加出谱，成为偏振0总谱和偏振1总谱，再把两个这偏振方向的总谱按照频率相同的bin各自加起来，这样就得到了一个5分钟on the target曝光文件里的，包含所有偏振信息的，谱了

需要注意的是，这里偏振0和偏振1的温度刻度系数是不一样的，所以要分开来做，以及建议把俩偏振谱的叠加谱留下来，各自画一下看看量级上差别多大，噪音结构是不是相似，也把俩偏振谱相减，那些起伏比较大的地方可能出现一个偏振方向强度过大的情况，有可能来自一些干扰信号，可以再对比另一次曝光的数据，鉴定是不是哪里不正常，这么说来，配置曝光时间的时候，尽量分成两个on，两个off，相互也好有个照应

对于off的文件也要做同样的事情，得到总得off的谱，相当于一个背景噪音

把on的总谱减掉off总谱，得到的几乎就是属于源和噪音的谱

这里好像还需要修正一下每个频率的响应，主要思路是：如果有一个平谱从天而降，用FAST收到后是不会真的得到一个平的谱的，各个频率bin里的流量相对来说会有个高低，这个效率也要修正一下，比如off的时候，FAST是对一个没有信号的地方观测噪音，此时对于一个略微窄一些的频率范围，比如做HI线关心的频率范围的谱来说，也许各个频率处的背景流量没什么大区别，可是FAST接收读出的off数据里的谱会显得不太平，就可以用来修正各个频率bin的效率，有点儿像CCD天文学的平场

这个谱只要换算成了温度，即使不是同一天的数据也可以叠加，如果对谱的信噪比不满意，可以继续补曝光时间

这就是FAST目前的数据的处理思路了，整体感觉还是挺直观的，下面是几个注记：

0，为什么温度可以相加：

对于射电波段来说，首先假设源的辐射是黑体辐射，其次射电波段对应的是黑体辐射的一个长波的尾巴，此时有1900年瑞利提出来的 瑞利-金斯定律：黑体辐射和温度正比，因此对于射电波段来说，温度相加也就是能量相加

1，驻波：

这么大个望远镜接受到的数据会有驻波，按照 五百米口径球面射电望远镜(FAST) 观测性能概述.pdf 的说法，大概是这么个感觉：

驻波示意图，摘自 五百米口径球面射电望远镜(FAST) 观测性能概述.pdf

图中的这个18.8K附近，一抖一抖的小震荡就是驻波了，其产生原理我至今还没听懂，模糊的印象是有个干扰信号，来来回回流窜于镜片和接收机之间，波长也被不断调整，终于各个频率都有一点儿，就有了驻波，参考市面上别的大的射电望远单镜，驻波的扣除方法可以是用正弦或者什么类似的函数拟合，或者干脆看点儿比较亮的源，远离驻波困扰

2，干扰源，比如卫星：

尽管FAST附近没有移动电话信号塔，并且有电磁静默规定，有一些卫星从天上路过的时候，会带来一些干扰，FAST园区有个监测卫星的天线，如果有干扰信号太强的卫星通过就赶紧调整FAST的工作状态，即使如此，平时仍然有各个频率的卫星过境，而且由于卫星在运动，所以FAST接受到的频率还有个宽度，并且在on的时候，off的时候频率可能有个缓慢的移动，这种卫星过境引起的干扰只在几分钟内会影响FAST的流量，如果曝光时间长一些的话，就会只能在几个on和off的时候看到这个信号，之后就消失了，这个特点可以用来判断是不是卫星信号的干扰

3，定标到Jy

FAST从温度到Jy的转换系数是16Jy/K，可以用望远镜性质估计出来，比如Nobeyama的望远镜官网给了大量的公式来暗示这个量，不像ESO一些望远镜直接在主页上提供这个数字，这个定标数字是可以通过观测已知流量的源来校准的，对于FAST来说现在就是16Jy/K，现在FAST每天都会观测标准源来做校准

4，最后的但是不是最不重要的

封面是半边莲的花，象征一下我对FAST半瓶子晃荡的理解水平