等离子体诊断方法-第二章:折射
本文相关代码与介绍可在PlasmaDiagnostic/Interferometer at main · DMCXE/PlasmaDiagnostic (github.com)中找到
写在前面 基本用法
计算文件为DataProcess.py
, 存在两个主要功能:
SignalAnalsis 构成了一个数据处理模块,最大程度模拟了我想象中可以单纯凭借电路实现的功能。输入采样频率fs,目标中频fm, 诊断信号,对照信号,可返回相位差等
class SignalAnalsis:
def __init__(self,fs,fm_target,signal_plasma,signal_ref):...
density_string_integral构成了相位信息至密度信息的换算模块
def density_string_integral(phi,f=650*1e9):
'''
计算弦密度积分
'''
e = 1.60217662e-19
m_e = 9.10938356e-31
c = 299792458
epsilon_0 = 8.854187817e-12
return phi*np.pi*c*epsilon_0*m_e*f/(e**2)
以Shot#232094中sig1为例:
from DataProcess import SignalAnalysis, density_string_integral
SA = SignalAnalysis(
fs = fs,
fm_target = fm_target,
signal_plasma = sig1,
signal_ref = ref
)
phi = SA.phi
phi_fix = SA.phi #经过零点漂移矫正后的相位
ndl = density_string_integral
目录
等离子体诊断方法-第二章:折射写在前面 基本用法目录1 多道偏振干涉仪诊断系统的应用1.1 弦密度积分推导1.2 探测器信号的采集与处理1.2.1 探测器采样信号的构成1.2.2 相位的提取1.2.2.1 处理流程1.2.2.2 离散傅里叶变换1.2.2.3 中频提取1.2.2.4 离散傅里叶逆变换1.2.2.5 过零追踪鉴相器1.2.2.6 零漂补偿1.2.3 测试集测试1.3 实验数据的处理结果Shot#232091Shot#232092Shot#2320941.4 误差分析与遗留问题2. Homodyne & Heterodyne3. 推导柱对称中的散射
1 多道偏振干涉仪诊断系统的应用
1.1 弦密度积分推导
要在冷等离子体假设下的折射系数的Appleton-Hartree公式:
其中,
在托卡马克诊断时,往往存在背景磁场,且波矢k与磁场B0相互垂直。对于寻常波入射,存在E//B0,此时Y=0, 这时,折射系数可以写为:
光束经过等离子体产生的相移变换表示为:
折射系数同样满足定义:
因此,相移可以进一步的化简为:
考察(1-n_e/n_c)^{1/2},当等离子体中电子密度ne远远小于截止密度nc时,将此表达式应用泰勒展开有:
截取到第二项,相移表达进一步简化为:
截止密度n_c表示为(即对应截止频率时等离子体的振荡频率):
移表达进一步简化为:
即得到了相位移与等离子体密度弦积分的关系。
1.2 探测器信号的采集与处理
1.2.1 探测器采样信号的构成
某干涉仪组成原理如下所示,是一个典型的外差干涉仪:
在经过等离子体诊断区域的Beam1信号E1和LO信号E2可以表示为:
未经过等离子体诊断区域的Beam1信号E0和LO信号E2可以表示为:
在探测器上分别进行混频,有
对于此两束信号,可取出中频角速度为ωm = ω1-ω2。在实际的诊断设备中,LO源与Beam源的频率可达到数百GHz,而两者的差频(中频)往往是数MHz,对探测器信号采样频率往往为数十MHz。
由Nyquist采样定理:
为了不失真地恢复模拟信号,采样频率应该不小于模拟信号频谱中最高频率的两倍。
因此在该混频信号中,可达到数百GHz的LO源与Beam源的频率在数十MHz的采样频率下,完全无法满足Nyquist定理,因而对应的信号全部失真。即原始混频信号中角速度为ω1,ω2,ω1+ω2的的信号将失真,而信号差频(中频)对应的角速度分量ωm=ω1-ω2由于只有数MHz,相对采样信号而言满足Nyquist定理,能够被不失真的被采集出来。
1.2.2 相位的提取
1.2.2.1 处理流程
整体流程为
鉴相器组成为
1.2.2.2 离散傅里叶变换
对于实信号x(t)的采样信号x[n],通过离散傅里叶变换DFT可以变换为X[t]:
对应的离散傅里叶逆变换为:
需要注意的是:
- 变换后序列X[k]是复数序列,前一半数据与后一半数据共轭对称
- 变换后序列X[k]中k点的含义是:f = k fs/N,即序列索引对应采样频率的1/N
这也说明,如果需要通过离散傅里叶变换对频谱进行处理,最好在处理时同时考虑共轭频率部分对原始信号的贡献。
在python中, 这一步可以通过使用np.fft
包中的fft()
函数进行快速傅里叶变换(FFT)
fft_signla = np.fft.fft(signal)
1.2.2.3 中频提取
对原始信号进行带通滤波,在中频fm左右设置一个带通滤波器,提取中频fm范围内的主要频率,过滤掉高频与基频范围内的噪音。
在对于本数据的处理中,由于受到多普勒效应等或源信号频率误差,会导致实际中频并不严格的处于2.125Mhz,由于缺乏对实际物理图像的了解,导致中频会产生一点的偏差,在充分的了解频移特性与范围前,尚且无法很好的确定带通滤波器的范围。在此采用了一种伪带通滤波:先找到真实的中频,再依据一定展宽提取。
在python中提供了多种便于计算处理DFT的工具:
根据np.fft.fftfreq()
可以方便地获得每个索引对应的频率,由此即可计算归一化功率谱
def get_power_spectrum(fs,fft_signal):
#通过fft获得结果以及频率分布
frequencies = np.fft.fftfreq(len(fft_signal), 1/fs)
#计算不同频率分量的功率谱密度
power_spectrum = np.abs(fft_signal)**2
#计算总功率
total_power = np.sum(power_spectrum)
#计算每个频率的占比
power_percentages = (power_spectrum/total_power)*100
#归一化
power_percentages = power_percentages/np.max(power_percentages)
return frequencies,power_percentages
在归一化功率谱中,先选出占有主要功率强度的频率,通过与目标中频距离确定真实的中频位置:
#确定主要功率强度
def get_main_frequencys():...
return main_frequencies
#获得真实中频
def get_fm_real(fm_target,freq_main):...
return fm_real
构建带通滤波器,获得中频分量
freq_range = range * fm_real
return filtered_fft
1.2.2.4 离散傅里叶逆变换
对于滤波后产生的频谱,进行离散傅里叶逆变换,即可获得中频信号
fm_signal = np.fft.ifft(filtered_fft)
1.2.2.5 过零追踪鉴相器
当对比信号和实际信号的中频被提取出来后,即可通过追踪在每个采样时间内正向通过零点的时间获得相位差。两信号正向过零的时间为:
相移可以表示为:
在粗略的认为相移不会超过2π的前提下,相移可以简写为:
正向过零点追踪能够被简写为:
def zero_detector(fs,signal):
time_array = np.arange(0,self.time_tot,1/fs)
#如果认为信号的相位差不会大于2pi,那么可以简化上述过程
zero_crossing = np.where(np.diff(np.sign(signal))>0)[0]
t = np.array([time_array[tt] for tt in zero_crossing])
return t
追踪出零点,即可给出在每个采样时间对应的相位差:
def phase_detector(fs,fm,signal_plasma,signal_ref):
t_p = self.zero_detector(fs,signal_plasma)
t_r = self.zero_detector(fs,signal_ref)
return (t_r-t_p)*2*np.pi*fm
1.2.2.6 零漂补偿
由于诊断系统在空间中具有不对称性,信号采集器自身也可能存在温度部分不均等问题,导致诊断系统整体产生零点漂移,表现为在系统中为进行放电时,仍然能够测出相位差。因此需要对零点漂移进行修正。这里非常简单的考虑将相位差整体平移消除直流分量。
def phase_zero_fix(self,phi):
#对相位做零漂补偿
phi_0 = phi[:int(0.1*len(phi))]
phi_average = np.sum(phi_0)/len(phi_0)
return phi - phi_average
1.2.3 测试集测试
为了确定上述算法的可用性,这里通过一个示例:对于模拟数据
对应参考频率信号为:
采样频率为20Khz, 信号持续时间30秒,目标中频频率为30Hz。
经过上述流程,可以提取出中频信号与过零点:
并可以正确的提取出中频信号的相位差sin(t):
验证了算法的正确性。
1.3 实验数据的处理结果
实验数据共有三组,分别为Shot#232091、Shot#232092和Shot#232094,每个数据文件中包含的信号种类如下图示
数据处理时,首先指定文件(以sig1为例):
import scipy.io as sio
shot = '232094'
data = sio.loadmat(shot+'.mat')
Ip =data['ip1'].reshape(1,-1)[0]
t_Ip = data['t_ip1'].reshape(1,-1)[0]
sig1 = data['sig1'].reshape(1,-1)[0]
指定采集频率、目标中频
fs = 60 *1e6
fm_target = 2.125*1e6
将数据输入鉴相器
SA = SignalAnalysis(
fs = fs,
fm_target = fm_target,
signal_plasma = sig1,
signal_ref = ref
)
将数据输入弦密度积分计算器
phi = SA.phi_fix #获取零点漂移修正后的相位差数据
dsi = density_string_integral(phi) #计算出弦密度积分
最后,能够得出Shot#232091、Shot#232091和Shot#232091的计算数据
Shot#232091
Shot#232092
Shot#232094
1.4 误差分析与遗留问题
- 在消除零点漂移的时候,采用的方法过于简单,不能产生很好的效果
- 实际上,注意到三炮数据在电流下降段,会出现多个密度峰值,在这些峰值处,实际对应的相位差大于2π,实际上无法使用之前我们假设的相位差小于2π的计算假设。但是如果归到[0,2π]中,这些峰值又会变成瞬时的第密度甚至零密度。如何解释这些数据需要进一步讨论。
- 本文数据处理的前提是:尽量模仿能够在FPGA或片上系统实现的信号处理方式。存在一些更先进的算法,例如Hilbert变换、EMD经验模态分解、小波变换等,具有很大的分析上述数据与解释的空间和准确性。为了更好的处理这些数据、获得更丰富的物理结果,希望能够在未来的工作中进一步发展。
2. Homodyne & Heterodyne
无论是零差干涉仪还是外差干涉仪,两臂信号混合之后均会产生两臂频率相加项和两臂频率相减相,主要关心频率相减项ω_m = ω1 -ω2
零差探测(Homodyne Detection)
当ωm ≪ ω1,ω2时,即ωm = 0, 对应的探测技术即为零差探测技术。
经过等离子体的诊断信号对应有表达式:
该信号的强度I主要贡献项为:
对于对比信号
该信号的强度I主要贡献项为:
主要能够通过接受信号的强度判断出判断出相位变换。但是由于强度项具有平方特征,因此只能得到相位变化的绝对值,无法通过强度判断相位是正向变换还是反向变换,也就无法直接判断出等离子体弦密度积分是增长还是下降。
外差探测(Heterodyne Detection)
当ωm≠0, 对应的探测技术即为外差探测技术。当ω1,ω2远远大于采样频率,而ωm$相比于ω1,ω2较小,且相对于采样频率而言能够满足Nyquist采样定理(采样频率大于fm两倍),此时就能够将混杂在信号中频率为fm的分量采集出来。即得到:
此时可以通过比较两个信号相位的相位差直接获得相位变化,且能够反应相位的相对变化,可以分辨出相位信号是正向的还是逆向的。即能够判断出等离子体弦密度积分是增长还是下降的。同时,外差探测不需要严格的保证两个源频率相同,只需要保证源本身输出频率是稳定的。同样的,如果只考虑相对相位的变化,通过合理的布置可以使其具有抗环境干扰能力,起到类似于差分补偿电路的作用。
3. 推导柱对称中的散射
#存在大量问题,还需要进一步修改
Consider a beam propagating along a chord of a refractive cylinder at a distance y from the axis. Suppose the cylinder has a refractive index N(r), where r is the radius. Obtain a general equation for the angular deviation of the beam θ due to refraction when θr≪y so that the chord can be approximated as straight. In the case where ω>>ωp and n = no(1-r^2/a^2), calculate the value of y at which 0 is greatest and prove that this maximum 0 is n0/nc
光线在介质中传播常见的方程形式为:
以圆柱圆截面圆心为原点建立极坐标系,由圆的的几何特性易得,射线方向与法线的夹角始终等于射线方向与径向方向的夹角。由射线折射的斯涅尔定律:
在球对称坐标系中,根据传播方程可以推出矢积关系:
其中,η为路径L的弧长S ,该方程说明了光线始终在这个平面内传播,标量形式满足:
其中,α为射线方向与径向的夹角。在光线入射边界处,N(R)=1, Rsin(α) = y ,则常数即为:
根据光线运动的几何关系,可以得到
与折射定律联立,可以得到
对偏转角进行积分, 其中 r0 是射线上距离原点最近的
折射率与等离子体密度直接具有关系式:
当ne = n0(1-r^2/a^2)时,有
因此θ(y)的表达式进一步写为
其中,V0 = n0/nc 。对上式求解,可以得到对于θ最大的值为V0
作者的编程能力令人佩服。柱散射问题可以只利用斯涅耳定律,对相关物理量做泰勒展开,取至一阶近似,可以得到散射角的线性微分关系,积分,进行数学处理即可,最后只需手算y*sqrt(a^2-y^2)的极值即可。请问圆(1).png这幅图作者是用什么软件画的?流程图是如何画的?请作者指教。