一种用于等离子体密度测量的自适应相位差计算方法

：
1.本发明涉及等离子体密度测量技术领域，具体涉及一种用于等离子体密度测量的自适应相位差计算方法。


背景技术：

2.目前受控核聚变实现的约束方式主要有惯性约束和磁约束，托卡马克是实现商用聚变堆的磁约束核聚变装置之一。托卡马克等离子体物理不仅是一门非常复杂的学科，而且随着等离子体参数的不断提高、新的运行模式的出现，有越来越多的聚变问题值得去研究和探索。通常等离子体密度参数通过相位的方式进行测量，由相位信息和转换公式得到电子密度的空间分布。所以提高等离子体密度测量的相位差检测技术能极大的完善对等离子体探究。
3.等离子体的是包含正负离子的中性电离气体，目前采用hcn(hydrogen cyanide)激光干涉仪对等离子体测量，其中hcn的参考光路与探针光路分别通过真空和等离子体，两路光信号因为介质不同而产生相位差，再由tgs(triglyceride sulfate)探测器将光信号转换为电信号。
4.在托卡马克装置中应用广泛的是硬件相位差计，计算机使用采集卡采集硬件相位差计转换后的相位差电压，但是该方法无法滤除直流干扰，测量分辨率低，并且只能适用于单频率测量。国内外研究使用傅里叶变换相位差比较法进行测量，该技术是通过快速傅里叶变换(fast fourier transform，fft)的方法来计算相位。该方法可以在抑制噪声和谐波干扰的同时获得高分辨率的相位，然而傅里叶变换中的栅极效应会影响相位测量精度，或者非周期采样带来频移，导致测量的结果误差较大。针对fft的不足，采用的解决方法有先对等离子体序列做fft和ifft得出相位，该方法可以不需要直接信号的频率，具有高时间分辨率和数值分辨率。另外，采用ap-fft(all phase fast fourier transform)，该方法抑制频率泄露能力强，并且具有相位不变性。
5.另一方面，目前托卡马克等离子体测量时为单频率载波测量，但不同介质对不同波长的吸收率不是相同，托克马克等离子体测量也会进行多频率多方向的测量。而目前对多频率相位测量的进行研究还在起步，一般已经设定信号采样率和fft转换的点数。工程上，信号的频率偏移和频谱泄露导致频谱重叠，进而测量的相位不准确。目前只研究适应了电力谐波的相位测量，有效频率为倍数关系，具有一定的特殊性，尚无完全的适用方法解决多频率相位测量频谱重叠方法。


技术实现要素：

6.本发明所要解决的技术问题在于提供一种用于等离子体密度测量的自适应相位差计算方法，解决信号处理中fft、dft只能人为提前设定参数的问题，并且减少托克马克中多频率测量、栅极效应的影响导致的误差不准确的问题。
7.本发明所要解决的技术问题采用以下的技术方案来实现：
8.一种用于等离子体密度测量的自适应相位差计算方法，包括以下步骤：
9.以下步骤s1-s7是以探针信号xd(n)的多频率相位计算为例，参考信号xr(n)的相位计算过程相同，只需将步骤中的xd(n)换成xr(n)即可。
10.s1计算探针信号xd(n)中的有效频率：对序列xd(n)进行2n点的fft，其中fs为信号的采样率，fft的点数由采样率决定。
[0011]2n-1
＜fs＜2nꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0012]
序列xd(n)的fft结果为yd(n)，yd(n)的实部与虚部取模结果为my(n)。取采样信号量程的1/2n为阈值，即归一化量程幅值的一半。如果序列my(n)中数值大于阈值，则存在有效频率，有效频率为估计频率fn(j)(j∈[0,n])。频率分辨率为fo＝fs/2n，频谱泄露会导致超过阈值的数值呈现连续性，导数为0则表示该频率是存在的，对于一段连续的频率，导数为零的频率点为有效的峰值频率。
[0013][0014]
本发明对估计频率的百位以下进行取零处理，一方面由于温度的影响，载波频率会产生频率偏移，而且载波频率会随着温度波动而产生微小变化；另一方面本发明方法的相位结果具有频率的抖动的抗干扰性。
[0015]
s2信号的滤波处理：将序列yd(n)以置零的方式将干扰频率滤除，n为序列的下标。
[0016][0017]
滤波结束采用ifft为新的序列数据xi(n)。
[0018]
s3根据有效频率确定序列dft的转换点数：df为相邻有效频率的最小间隔差，dft转换点数nd由最小间隔差、信号频率和采样率的最大公约数决定。
[0019]df
＝min(fn(2)-fn(1),fn(3)-fn(2),fn(j)-fn(j-1))
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0020]
nd＝fs/min(fn(j),fs,df)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0021]
s4根据dft的点数自适应dft混合基的基数：假设存在nd＝n1×
n2…×
nr(r为有限项)，则n点的dft拆分为短序列nr(r∈[0,n])。nd点从小到大依次取质数表中1000以内大于1的质数进行整除，nr为质数的幂指数，p为质数。
[0022][0023]
s5子序列的预处理：每次以步进(2
×nd-1)的长度截取xi(n)为子序列xs(ns)(ns＝2
×nd-1，s∈[1,n])，每一个子序列xs(ns)对应一个相位结果。取子序列xs(ns)中点位置为xs(n0)，xs(ns)从左到右以步进1依次截取nd长度的子序列x
ss
(nd)(ss∈[1,n/nd])，共nd组x
ss
(nd)。nd组x
ss
(nd)都包括了xs(n0)数据，每组x
ss
(nd)首尾相连循环右移，达到设置xs(n0)为序列的第一位。将nd组旋转序列x
ss
(nd)对应下标相加，并对结果乘以1/nd得到预处理序列xa(ns)。
[0024]
s6计算序列xa(ns)的dft结果ya(ns)：s4已经确定n1，n2…
nr的值，根据式(7)将信号xi(n)排列为n1行n2列的数组，对每一行进行dft，结果为y1。y1乘以旋转因子得
最后对y2的每一列进行dft，得到ya(ns)。
[0025][0026]
其中，为旋转因子，k0，k1为输出的序列的下标，式中k0，k1，m0，m1由下面的式(8)决定。
[0027][0028]
s7计算有效频率的相位：根据fn(j)的结果，依次从ya(ns)计算出多个频率的相位，为有效频率的相位。
[0029][0030]
s8计算相位差：由s1-s7的计算，得到探针信号xd(n)的多频率相位参考信号xr(n)的相位两序列数据对应相减，得到相位差pd(j)。探针信号与参考信号各频率的相位进行相差得到探针信号因为透过等离子体而发生的相移。从相位差信号来看，初始情况等离子体未放电，相位差为零，但是实际上存在偏差数值，该相位定义为p0(j)，所以这个“零漂”需要去掉，选取初始状态前t时间内的相位取平均值为零漂相位。一般情况下，t取10-20ms中的值。对应等离子体密度的相位需要减去初始偏差相位p0(j)。
[0031][0032][0033]
p(j)＝pd(j)-p0(j)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)
[0034]
s9相位差叠加：由于反正切的计算相位限制在-π/2～π/2内，而等离子体放电过程实际相位差超过此范围，则根据信号变化率和相位差变化阈值判断信号发生翻转是正翻转、负翻转或者干扰。对应等离子体密度的相位输出进行叠加计算，从而得到对应真实的相位。
[0035][0036]
其中，k为翻转次数。当相位差为正翻转则k为正数，相位差为负翻转则k为负数。
[0037]
本发明的有益效果是：
[0038]
1)本发明采用离线数据的方式来计算等离子体密度，不需要硬件平台，采用云计算的方式实现等离子体密度的测量。本发明以偏移的实际频率、有效频率成份和相邻频率差作为dft转换点数的因数，进行自适应的点数调节，能更有效的减少频谱泄露和栅极效应求取有效相位的误差，提高相位准确度。
[0039]
2)本发明针对实际上混合基dft分裂需要提前设定的情况，根据互质方式自适应选择确定混合基的dft的分裂项。
[0040]
3)本发明对傅里叶变换的数据进行截取、翻转、平移、叠加和归一化操作以减少数
据的频谱泄露影响。
附图说明：
[0041]
图1为本发明计算方法的流程示意图；
[0042]
图2为托卡马克等离子体密度数据采集过程的结构图；
[0043]
图3为离散hcn中探针信号数据的频谱图；
[0044]
图4为子序列xs(ns)预处理过程；
[0045]
图5为离散hcn中两路信号相位数据；
[0046]
图6为离散hcn信号数据的计算结果图；
[0047]
图7为相位去除初始零漂后的相位差图；
[0048]
图8为离散hcn信号数据的相位叠加结果图。
具体实施方式：
[0049]
为了使本发明实现的技术手段、创作特征、达成目的与功效易于明白了解，下面结合具体实施例和图示，进一步阐述本发明。
[0050]
符号说明
[0051]
探针信号：携带等离子体信息的光路
[0052]
参考信号：不携带等离子体信号的初始光路
[0053]
xd(n)：探针离散数据中等离子体对应的序列
[0054]
xr(n)：参考离散数据中等离子体对应的序列
[0055]
n：自然数
[0056]
j：为序列中有效频率的下标
[0057]fs
：信号的采样率
[0058]
fft：快速傅里叶变换(fast fourier transform)
[0059]
ap-fft：全相位傅里叶变换(all phase fast fourier transform)
[0060]
ifft：傅里叶逆变换(inverse fast fourier transform)
[0061]
dft：离散傅里叶变换(discrete fourier transform)
[0062]
yd(n)：探针(参考)等离子体放电序列fft的结果
[0063]my
(n)：yd(n)取模运算的结果
[0064]fn
(j)：有效频率的估计频率
[0065]df
：相邻的有效频率的最小间隔差
[0066]
n1×
n2…×
nr：混合基拆分后的整数
[0067]
xs(ns)(ns＝2
×nd-1,s∈[1,n/nd])：以步进(2
×nd-1)的长度截取序列xd(n)的子序列
[0068]
xs(n0)：xs(ns)中点位置
[0069]
x
ss
(nd)：以步进1，宽度nd截取xs(ns)的子序列
[0070]
xa(ns)：预处理最终结果序列
[0071]
有效频率的相位
[0072]
t：系统开始后初始阶段的一段时间
[0073]
pd(j)：两路信号相位差结果序列
[0074]
p0(j)：零漂相位
[0075]
p(j)：去除零漂后的相位差
[0076]
p
p
(j)：等离子体密度放电过程中相位差真正的相位
[0077]
k：正翻转或者负翻转的次数
[0078]
本实施例以j-text装置1011685次放电试验的数据进行说明。
[0079]
如图2所示，等离子体电子密度是通过两路光信号进行相位差得到的结果，其中携带等离子体密度信息的光路为探针信号，未携带任何信号、做比较的光路为参考信号。两路光信号由hcn激光干涉仪发出。以下s1-s7过程为以探针信号xd(n)的多频率相位计算为例，参考信号xr(n)的相位计算过程相同，只需将步骤中的xd(n)换成xr(n)即可。
[0080]
s1计算探针信号xd(n)中的有效频率：对xd(n)序列进行2n点的fft，其中fs为信号的采样率，fft的点数由采样率决定。实施例中等离子体密度的信号采样频率为fs＝250khz，所以n为18。
[0081]2n-1
＜fs＜2n[0082]
序列xd(n)的fft结果为yd(n)。yd(n)的实部与虚部取模结果为my(n)。取采样信号量程的1/2n为阈值，即归一化量程幅值的一半。如果序列my(n)中数值大于阈值，则存在有效频率，有效频率为估计频率fn(j)(j∈[0,n])。频率分辨率为fo＝fs/2n，则分辨率精度小于1hz，并且在栅极效应影响下，有效频率不在频率分辨点上，进而发生频谱泄露情况。频谱泄露会导致超过阈值的数值呈现连续性，而偏导数为0则表示该频率是存在的，采用的比较的方式来确定有效信号的存在，对于一段连续的频率，导数为零的频率点为有效的峰值频率。
[0083][0084]
本发明对估计频率的百位以下进行取零处理，一方面由于温度的影响，载波频率会产生频率偏移，而且载波频率会随着温度波动而产生微小变化；另一方面本发明方法的相位结果具有频率的抖动的抗干扰性。图3是离散hcn中信号数据的频谱图，如图3.b所示，my(n)的频率在10.4khz附近,并且频谱成分复杂，10.3khz至10.5khz有多个有效频率，其中10.38khz为这段区域的峰值频率，所以有效频率估计为10.4khz。参考信号的频谱入3.a所示，其有效估计频率和针信号相同。
[0085]
s2信号的滤波处理：将序列yd(n)以置零的方式将干扰频率滤除。n为序列的下标。
[0086][0087]
滤波结束采用ifft为新的序列数据xi(n)。
[0088]
s3根据有效频率确定序列dft的转换点数：df为相邻有效频率的最小间隔差。dft转换点数nd由最小间隔差、信号频率和采样率的最大公约数决定。
[0089]df
＝min(fn(2)-fn(1),fn(3)-fn(2),fn(j)-fn(j-1))
[0090]
nd＝fs/min(fn(j),fs,df)
[0091]
由于目前载波的等离子体信息都是单频率信号，并且由图3可知信号只有10.4k的
有效频率，通过计算可得dft的转化点数约为24。
[0092]
s4根据dft的点数自适应dft混合基的基数：假设存在nd＝n1×
n2…×
nr(r为有限项)，则n点的dft拆分为短序列nr(r∈[0,n])。nd点从小到大依次取质数表中1000以内大于1的质数进行整除，nr为质数的幂指数，p为质数。
[0093][0094]
s5子序列的预处理：每次以步进(2
×nd-1)的长度截取xi(n)为子序列xs(ns)(ns＝2n
d-1,s∈[1,n])，每一个子序列xs(ns)对应一个相位结果。如图4所示，令xs(ns)中点位置为xs(n0)，xs(ns)从左到右以步进1分别截取nd长度的序列x
ss
(nd)(ss∈[1,n/nd])，共nd组x
ss
(nd)。nd组x
ss
(nd)都包括了xs(n0)数据，每组x
ss
(nd)首尾相连循环右移，达到设置xs(n0)为序列的第一位。将nd组旋转序列x
ss
(nd)对应下标相加，并对结果乘以1/nd得到预处理序列xa(ns)。
[0095]
s6计算序列xa(ns)的dft结果ya(ns)：s4已经确定n1，n2…
nr的值，根据式(7)将信号xi(n)排列为n1行n2列的数组，对每一行进行dft，结果为y1；y1乘以旋转因子得最后对y2的每一列进行dft，得到ya(ns)。
[0096][0097]
其中，为旋转因子，k0，k1为输出的序列的下标，式中k0，k1，m0，m1由下面的式子决定。
[0098]
k＝n2k0+k1(k0∈[0,n
1-1],k1∈[0,n
2-1])
[0099]
n＝n1m1+m0(m0∈[0,n
1-1],m1∈[0,n
2-1])
[0100]
s7计算有效频率的相位：根据fn(j)的结果，依次从ya(ns)计算出多个频率的相位，为有效频率的相位。
[0101][0102]
探针信号和参考信号两路求的有效频率的相位如图5所示，其中a为参考信号的相位，b为探针信号的相位。
[0103]
s8计算相位差：由s1-s7的计算，得到探针信号xd(n)的多频率相位参考信号xr(n)的相位两序列数据对应相减，得到相位差pd(n)。探针信号与参考信号各频率的相位进行相差得到探针信号因为透过等离子体而发生的相移。从相位差信号来看，初始情况等离子体未放电，相位差为零，但是实际上存在偏差数值，该相位定义为p0(j)，所以这个“零漂”需要去掉，选取初始状态前t时间内的相位取平均值为零漂相位，对应相位需要减去初始偏差相位p0(j)。在本实施例中，t取20ms。
[0104]
[0105][0106]
p(j)＝pd(j)-p0(j)
[0107]
探针信号与参考信号有效频率的相位差结果如图6所示，相位开始存在的零漂现象，初始相位不为零。通过去除初始零漂操作后，其结果如图7所示，相位差大小随着等离子体放电而相位持续增加并翻转。当增加到最大后停止上升并保持一段时间，随后相位开始向下减少并翻转，最后回归至初始相位。
[0108]
s9相位差叠加：由于反正切的计算相位限制在-π/2～π/2内，而等离子体放电过程实际相位差超过此范围，则根据信号变化率和相位差变化阈值判断信号发生翻转是正翻转、负翻转或者干扰。对应等离子体密度的相位输出进行叠加计算，从而得到对应真实的相位，其中k为翻转次数。当相位差为正翻转则k为正数，相位差为负翻转则k为负数。
[0109][0110]
对计算的相位结果经过去零、翻转步骤得到等离子体放电时的相位信息，如图7所示，其中在100ms时进行放电，约800ms放电结束，相位最高达到14.8948rad。该过程符合等离子体放电过程，其方法有效。
[0111]
以上显示和描述了本发明的基本原理和主要特征和本发明的优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。