等离子体诊断方法-第二章：折射

本文相关代码与介绍可在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公式：

μ^{2} = 1 - \frac{X (1 - X)}{1 - X - \frac{1}{2} Y^{2} \sin^{2} θ \pm [(\frac{1}{2} Y^{2} \sin^{2} θ)^{2} + (1 - X)^{2} Y^{2} \cos^{2} θ]^{1 / 2}}

其中，

X \equiv ω_{p e}^{2} / ω^{2}, Y \equiv ω_{c e} / ω, μ \equiv k c / ω

在托卡马克诊断时，往往存在背景磁场，且波矢k与磁场B0相互垂直。对于寻常波入射，存在E//B0，此时Y=0, 这时，折射系数可以写为：

μ^{2} = 1 - X = 1 - \frac{ω_{p e}^{2}}{ω^{2}} = 1 - \frac{n_{e}}{n_{c}}

光束经过等离子体产生的相移变换表示为：

ϕ_{p} = \int (k_{p l a s m a} - k_{0}) d l

折射系数同样满足定义：

μ \equiv \frac{k c}{ω} = \frac{k}{k_{0}}

因此，相移可以进一步的化简为：

\begin{aligned} ϕ_{p} = \int (k_{p l a s m a} - k_{0}) d l = - \int \frac{ω}{c} (1 - μ) d l \\ = - \frac{ω}{c} \int (1 - (1 - \frac{n_{e}}{n_{c}})^{1 / 2}) d l \end{aligned}

考察(1-n_e/n_c)^{1/2}，当等离子体中电子密度ne远远小于截止密度nc时，将此表达式应用泰勒展开有:

(1 - \frac{n_{e}}{n_{c}})^{1 / 2} = 1 - \frac{1}{2} \frac{n_{e}}{n_{c}} - \frac{1}{8} (\frac{n_{e}}{n_{c}})^{2} + o (. . .)

截取到第二项，相移表达进一步简化为：

\begin{aligned} ϕ_{p} = - \frac{ω}{c} \int (1 - (1 - \frac{n_{e}}{n_{c}})^{1 / 2}) d l = - \frac{ω}{2 n_{c} c} \int n_{e} d l \end{aligned}

截止密度n_c表示为(即对应截止频率时等离子体的振荡频率)：

n_{c} = \frac{ε_{0} m_{e}}{e^{2}} ω^{2} = \frac{ε_{0} m_{e}}{e^{2}} f^{2} [H z] = \frac{ε_{0} m_{e} c^{2}}{e^{2}} λ^{- 2} [m]

移表达进一步简化为：

\begin{aligned} ϕ_{p} = - \frac{ω}{2 n_{c} c} \int n_{e} d l = \frac{e^{2}}{π c ε_{0} m_{e}} f^{- 1} \int n_{e} d l \end{aligned}

即得到了相位移与等离子体密度弦积分的关系。

1.2 探测器信号的采集与处理

1.2.1 探测器采样信号的构成

某干涉仪组成原理如下所示，是一个典型的外差干涉仪：

在经过等离子体诊断区域的Beam1信号E1和LO信号E2可以表示为：

\begin{aligned} E_{1} = A_{1} \cos (ω_{1} t + ϕ (t)) \\ E_{2} = A_{2} \cos (ω_{2} t) \end{aligned}

未经过等离子体诊断区域的Beam1信号E0和LO信号E2可以表示为：

\begin{aligned} E_{0} = A_{0} \cos (ω_{1} t) \\ E_{2} = A_{2} \cos (ω_{2} t) \end{aligned}

在探测器上分别进行混频，有

\begin{aligned} s i g n a l_p l a s m a = (E_{1} + E_{2})^{2} & = A_{1}^{2} \cos^{2} (ω_{1} t + ϕ) + A_{2}^{2} \cos^{2} (ω_{2} t) \\ + A_{1} A_{2} \cos [(ω_{1} + ω_{2}) t + ϕ] + A_{1} A_{2} \cos [(ω_{1} - ω_{2}) t + ϕ] \\ = \frac{1}{2} (A_{1}^{2} + A_{2}^{2}) + \frac{1}{2} A_{1}^{2} \cos (2 ω_{1} t + 2 ϕ) + \frac{1}{2} A_{2}^{2} \cos (2 ω_{2} t) \\ + A_{1} A_{2} \cos [(ω_{1} + ω_{2}) t + ϕ] + A_{1} A_{2} \cos [(ω_{1} - ω_{2}) t + ϕ] \end{aligned}

\begin{aligned} s i g n a l_r e f e r = (E_{0} + E_{2})^{2} & = A_{0}^{2} \cos^{2} (ω_{1} t) + A_{2}^{2} \cos^{2} (ω_{2} t) \\ + A_{0} A_{2} \cos [(ω_{1} + ω_{2}) t] + A_{0} A_{2} \cos [(ω_{1} - ω_{2}) t] \\ = \frac{1}{2} (A_{1}^{2} + A_{2}^{2}) + \frac{1}{2} A_{1}^{2} \cos (2 ω_{1} t) + \frac{1}{2} A_{2}^{2} \cos (2 ω_{2} t) \\ + A_{1} A_{2} \cos [(ω_{1} + ω_{2}) t] + A_{1} A_{2} \cos [(ω_{1} - ω_{2}) t] \end{aligned}

对于此两束信号，可取出中频角速度为ω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 处理流程

整体流程为
整体流程图.png

鉴相器组成为
鉴相器组成.png

1.2.2.2 离散傅里叶变换

对于实信号x(t)的采样信号x[n]，通过离散傅里叶变换DFT可以变换为X[t]：

X [k] = \sum_{n = 0}^{N - 1} x [n] e^{- j (2 π / N) k n}, (0 \leq k \leq N - 1)

对应的离散傅里叶逆变换为：

x [n] = \frac{1}{N} \sum_{k = 0}^{N - 1} X [k] e^{j (2 π / N) k n}, (0 \leq n \leq N - 1)

需要注意的是：

变换后序列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，由于缺乏对实际物理图像的了解，导致中频会产生一点的偏差，在充分的了解频移特性与范围前，尚且无法很好的确定带通滤波器的范围。在此采用了一种伪带通滤波：先找到真实的中频，再依据一定展宽提取。
中频提取.png
在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 过零追踪鉴相器

当对比信号和实际信号的中频被提取出来后，即可通过追踪在每个采样时间内正向通过零点的时间获得相位差。两信号正向过零的时间为：

\begin{aligned} ω_{m} t_{R} = 2 π m_{R} + \frac{3}{2} π \\ ω_{m} t_{S} + ϕ_{p} = 2 π m_{S} + \frac{3}{2} π \end{aligned}

相移可以表示为：

ϕ_{p} = ω_{m} (t_{R} - t_{S}) + 2 π (m_{S} - m_{R})

在粗略的认为相移不会超过2π的前提下，相移可以简写为：

ϕ_{p} = ω_{m} (t_{R} - t_{S})

正向过零点追踪能够被简写为：

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 测试集测试

为了确定上述算法的可用性，这里通过一个示例：对于模拟数据

f (t) = \cos (2 π \cdot 10 \cdot t + \sin (t)) + \cos (2 π \cdot 30 \cdot t + s i n (t)) + \cos (2 π \cdot 1000 \cdot t) + n o i s e (t)

对应参考频率信号为：

f (t) = \cos (2 π \cdot 30 \cdot t)

采样频率为20Khz, 信号持续时间30秒，目标中频频率为30Hz。
经过上述流程，可以提取出中频信号与过零点：
过零点.png
并可以正确的提取出中频信号的相位差sin(t)：
相位差.png
验证了算法的正确性。

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

\begin{aligned} s i g n a l_p l a s m a = (E_{1} + E_{2})^{2} & = A_{1}^{2} \cos^{2} (ω_{1} t + ϕ) + A_{2}^{2} \cos^{2} (ω_{2} t) \\ + A_{1} A_{2} \cos [(ω_{1} + ω_{2}) t + ϕ] + A_{1} A_{2} \cos [(ω_{1} - ω_{2}) t + ϕ] \\ = \frac{1}{2} (A_{1}^{2} + A_{2}^{2}) + \frac{1}{2} A_{1}^{2} \cos (2 ω_{1} t + 2 ϕ) + \frac{1}{2} A_{2}^{2} \cos (2 ω_{2} t) \\ + A_{1} A_{2} \cos [(ω_{1} + ω_{2}) t + ϕ] + A_{1} A_{2} \cos [(ω_{1} - ω_{2}) t + ϕ] \end{aligned}

\begin{aligned} s i g n a l_r e f e r = (E_{0} + E_{2})^{2} & = A_{0}^{2} \cos^{2} (ω_{1} t) + A_{2}^{2} \cos^{2} (ω_{2} t) \\ + A_{0} A_{2} \cos [(ω_{1} + ω_{2}) t] + A_{0} A_{2} \cos [(ω_{1} - ω_{2}) t] \\ = \frac{1}{2} (A_{1}^{2} + A_{2}^{2}) + \frac{1}{2} A_{1}^{2} \cos (2 ω_{1} t) + \frac{1}{2} A_{2}^{2} \cos (2 ω_{2} t) \\ + A_{1} A_{2} \cos [(ω_{1} + ω_{2}) t] + A_{1} A_{2} \cos [(ω_{1} - ω_{2}) t] \end{aligned}

零差探测（Homodyne Detection）

当ωm ≪ ω1,ω2时，即ωm = 0, 对应的探测技术即为零差探测技术。
经过等离子体的诊断信号对应有表达式：

\begin{aligned} s i g n a l_p l a s m a = (E_{1} + E_{2})^{2} & = \frac{1}{2} (A_{1}^{2} + A_{2}^{2}) + \frac{1}{2} A_{1}^{2} \cos (2 ω_{1} t + 2 ϕ) + \frac{1}{2} A_{2}^{2} \cos (2 ω_{2} t) \\ + A_{1} A_{2} \cos [(ω_{1} + ω_{2}) t + ϕ] + A_{1} A_{2} \cos (ϕ) \end{aligned}

该信号的强度I主要贡献项为:

I_{p l a s m a} \propto \frac{1}{2} (A_{1}^{2} + 2 A_{1} A_{2} \cos (ϕ) + A_{2}^{2})

对于对比信号

\begin{aligned} s i g n a l_r e f e r = (E_{0} + E_{2})^{2} & = \frac{1}{2} (A_{1}^{2} + A_{2}^{2}) + \frac{1}{2} A_{1}^{2} \cos (2 ω_{1} t) + \frac{1}{2} A_{2}^{2} \cos (2 ω_{2} t) \\ + A_{1} A_{2} \cos [(ω_{1} + ω_{2}) t] + A_{1} A_{2} \cos (0) \end{aligned}

该信号的强度I主要贡献项为:

I_{r e f e r} \propto \frac{1}{2} (A_{1}^{2} + 2 A_{1} A_{2} + A_{2}^{2}) = \frac{1}{2} (A_{1} + A_{2})^{2}

主要能够通过接受信号的强度判断出判断出相位变换。但是由于强度项具有平方特征，因此只能得到相位变化的绝对值，无法通过强度判断相位是正向变换还是反向变换，也就无法直接判断出等离子体弦密度积分是增长还是下降。

外差探测（Heterodyne Detection）
当ωm≠0, 对应的探测技术即为外差探测技术。当ω1,ω2远远大于采样频率，而ωm$相比于ω1,ω2较小，且相对于采样频率而言能够满足Nyquist采样定理（采样频率大于fm两倍），此时就能够将混杂在信号中频率为fm的分量采集出来。即得到：

\begin{aligned} s i g n a l_p l a s m a & = A_{1} A_{2} \cos [(ω_{1} - ω_{2}) t + ϕ] \\ s i g n a l_r e f e r e n & = A_{1} A_{2} \cos [(ω_{1} - ω_{2}) t] \end{aligned}

此时可以通过比较两个信号相位的相位差直接获得相位变化，且能够反应相位的相对变化，可以分辨出相位信号是正向的还是逆向的。即能够判断出等离子体弦密度积分是增长还是下降的。同时，外差探测不需要严格的保证两个源频率相同，只需要保证源本身输出频率是稳定的。同样的，如果只考虑相对相位的变化，通过合理的布置可以使其具有抗环境干扰能力，起到类似于差分补偿电路的作用。

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

圆(1).png

光线在介质中传播常见的方程形式为：

\nabla n = \frac{d}{d s} (n \frac{d r}{d s})

以圆柱圆截面圆心为原点建立极坐标系，由圆的的几何特性易得，射线方向与法线的夹角始终等于射线方向与径向方向的夹角。由射线折射的斯涅尔定律：

n_{1} \sin (θ_{1}) = n_{2} \sin (θ_{2})

在球对称坐标系中，根据传播方程可以推出矢积关系：

r \times \frac{d r}{d η} = c o n s t

其中，η为路径L的弧长S ，该方程说明了光线始终在这个平面内传播，标量形式满足：

N (r) \cdot r \cdot \sin (α) = c o n s t

其中，α为射线方向与径向的夹角。在光线入射边界处，N(R)=1, Rsin(α) = y ，则常数即为:

N (r) \cdot r \cdot s i n (α) = y

根据光线运动的几何关系，可以得到

(\frac{r \cdot d θ}{d r})^{2} = \tan^{2} (α)

与折射定律联立，可以得到

(\frac{r \cdot d θ}{d r})^{2} = \frac{1}{N^{2} (r) \cdot r^{2} / y^{2} - 1}

对偏转角进行积分, 其中 r0 是射线上距离原点最近的

θ (y) = π - 2 \int_{r_{0}}^{\infty} | \frac{d θ}{d r} | d r = π - 2 \int_{r_{0}}^{\infty} (N^{2} (r) r^{4} y^{- 2} - r^{2})^{- 1 / 2} d r

折射率与等离子体密度直接具有关系式：

N (r) = (1 - \frac{n_{e}}{n_{c}})^{1 / 2}

当ne = n0(1-r^2/a^2)时，有

N (r) = \frac{r}{a}

因此θ(y)的表达式进一步写为

\begin{aligned} θ (y) = \sin^{- 1} \frac{V_{0} {[(\frac{y}{a})^{2} - (\frac{y}{a})^{4}]}^{1 / 2}}{{[V_{0} (\frac{y}{a})^{2} + (\frac{1 - V_{0}}{2})^{2}]}^{1 / 2}} \end{aligned}

其中，V0 = n0/nc 。对上式求解，可以得到对于θ最大的值为V0

等离子体诊断方法-第二章：折射

等离子体诊断方法-第二章：折射

写在前面基本用法

目录

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#232091

Shot#232092

Shot#232094

1.4 误差分析与遗留问题

2. Homodyne & Heterodyne

3. 推导柱对称中的散射

等离子体诊断方法-最终章：实验仿星器装置

等离子体诊断方法-第四章：等离子体辐射

数值计算三——微分方程组的求解与Modelica

传热学二：利用稀疏矩阵快速求解

评论 (1)

等离子体诊断方法-第二章：折射

等离子体诊断方法-第二章：折射

写在前面 基本用法

目录

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#232091

Shot#232092

Shot#232094

1.4 误差分析与遗留问题

2. Homodyne & Heterodyne

3. 推导柱对称中的散射

等离子体诊断方法-最终章：实验仿星器装置

等离子体诊断方法-第四章：等离子体辐射

数值计算三——微分方程组的求解与Modelica

传热学二：利用稀疏矩阵快速求解

评论 (1)

写在前面基本用法