现在，取代污染环境并且只能利用几十年的陈旧选择，受控的 核聚变似乎是用于制造无污染的、几乎无限的能源的最优选择。
它的原理简单，与发生在太阳中的机制相同，以在被称作 Tokamak的机器中最终再生为目的。通过强磁场和几百万安培的强 电流的结合作用而被强约束的在非常高的温度处电离的气环(被称 作等离子体)在其中心进行生成中子的氚化氘聚变反应，该中子传 送能量：这是Tokamak的理论基础。
聚变装置的物理和技术约束的最优化引出了“先进Tokamak” 概念的定义。它与提供固定的改进约束模式一致，其中，是以非感 应方式生成电流总量(由自动生成的等离子体自举电流给定大部分 电流)。这种“先进Tokamak”模式的实现需要控制电流密度分布 的能力，电流密度分布仅可以借助于其他的非感应电流 (non-inductive current)生成方法来得到。在已知的关于Tokamak 的各种方法中，大功率电磁波的注入构成了用于在等离子体中生成 附加非感应电流的理想选择。由于这个原因，能够控制混合型波的 功率沉积分布是关键，这是当前用于生成附加电流的最有效方法。 借助于光谱测定HXR诊断，在Tokamak Tore Supra(TS)中研究这 种波的传播和吸收。测量由混合型波加速的SE在HXR范围内发射 的辐射是获得关于波功率沉积分布的信息的最有效方法。在1996 年1月，该事实判定在TS中具有极好的空间和时间分辨率的光谱 测定的新HXR诊断设备中是合适的，该设备特别适于研究在长周 期上的电流分布的控制(见参考文献[1])。光谱测定系统包括分布 在2个相机内的59个检测器，其中，一个相机水平设置，而另一 垂直设置，从而通过使等离子体的截面与斜率大不相同的检测器的 视线相交，使其能够增大测量的空间冗余。该诊断法测量沿着每条 视线积分的等离子体发射率，其主要目的在于根据原始的积分测量 结果来确定等离子体发射率的径向分布。这可以通过传统Abel反 演法在一定假设条件下来实现。
图1示出了Tore Supra真空管V以及HXR诊断的环形等离子 体CP和视线L(弦)的实例的截面图。图1所示的结构包括两个 相机1a、1b、以及环形等离子体CP位于其中的管V。相机1a是 具有2 1条视线L的垂直相机，而相机1b是具有38条视线L的水 平相机。该结构并不是唯一可能的结构。例如，这些相机可以交换 它们的位置，或者可以在剖面的框架中改变它们的相对位置，或者 也可以只使用一个相机。当由等离子体发射的光子撞击置于相机中 的检测器时，这些光子释放其全部或部分能量。在那时释放的能量 被转换为脉冲。检测器由镉碲(CdTe)半导体制成。是通过专门对 CdTe进行最优化的处理电子通路来实现对来自检测器的脉冲的处 理。图2示出了根据现有技术的这种处理电子通路。
处理电子通路包括相机1、接收结构2、偏振电路3、电源电路 4、校准电路5、处理电路6、以及数据存储单元7。开关8将接收 结构2的输出连接至处理电路6的输入(在测量阶段中)或者将其 连接至校准电路5的输入(在校准阶段中)。相机1包括基于CdTe 半导体的检测器9、前置放大器10、以及微分发射器(differential emitter)11。接收结构2包括微分接收机12和线性放大器13。例 如，偏振电路3利用等于-100V的偏振电压来使检测器偏振。电 源电路4将电能提供给相机1的电子电路10和11、以及接收结构 2的12和13，例如，使用+/-12V，40mA的电源。处理电路6 包括一组鉴别器D1～D8、一组计数器C1～C8、以及数据采集单元 14。
检测器9是等离子体发射的光子P将其所有或部分能量传递到 其中的材料介质(material medium)。在检测器中传递的能量被转换 为电脉冲，然后，通过专门对CdTe进行最优化的电子计数系统来 处理该脉冲。通过前置放大器10来收集半导体中的电荷载体。微 分发射器11通过微分接收机12将由前置放大器10输出的信号传 输至线性放大器13(通常被称作“整形器”)。整形器的作用是将所 接收的脉冲(所接收的脉冲通常具有相当长的张驰时间(relaxation time)因此如果计数速度过快则存在重叠的风险)转换为可以在采 集系统的最后部分中轻松地计数的相当短的脉冲。可以人工调节整 形器的增益，以校准能量域内的信号。
在测量阶段内，开关8将接收结构2的输出连接至处理电路6 的输入。接着，通过八个积分鉴别器D1～D8来分析所接收的脉冲 幅度。当脉冲上升沿的幅度大于鉴别阈值时，积分鉴别器D1～D8 将逻辑信号发送至与积分鉴别器D1～D8相连的计数器C1～C8。 通过计数器Ci(i＝1，2，...，8)接收逻辑信号，将1加到计数器 Ci的缓冲存储器，因此，该存储器中存有对能量大于鉴别阈值的脉 冲进行计数的数量。在每个采样步骤中(即，获得数据的时间，例 如，16ms)，读取每个计数器的缓冲存储器，然后通过将八个计数 结果传送数据存储单元7中的数据采集单元14将缓冲存储器重置 为0。
该系统具有几个缺点。
首先，与输入信号相关的信息是不可用的，这就阻碍了任一成 形脉冲的显示，从而不能分辨由于两个光子同时达到检测器而造成 的叠加。同样，测量信号不是实时可用的，这阻碍了实时进行任一 分布反演，从而阻碍混合型波功率沉积的任一反馈控制和电流分布 的任一反馈控制。
为了获得可靠的测量结果，校准步骤是必须的。然后，将接收 结构2的输出连接至校准结构5的输入。
校准包括调节整形器电路的增益，以使其在由接收结构2输出 的脉冲幅度与入射光子的能量之间具有良好的一致性。如上所述， 根据现有技术的光谱测定系统包括两个相机，一个是垂直的，而另 一个是水平的，垂直的相机包括21个检测器，而水平的相机包括 38个检测器，从而总共给出了59个检测器。随后，对每个检测器 进行校准。
校准是必须的，以获得不同能量通道中发射率分布的精确重 构。可以利用具有1024个通道的数字光谱仪并利用三个放射源来 进行校准。然后，调节整形器的增益，以将每个放射源的主峰值设 置在适当的能量处。
校准步骤也具有缺点。它需要将采集系统中的不包括在校准中 的那部分电子装置(electronics)断开。这样做的结果将导致校准误 差。此外，这种断开增加了对系统所做的操作，从而增加了系统劣 化的风险。此外，相机1远离与校准台相连的采集系统。则这意味 着，当操作员不得不更改源相对于相机的位置时，他必须往返多次。
借助于采集电子电路来消除上述缺点，这是以申请人的名义于 2004年2月24日向法国专利局提交的题为“Circuit électronique de diagnostic spectrométrie”的法国专利申请的主题，其申请号为04 50388。下文中，参照图3到图7描述该采集电子电路。
根据本发明的方法(process，过程)处理从采集电子电路(如 在上述法国专利申请中描述的采集电子电路)输出的数据。根据本 发明的过程的优点在于获得实时发射率分布(例如，在前述的16ms 内执行的处理)，从而能够监测和控制实时发射率分布。
根据本发明的方法包括一些计算步骤。应当注意，在这些计算 中，符号“*”用于表示任何种类的相乘。

本发明涉及一种用于利用至少包括一个检测器的光谱测定系 统来确定超热电子的局部发射率分布的过程，超热电子来自位于环 形管(toric vessel)中的被称作等离子体的电离气环，该检测器相 对于环形管设置以使检测器的视线与环形管的横截面及半径为“a” 的电离气环的横截面相交，电离气环的横截面的中心偏离于环形管 的横截面的中心，其特征在于，包括以下步骤：
-计算0阶贝塞尔函数Jo的前NB个零点Z1，Z2，...，ZNB，
-构建矩阵Jρ，其k行j列的元素为Jo(ρk*Zj)，其中，Jo(ρk*Zj) 是自变量(ρk*Zj)(k＝1，2，...，NM以及j＝1，2，...，NB)的 0阶函数Jo，其中，ρk是等离子体的点PK与等离子体中心之间相对 于半径“a”的归一化距离，矩阵Jρ如下：
Aρ＝Jρ*C，
其中，Aρ是数组，其元素代表沿归一化等离子体半径ρ的发射率分 布，以及C是矩阵的系数，
-读取测量数据，该测量数据包括表示沿视线通过检测器测量 的等离子体积分发射率的等离子体发射率数据Y，以及包括大半径 值Rp和垂直偏移值Zp、以及等离子体小半径“a”的等离子体中 心坐标，
-相对于以(Rp，0，Zp)为中心的坐标系计算与电离气环的 横截面相交的视线线段的几何位置，
-计算在该视线线段上的NL个连续点P1，P2，...，PNL的位 置，P1和PNL是视线线段与电离气环横截面的边界相交的点，
-计算点Pi(i＝1，2，...，NL)和等离子体中心之间的距离 ri，以及计算归一化距离ρi＝ri/a，
-构建矩阵Jρi，其i行j列的元素是Jo(ρi*Zj)，其中，Jo(ρi*Zj) 是自变量ρi*Zj(i＝1，2，...，NL以及j＝1，2，...，NB)的0阶 函数Jo，矩阵Jρi如下：
Aρi＝Jρi*C，
其中，Aρi是表示沿归一化距离ρi的发射率分布的数组，以及C是 该矩阵的系数，
-计算Jρi的每列j(j＝1，2，...，NB)的积分Fj，如下：
Fj＝δ*∑i Jo(ρi*Zj)
其中，δ是几何常数，以及i从1到NL，
-计算矩阵F，如下：
F＝[F1 F2...FNB]
-计算矩阵F-1，其为矩阵F的伪逆矩阵，
-计算矩阵M，如下：
M＝(Jρ*F-1)/EG
其中，EG是检测器的几何延拓，
-计算超热电子的局部发射率分布数组Aρ，如下：
Aρ＝M*Y/1000
根据附加特征，本发明的过程包括检验数组Aρ的元素的检验 步骤，该检验步骤包括以下步骤中的至少一个：
-a)检验表示在电离气环的圆周上局部发射率分布的数组Aρ 的元素是否不等于0，
-b)检验数组Aρ中的任一元素是否为负数或者是否是超过阈 值的值，
-c)检验数组Aρ的局部最大数是否大于所允许的预定最大 数。
根据另一附加特征，本发明的过程包括：
-计算矩阵Jρ的伪逆矩阵(Jρ)-1，
-计算矩阵N，如下：
N＝(F*(Jρ)-1)  *EG，
-计算对应于积分等离子体发射率数据Y的重构线积分测量 值YR，如下：
YR＝N*Aρ*1000
-计算值x2，如下：
${\chi }^2=\underset{\mathrm{ri}=1}{\overset{{\mathrm{NC}}_F}{\Sigma }}{(Y_n-Y_{\mathrm{Rn}})}^2/{\mathrm{NC}}_F$
其中，Yn是积分等离子体发射率数据，其表示通过秩为n的检测器 测量的积分等离子体发射率，YRn是对应于等离子体发射率数据Yn 的重构线积分测量值，以及NCF是检测器的数量，
-检验x2是否大于固定阈值。
根据另一附加特征，在计算x2值之前，本发明的过程包括检验 数组YR的至少一个元素是否为负数的步骤。
根据本发明过程的另一附加特征，如果数组YR的至少一个元 素为负数，则过程停止，而如果数组YR中的元素均为非负数，过 程继续。
根据本发明过程的另一附加特征，如果x2大于固定阈值，则过 程停止，否则，过程继续。
根据本发明过程的另一附加特征，如果点Pi间隔相等，则仅对 位于P1和线段P1PNL的中心之间的点Pi(包括P1)、或者位于线段 P1PNL的中心和PNL之间的点Pi(包括PNL)进行点Pi(i＝1，2，...， NL)与等离子体中心之间的距离ri的计算。
根据另一附加特征，本发明的过程包括求平均值阶段，以沿着 预定数量的连续时间采样计算的原始积分发射率数据的平均数的 形式来计算该测量发射率数据Y。
根据另一附加特征，本发明的过程包括用于过滤原始测量数据 的步骤。
根据另一附加特征，本发明的过程包括初始步骤(preliminary step)，以将可用视线的数量与允许的最小值进行比较，使得如果该 可用视线的数量小于该允许的最小值，则过程停止。
本发明也涉及一种用于实时确定超热电子的局部发射率分布 的过程，该超热电子来自位于环形管内的被称作等离子体的电离气 环，该过程包括：
-读取沿相对于环形管设置的至少一个检测器的视线积分得 到的等离子体发射率的实时测量值Y，以使检测器的视线与环形管 的圆形横截面和半径为“a”的电离气环的圆形横截面，
-实时读取包括大半径值RP、垂直偏移值ZP、以及等离子体 小半径“a”的等离子体中心坐标，
-基于根据本发明的过程，实时确定局部发射率分布。
根据本发明的过程允许在Tore Supra(TS)中根据积分硬X射 线(HXR)测量数据来重构超热电子(SE)的局部发射率分布。
该过程以层析成象反演为基础，更具体地，以借助于0阶贝塞 尔函数的Abel反演为基础，该过程采用了通过当前的实时HXR TS 诊断得到的原始数据(线积分测量值)。
有利地，可以通过反馈来控制重构的局部发射率分布，以直接 影响总的电流密度分布。该算法也允许根据其形状和重构的原始数 据和测量数据之间的比较(x2比较)来检验重构的分布的有效性。
该算法也可以采用原始数据的过滤技术，以避免统计噪声。最 后，对重构的SE分布进行归一化，然后将其发送(如果满足所有 的有效性条件)至用于反馈控制的控制系统。
有利地，根据本发明的方法具有以下性能：
○鲁棒性；
○可靠性；
○为了考虑实时约束的短的计算时间。
附图说明
当读取参照附图进行描述的本发明的优选实施例时，本发明的 其他特征和优点将变得更加显而易见，其中：
已描述的图1示出了Tore Supra真空管以及环形等离子体实例 的横截面；
已描述的图2示出了根据现有技术的HXR辐射诊断测量通路；
图3示出了用于输出根据本发明的过程处理的测量数据的实时 采集电子电路；
图4示出了进入图3的实时采集电子电路的典型信号；
图5示出了图3的实时采集电子电路的具体部件；
图6示出了直方图存储器内容的实例；
图7示出了图3的实时采集电子电路的改进；
图8示出了执行根据本发明的过程的相关单元的光谱测定诊断 系统的实例；
图9A示出了位于以Tore Supra真空管的中心为中心的笛卡儿 坐标(Cartesian coordinate)系统中的Tore Supra真空管的俯视图；
图9B示出了图9A中的Tore Supra真空管的截面图；
图10示出了NB个0阶贝塞尔函数Jo的值(NB＝6)，这些值 是沿与函数Jo的前NB个零点Zi(i＝1...NB)相乘的归一化等离子 体小半径ρ的NM个等间距的点进行计算。
图11叠加地示出了NB个贝塞尔函数Jo的值，其自变量为与 NB个零点Zi相乘的从等离子体的几何中心(相对于等离子体小半 径“a”进行了归一化)到包含在等离子体中的视线L的一部分的 NL个点的距离ri；
图12示出了根据本发明的整个过程的简化方案；
图13根据本发明的过程的基本步骤；
图14-26示出了对应于图13的基本步骤的流程图；
图27示出了在TS放电(当时间t＝10s时发射数为32570)期 间使用所采用的反演法得到的SE局部发射率分布的实例。
在所有的附图中，相同的标号表示相同的元件。

根据本发明的过程处理从实时采集电子电路输出的测量数据， 该实时采集电子电路是以申请人的名义于2004年2月24日向法国 专利局提交的题为“Circuit électronique de diagnostic de spectrométrie”的法国专利申请的主题，其申请号为04 50338。
图3示出了输出在本发明的结构中处理的测量数据的实时采集 电子电路的实例。电路15a包括两个数据处理模块21、22以及可 编程逻辑接口和控制部件23。每个数据处理模块21、22通过总线 Bi连接至可编程逻辑接口和控制部件23。例如，数据处理模块21、 22可包括四个并联的输入放大器A、与四个输入放大器串联设置的 四个模/数(A/D)转换器、以及可编程逻辑脉冲处理部件PROG-I。 通过控制数据采集率的命令K1来控制可编程逻辑接口和控制部件 23。VME(VERSAModule Eurocard)总线B输出将通过根据 本发明的过程处理的测量数据。
每个可编程逻辑脉冲处理部件PROG-I对它接收的数字数据施 加一组操作，下面将结合图5的描述更详细地介绍这组操作。
图4示出了进入实时电子电路15a的典型信号。
图4中的曲线表示作为时间t的函数的信号能量E。能量曲线 E包括正脉冲成形部分和负部分。信号的“有用”部分是正部分。 正部分的持续时间是一微秒级的。由处理电子电路造成了具有几微 秒级(一般为3或者4μs)的持续时间的负部分。在图4中出现了 几个时间参数(ta，tb，tc，td，T1，T2，T3)，它们将在说明书的 其他部分中进行详细描述。
图5示出了处理模块21、22的详细示图。
处理模块21、22包括几个处理通道。为了简明并避免使图复 杂化，图5仅示出了由一个输入放大器A、一个模/数(A/D)转换 器、该转换器的增益调整电路G、以及与模/数转换器相结合的可编 程逻辑脉冲处理部件PROG-I的一部分所组成的一个处理通道。
部件PROG-I包括以下功能模块：
-脉冲检测和幅度测量模块24，
-堆叠去除(stack rejection)模块25，
-两个能量段(energy slot)分类模块26、28，
-两个数字计数模块27、29，以及
-直方图存储器30。
除了放大功能之外，输入放大器A执行阻抗匹配功能并且消除 接收信号的负部分(参见图4)。模/数(A/D)转换器量化放大器A 的信号输出。增益调整电路G用于通过VME，总线设计转换器增益。 在校准步骤期间设计转换器增益。处理模块24首先检测脉冲，然 后测量脉冲幅度。根据本发明的一个优选实施例，为了消除测量的 噪音，脉冲能量的阈值Es被用于检测(参见图4)。接收能量级大 于或等于阈值Es的脉冲，而消除能量级较低的脉冲。当脉冲已被 接收时，测量其宽度T1(参见图4)。测量脉冲宽度的起始时间是 脉冲能量升至超过阈值Es的时间ta。然后，脉冲幅度降至低于阈 值Es的时间tb被用于限定脉冲宽度T1，即写为：
T1＝tb-ta
脉冲宽度时间阈值tc被用作脉冲宽度的函数，来分类脉冲。例 如，最大脉冲宽度T2(T2＝tc-ta)可以等于1.5μs。
测量脉冲宽度的起始时间ta也是可设计时间T3的起点，在时 间T3期间，任何新的脉冲都不会被计数。例如，时间T3可以等于 5μs。例如，限定时间段T3(T3＝td-ta)的可设计时间td可以与 原始脉冲(换句话说在其负部分消除之前的脉冲)返回到0附近时 的时间(见图4)相对应。
堆叠去除模块25去除其宽度超过脉冲宽度阈值tc的任一脉冲， 并且在可设计时间间隔(例如，间隔T3)期间，如果已检测到第一 个脉冲，则去除任一新脉冲。使用未被堆叠去除模块25去除的脉 冲并通过可设计能量段来对其分类(分类模块26)对其进行分类。 例如，能量段可等于以下值：
-[20keV-40keV]，
-[40keV-60keV]，
-[60keV-80keV]，
-[80keV-100keV]，
-[100keV-120keV]，
-[120keV-140keV]，
-[140keV-160keV]，
-≥160keV。
然后，在计数模块27中对每个能量段的脉冲进行计数。例如， 在存在上述八个能量段的情况中，计数模块27可以包括八个12位 计数器，换句话说，每个能量段有一个计数器。只对与检测当前脉 冲的能量段相关的计数器加一。
同样，也通过能量段对被去除的检测脉冲进行分类，从而对所 有检测脉冲进行分类(分类模块28)并计数(计数模块29)。
然后，直方图存储器30被用于校准测量。然后，电子HXR诊 断电路进入校准模式。
现在，我们将描述校准过程。开始从已知的外部激励(标准源) 采集数据。直方图存储器30通过校准能量段来分类信号。例如， 校准能量段可以是1keV级的。其中只包括在堆叠去除后的经分类 的脉冲。进入直方图存储器的每个脉冲都对应于其能量的最大幅度 累加到一个存储单元。然后，可以搜索数量最多的脉冲设置于其中 的单元或者单元组。接着，通过VME总线，可将增益调整的操作 用于自动使该最大值与来自标准源的期望和已知的能量一致。
图6示出了直方图存储器的内容的实例。横坐标表示不同的能 量级E，而纵坐标表示对每个能量级所采集的脉冲数NI。
图7示出了对实时电子电路15a的改进。
除了上述参照图3的部件之外，实时电子电路15b包括两个下 拉(pull down)缓冲存储器M1和M2，其在它们的输入处分别接收由 处理模块21和22输出的数字数据。总线Bi将每个下拉存储器M1、 M2连接至可编程逻辑接口和控制部件23。应用于该可编程逻辑部 件23的命令K2触发将处理模块21和22输出的数据存储到对应下 拉存储器M1和M2中。例如，下拉存储器M1和M2可以存储通 过VME总线B以可配制的速率从包括在对应处理模块21和22中 的A/D转换器输出的数据的历史，或者通过总线B以可配制的速率 (该速率可能大于基本采集率)改变存储计数器27、29的状态的 历史，以便可以由此观察在两次采集之间的计数器改变。
图8示出了执行根据本发明的过程的相关单元的光谱测定诊断 系统的实例。
光谱测定诊断系统包括相机1、接收结构2、偏振电路3、电源 电路4、数据处理电路16、以及数据存储单元7。根据本发明的光 谱测定诊断系统与根据现有技术的光谱测定诊断系统的区别在于 数据处理电路16。根据本发明的数据处理电路包括与数据采集和处 理单元17以及管理单元18串联的实时采集电路15，所有这些部分 之间均串联连接。例如，实时采集电路15是图3中所示的电路15a 或者图7中所示的电路15b。处理单元17和管理单元18执行根据 本发明的过程。数据处理电路16包括共享RAM 19(RAM：随机 存取存储器)。诸如SCRAMNET卡(SCRAMNET：共享通用随机 存取存储网络)的共享RAM 19通过通信网络20方便地与其他采 集单元共享数据。特别地，共享RAM 19实时读取包括大半径值 Rp、垂直偏移Zp、以及等离子体小半径“a”的等离子体坐标。
图9A示出了位于以Tore Supra真空管的中心为其中心(0，0， 0)的正交系(orthonormal base)(x，y，z)上的Tore Supra真空管 的俯视图，而图9B示出了图9A的Tore Supra真空管的环形横截面。 Tore Supra真空管的环形横截面包括在平面(x，0，z)中。
真空管的俯视图示出了两个同心圆C1、C2(参见图9A)，而 图9B的环形横截面是由一个圆C3组成。x轴与圆C3相交于点x1 和x2，并且圆C3在正交系(x，y，z)中的中心坐标是(Cv，0， 0)。等离子体的环形横截面是圆CP。圆CP偏离圆C3的中心。在 圆CP正交系(x，y，z)中的中心O的坐标是(Rp，0，Zp)。坐 标Rp和Zp分别被称作“等离子体大半径”和“等离子体垂直偏移”。 圆CP的半径(通常被称作“等离子体小半径”)是“a”。
例如，在图9B中只示出了一条视线L。视线L与圆CP的圆 周相交于点A和B。位于视线L上的点Pi和圆CP的中心之间的距 离为ri。线段AB的中点为H。视线L与z轴相交于点I。将视线L 和x轴之间的夹角称为α，则视线的斜率T可以定义为T＝tan(α)， 其中，“tan”是正切且α具有正负号(根据其位置在x轴的下方或 上方来确定其为正还是为负)。
作为实例，现在，将在NL和NM个点是等间隔的特殊情况下 描述根据本发明的过程。为了正确理解本发明，首先必需回顾一些 基本原理和公式。从而，可以实现后续假设：
a)仅以一条视线来进行该方法的描述；
b)将视线L在点A和B之间的部分划分为NL个等间隔的点 (根据本发明的优选实施例，NL在20与30之间)；
c)将归一化的等离子体半径ρ划分为NM个等间隔的点；
d)从c)中可以得出，根据归一化等离子体半径ρ计算得到的 最终分布也划分为NM个等间隔的点；
e)贝塞尔(Bessel)函数(见下文)的数量为NB个(根据优 选实施例，NB等于6)。
图10叠加地示出了NB个0阶贝塞尔函数Jo的值(NB＝6)， 这些值是沿与Zi(i＝1，2，...，NB)相乘的归一化等离子体小半 径ρ的NM个等间隔的点进行计算，这些Zi是函数Jo的前NB个 零点(31为Jo(ρ*Z1)，32为Jo(ρ*Z2)，33为Jo(ρ*Z3)，34为Jo(ρ*Z4)， 35为Jo(ρ*Z5)，36为Jo(ρ*Z6))。在算法说明中，包含该值的矩 阵被称作Jρ。
图11叠加地示出了NB个Jo的值，其自变量为与NB个零点 Zi相乘的从等离子体的几何中心(相对于等离子体小半径“a”进行 了归一化)到包含在等离子体中的视线部分的NL个点的距离ri。 该图与前一幅类似的地方在于：实际上，如果视线与圆周(等离子 体边界)的直径一致并且如果认为其仅是直径的一半(即，半径 “a”)，则图11与图10(在相对于自身进行了半径归一化之后)相 同。由于距离相对于点H对称(见图9B)，因此，只显示一半的值。 以与图10相同的方式，41是与Jo(ρi*Z1)相关的曲线，42是与Jo(ρi*Z2) 相关的曲线，...，46是与Jo(ρi*Z6)相关的曲线，i从1到NL。在算 法说明中，包括该函数的矩阵被称作Jρi。
接下来，开始介绍公式：
Y＝N*Aρ    (1)
$N=F·J_{\rho }^{-1}·\mathrm{EG}---\left(2\right)$
其中：
-这些Y是通过相机的检测器得到的原始线积分测量值(在公 式(1)中，Y是标量，为了简明，只考虑一条视线)；
-Aρ是数组，其元素表示沿归一化等离子体半径ρ的局部发射 率分布；
-N是转换矩阵(transfer matrix)；
-EG表示检测器的几何延拓；
-Jρ是在归一化等离子体半径ρ(即，它在0和1之间变化) 上计算的0阶贝塞尔函数Jo的矩阵；
-F是矩阵，其功能将进一步说明。
除上述假设外，必须紧记以下假设：
f)作为贝塞尔函数自变量的所有长度量相对于等离子体小半径 a进行归一化；
g)Shafranov偏移(其是等离子体几何中心和等离子体磁场中 心之间的位移)可以忽略并且不用考虑；
h)局部发射率分布Aρ近似为NB个0阶贝塞尔函数Jo的叠加。
假设g)是重要的，因为它避免了计算从每个点所属的磁通量 面的中心到该点的距离(见参考文献[2])。实际上，该计算应该通 过非常耗时的迭代法来完成，而从几何中心到每个点的距离的计算 非常简单。此外，由于在系统中引入了对称性，因此，假设g)进 一步减少了计算。借助于该对称性，可以考虑将线段AH或线段 HB替代AB，如在下面将进一步说明的。
由于从中心到点Pi的距离ri被用作贝塞尔函数的自变量，以沿 着弦计算贝塞尔函数的值，因此，必须计算从中心到点Pi的距离ri。
用来保证采用假定g)的可能性的证据如下。考虑重构的局部 发射率分布Aρ。根据Aρ，可以利用公式(1)来确定重构的线积分 分布YR：
YR＝N * Aρ
因而，容易计算原始数据Y和重构的测量值YR之间的卡方 (chi-square)x2。在对不同脉冲的几种模式进行测试后，可以看到， 未考虑Shafranov位移而得到的x2平均起来相当于考虑了Shafranov 位移得到的x2：对于来自TS数据库的大约1000个采样，在第一种 情况下，x2为0.00155，而在第二种情况下，x2为0.00157(请注意 这些x2是利用归一化分布来计算的)。
为了采用假设h)，必须对NB个不同的自变量Args来计算NB 个贝塞尔函数。例如，如果我们认为Args等于ρ*Z，其中，ρ是划 分为NM个点的归一化等离子体半径，而Z是一个数组，其元素是 0阶贝塞尔函数Jo的前NB个零点，则：
Args＝ρ*Z    (3)
在图10中示出了对Args计算得到的NB个贝塞尔函数，用这 些函数构成维数为NM*NB的矩阵Jρ。Z为按照从1到NB(在本 实施例中为6)的升序的Jo的前NB个的零点。
如果设A(r/a)为距等离子体中心的距离为r的等离子体的一 个普通点P处的局部发射率分布，由于假设g)，则可认为分布A 在整个极向截面呈径向对称，以及服从假设h)，可以认为极向平面 上的每个点Pi(ri/a)上的局部发射率的每个值是以Pi(ri/a)来计 算的NB个贝塞尔函数Jo值的线性组合。因此，在极向截面上，存 在NB个为常数的线性组合系数。
对一条弦测量的线积分发射率可以表示为：
$Y=\mathrm{EG}*\underset{A}{\overset{B}{\int }}A(r/a){\mathrm{dr}}_a\cong \mathrm{EG}*\frac{D}{\mathrm{NL}-1}*\underset{i=1}{\overset{\mathrm{NL}}{\Sigma }}A(r_i/a)=\mathrm{EG}*{\delta }^{\prime }\underset{i=1}{\overset{\mathrm{NL}}{\Sigma }}A(r_i/a)---\left(4\right)$
其中，δ′＝D/(NL-1)，以及D是线段AB的长度。
由于点的数量是默认值，所以公式(4)对于合理多数量的点 特别有效，并且公式还可以改写(记住Aρi的定义)为：
Y＝Δ*Aρ     (5)
其中，Δ是维数为1*NL的行矩阵，如下：
Δ＝[δ δ.....δ]，其中，δ＝δ′*EG
由于假设h)，也得出以下公式：
$A(r_i/a)=\underset{j=1}{\overset{\mathrm{NB}}{\Sigma }}{c}_j*J_o((r_i/a)*z_j)=\underset{j=1}{\overset{\mathrm{NB}}{\Sigma }}{c}_j*J_o({\rho }_i*z_j)=J\left({\rho }_i\right)*c,---\left(6\right)$
其中，J(ρi)是根据ρi计算的NB个贝塞尔函数值；在图11中示 出了它们沿视线的典型延伸。此处，横轴上的值1意味着ri/a＝1， 从而计算在等离子边界上的J(ρi)。可以看出，J(ρi)在等离子体边界 上的所有值均趋于0(ri/a＝1乘以这些零点Zj)，并且只计算半条弦 的值。
接着，考虑到Jρi＝[J(ρ1)，J(ρ2)，...，J(ρNL)]，并且从(5)可 以得到Aρi＝[J(ρ1)*C，J(ρ2)*C，...，J(ρNL)*C]＝Jρi*C，将该公式带 入公式(5)中，可推出：
Y＝Δ*Jρi*C    (7)
考虑到公式(4)、(5)、和(7)，可以得到，数组Δ*Jρi的值就 是Jρi沿着视线的积分。上述数组Δ*Jρi在上文中称作F。因此：
Y＝F*C    (8)
根据假设h)，类似于公式(6)，可以写出：
Aρ＝A(ρ)＝Jρ*C    (9)
其中，Aρ就是期望的沿ρ的分布。根据公式(9)，可以得到：
$c=J_{\rho }^{-1}*A_{\rho }---\left(10\right)$
将公式(10)带入公式(8)后，得到以下结果：
$Y=F*C=F*J_{\rho }^{-1}*A_{\rho }=N*A_{\rho }---\left(11\right)$
记住EG已经包括在Δ中，即，在F中，可以得到，公式(11) 即是公式(1)。
因此，由下式给出未知的局部发射率分布Aρ：
$A_{\rho }=\frac{J_{\rho }*F^{-1}}{\mathrm{EG}}*Y=M*Y---\left(12\right)$
其中，从F中提取EG。
实际上，公式(12)不正确的，因为必须调整其量纲。因此， 局部发射率分布的最终公式是：
$A_{\rho }=\frac{J_{\rho }*F^{-1}}{\mathrm{EG}*1000}*Y=\frac{M}{1000}*Y---\left(13\right)$
得到公式(13)的方法由以下说明给出。
原始线积分测量值Y可以被定义为每单位能量和每单位时间 上撞击检测器的光子数量：
$Y=\frac{N_{\mathrm{PH}}}{E*t}=\frac{\mathrm{number}\_\mathrm{of}\_\mathrm{photons}\left[\right]}{\mathrm{energy}\left[\mathrm{eV}\right]*\mathrm{time}\left[s\right]}---\left(14\right)$
其中，括号中的量是物理量纲(eV是电子伏特且s是秒)。代 替地，可将局部发射率分布Aρ定义为每单位能量、每单位时间、 每单位体积、以及每单位立体角上撞击检测器的光子数量：
$A_{\rho }=\frac{N_{\mathrm{PH}}}{E*t*V*\Omega }=\frac{\mathrm{number}\_\mathrm{of}\_\mathrm{photons}\left[\right]}{\mathrm{energy}\left[\mathrm{eV}\right]*\mathrm{time}\left[s\right]*\mathrm{volume}\left[{\mathrm{mm}}^3\right]*\mathrm{solid}\_\mathrm{angle}\left[\mathrm{str}\right]}---\left(15\right)$
此处，具体地，V是以mm3为单位而Ω是以str(球面度)为 单位。考虑到公式(1)、(2)、以及(15)，EG是以mm2*str为单 位，可以将Y写为以下公式：
$Y=\frac{N_{\mathrm{PH}}}{R*t}=\mathrm{EG}\left[{\mathrm{mm}}^2\right]\left[\mathrm{str}\right]*F*J_{\rho }^{-1}\left[m\right]*A_{\rho }---\left(16\right)$
由于F*Jρ-1的量是无量纲量(贝塞尔函数)沿着以米为单位的 距离的积分(F)与无量纲量(Jρ-1)的乘积，因此，F*Jρ-1是以米为 单位。公式(16)在量纲上是不正确的，因为牢记公式(15)， mm2」*[m]≠mm3」。实际上，应该用毫米来代替米。这就是为什么 必须要使Y乘以1000：mm*1000＝m。
可以得到，Y的最终公式是：
$Y=1000*\mathrm{EG}*F*J_{\rho }^{-1}*A_{\rho }=1000*N*A_{\rho }.---\left(17\right)$
如果将公式(17)和公式(1)和(2)进行比较，可以假定在 公式(2)中因子1000已包括在EG中。另一方面，Aρ的最终公式 是：
$A_{\rho }=\frac{J_{\rho }*F^{-1}}{\mathrm{EG}*1000}*Y=\frac{M}{1000}*Y---\left(18\right)$
因此，公式(18)与公式(13)相同。
用于获得局部发射率分布的实时计算过程的总方案可如图12 所示；五个主要逻辑部分构成该方案：
●实时数据网络，用于交换数据(图8中的通信网络20)；
●原始输入数据；
●该过程必需的参数(见下文)；
●数据处理，其为实时处理的核心；
●计算的输出数据；
然后，发送至实时网络的输出数据将用于通过控制系统来实现 等离子体的反馈控制。
按照其顺序，实时算法可以划分为两个主要部分：
●参数设置，在每个实验期间之前完成；
●原始数据处理，能够根据之前定义的参数和从来自电子卡的 原始数据来获得实时发射率分布。
对于参数设置部分，操作员必须提供以下数据：
●用于选择将使用的弦的掩码(mask)(在59条弦不是都可用 的情况下)；
●截止频率，用于在数据处理之前滤除积分的原始数据；以及 阈值，用于借助于x2的计算结果来估计分布的优度(goodness)；
●每秒的光子计数(photo counts/s)的最小数，其被认为是对 于每条弦都可接受的；
●标记(flag)，选择用于滤除积分的原始数据的选项；
●标记，选择用于计算输入图案平均值(基于固定数量的时间 采样)的选项(见下文)；
●最终计算的分布的可接收局部最大值的最大数(见下文)。
对于数据处理部分，算法可以划分为三个主要部分：
●用于在每次放电的最开始初始化代码的部分；
●用于计算SE局部发射率分布的部分，按照其顺序，通过以 下主要步骤来计算该SE局部发射率分布：
○主要程序变量的初始化；
○输入数据的读取；
○数据预加工；
○坐标中心的改变；
○弦与等离子体边界交点的计算；
○矩阵F的计算；
○局部发射率分布Aρ的确定
○局部发射率分布Aρ的检验；
○重构的积分数据YR的计算；
○ Y与YR之间的卡方x2的确定。
●专门用于写入输出数据的部分。
上述内容如图13所示。
在实际执行放电之前，如果必要，通过图形界面来改变与该过 程相关的一些参数(即使一旦调整这些参数，它们将不必再改变) 使整个过程开始。在脉冲刚开始处，通过控制系统(CS)来读取这 些参数。在脉冲的开始时期，控制系统也安排初始化在该处理中使 用的主要量：这通过调用FUN_1(见图14)来实现。接着，当放 电实际开始时，CS调用FUN_2，其为在每个控制过程循环中执行 的主程序。
通过检验主要初始化阶段是否由于出错而结束，算法的主体开 始；如果主要初始化阶段由于出错而结束，那么登记事件并且 FUN_2结束(循环出错程序-CEP)，否则继续执行并且安排初始 化在循环期间使用的变量(FUN_4)。由于假设除FUN_4之外在调 用每个函数后都要进行出错检验(因此能够进行CEP)，因此从该 点直到FUN_2的主体结束，省略涉及CEP的调用。随后的步骤是 从输入缓冲器读取算法数据(FUN_5)；此后，程序必须执行一些 数据预处理(FUN_6)。然后，必须进行系统坐标的改变，其尤其 对视线与z轴的交点I产生影响(FUN_7)。此处，可以计算弦与等 离子体边界的交点IP(FUN_8)。需要用它们来计算贝塞尔函数 Jo(ρi*Zj)的自变量(见图11和公式(6))，以标准程序计算的该函 数的积分将构成矩阵F(FUN_9)。但是，为了确定分布Aρ，需要F 的伪逆矩阵F-1(见公式(13))：这可以利用标准伪逆矩阵程序来计 算。接下来，可以计算对当前循环的分布Aρ(FUN_10)，但是必须 验证刚计算的分布的连续性(FUN_11)。得到Aρ后，接着，可以 确定重构的线积分测量值YR(FUN_12)并且使YR和Y一起用于 计算它们之间的卡方x2(FUN_13)。接着，该量被用作对Aρ的优 度的最终估计：如果也通过了该检验，则可将SE局部发射率分布 用于实时反馈控制。这在构成FUN_2的写入体的FUN_3中完成， 其中，所有的数据都放置于输出缓冲器中，以将其记录在数据库中 并将其用于反馈系统。
在图27中可以看到该算法的具体结果，其中，示出了根据给 定时间(例如，t＝10s)内一次TS放电(例如，放电n°32570)的 HXR数据计算的SE局部发射率分布。
以下部分旨在说明在实时系统中执行的代码。代码由13个函 数构成；并以伪语言的形式，使用基本流程图在图13到图26中呈 现。在这些图中，函数用符号名(FUN_N，其中，N＝1，...，13) 来表示。除了图13，每个图均示出了一个函数的流程图(FUN_2 分为图15和图16)。只有FUN_4未示出：原因是该函数仅是在每 个循环的开始(即，每次调用FUN_2)简单地初始化变量。
按照以下方式将符号名与实际的执行函数相联系：
FUN_1＝initSPX(fig.14)；
FUN_2＝calSPX(fig.15-16)；
FUN_3＝writeDataSPX(fig.26)；
FUN_4＝cycleInitSPX(no figure)；
FUN_5＝readDataSPX(fig.17)；
FUN_6＝preElaborationSPX(fig.18)；
FUN_7＝newCoordCenterSPX(fig.19)；
FUN_8＝lcfsInterceptSPX(fig.20)；
FUN_9＝computeRFMatricesNoMagSPX(fig.21)；
FUN_10＝computeProfileSPX(fig.22)；
FUN_11＝checkProfileSPX(fig.23)；
FUN_12＝computeIntegratedDataSPX(fig.24)；
FUN_13＝computeChiSquareSPX(fig.25)；
图12示出了实时计算过程的概括图：需要用于交换一部分输 入数据和所有输出数据的实时数据网络。然后，标出具有原始输入 数据的块：它包含通过光谱测定诊断提供的原始HXR数据和与等 离子体几何形状相关的数据。显然，在实时数据处理阶段，通过算 法，将原始数据用作输入，但是在此之前，在每次放电刚开始时存 在可以改变某些算法参数的阶段：例如，可以更改用来计算SE局 部发射率分布的弦。最后，标出具有计算的输出数据的块。
图13示出了算法的全局视图并且示出了其所有基本步骤。在 第一步骤中，在实际执行代码的主要部分之前(即，在等离子体放 电之前)，通过图形界面将一些主要算法参数写入硬件存储器中。 第二步骤表示在等离子体放电之前的代码初始化阶段，并在图14 中对其进行了说明。第三步骤是算法的主要部分，该步骤执行实时 控制过程的每个循环(在图15和图16中对其进行详细分析)。该 步骤分为两个部分：一个部分用于计算SE发射率分布，另一个部 分用于写入将用于反馈控制的输出数据(图26)。按照顺序，用于 计算分布的部分可以分为几个分支，其中每一个均执行不同的任务： 主要程序变量的初始化，输入数据的读取(图17)，数据预加工(图 18)，坐标中心的改变(图19)，弦与等离子体边界交点的计算(图 20)，贝塞尔函数沿着弦的积分的确定(矩阵F-图21)，SE局部 发射率分布的确定(图22)，发射率分布的检验(图23)，重构的 积分数据YR的确定(图24)，Y和YR之间的卡方x2的计算以及局 部发射率分布的最终估计(图25)。
在图14中，示出了FUN_1，其被称作代码initSPX。这是在放 电前调用的函数，用来初始化一些重要的量。在其开头部分，用于 求输入平均值阶段(见图18)的一些全局变量被设置为0。然后， 假设每条视线都为直线，计算弦的正切T(斜率)和交点I。实际 上，需要这些量来计算弦与等离子体边界的交点，从而计算从等离 子体几何中心到这些点的距离(同样见图9A/9B)。在求平均值阶 段，需要存储矩阵，以防止数据成为最终N_AV图形(见用于说明 求平均值阶段的图18)。下一步骤是计算归一化等离子体半径ρ的 NM个点的值(假设c)和计算0阶贝塞尔函数Jo的前NB个零点。 实际上，需要这些零点来为Jρ(见公式(9)和图10)和Jρi(见公 式(6)和图11)准备自变量。如果该计算出现错误，则将标记 abort_from_initialization设置为1，登记事件并结束FUN_1的执行 (初始化出错程序-IEP)，否则下一步骤为计算Jρ。由于在所有脉 冲期间该矩阵都是不变的，因此，该步骤在初始化阶段完成。如果 该计算出现错误，则执行IEP，否则，计算Jρ-1，作为计算重构的线 积分测量值YR(见公式(2)和(11))所需的(常数)矩阵。如果 执行了该计算引起任何错误，则执行IEP，否则在没有出错的情况 下FUN_1终止。
图15示出了作为整个处理主要程序的FUN_2的第一部分，在 某种意义上，确定SE局部发射率分布所需的所有子程序都要通过 它来调用。如果在没有出错的情况下从FUN_1退出，则每个循环 都执行了该函数。
FUN_2可分为两个部分：第一部分涉及SE发射率分布的计算 和与其相关的量的计算(主体)；第二部分执行输出数据的写入(写 入体)。在可以或不可以执行主体、或者仅可以部分地执行主体的 同时，始终执行写入体。实际上，FUN_2开始检验标记 abort_from_initialization，并且如果该标记等于1，则FUN-2也将标 记abort_cycle设置为1，登记该事件，并且结束主体(循环出错程 序-CEP)。相反地，如果没有出错地从FUN_1退出，则存在用于 计算发射率分布的主体的第一步骤，即，初始化每个循环均更新的 参数，例如，包括通过诊断得到的原始数据的数组(FUN_4)。由 于FUN_4不能给出错误，因此在FUN_4之后，接着对原始数据进 行读取(FUN_5)，原始数据包括HXR原始数据、等离子体大半径 和小半径(RP和a)以及等离子体垂直偏移ZP。在这点之前，由于 图15中调用的每个函数可以给出CEP，所以，在该段中不再考虑 该问题。在FUN_5调用后，下一步骤是预加工HXR原始数据 (FUN_6)。后一函数的主要目的在于滤除原始数据，以减少由于 噪音引起的副作用。然后，考虑从(Cv，0，0)到等离子体几何中 心(Rp，0，Zp)的坐标系中心偏移(FUN_7)并且确定视线与等 离子体边界的交点(图9B中的点A和B)。通过FUN_8来实现该 过程，以下将详述其基本原理。FUN_2的最后步骤(参见图15)是利 用FUN_9来计算矩阵F(见公式(8))。
图16呈现了FUN_2的第二部分。同样，这里，主体的每个函 数都可以引发CEP。在计算F之后，必须计算其伪逆矩阵F-1(见 公式(12)和(13))。接着，可利用确定SE局部发射率分布的所 有手段，因此调用FUN_10。必须分析刚计算的发射率分布，以检 验其优度，这就是FUN_11的目的，而FUN_12的目的是确定重构 的线积分测量值YR。接着，该量可以与Y进行比较(FUN_13)， 以对发射率分布的优度提供进一步检验。在FUN_13之后，FUN_2 的主体结束；不论FUN_13的结果如何(是否引发了CEP)，调用 写入体的程序FUN_3，其目的在于将所有有用数据写入输出缓冲器 中。
在通过FUN_4初始化循环变量之后，在FUN_2中执行的最重 要的操作是从输入缓冲器读取算法数据。这是FUN_5的功能，在 图17中示出了其流程图。首先，代码从输入缓冲器获得原始检测 器数据，并将数据放入Y中。此处，Y包括在一个采样周期内获得 的计数值，但是重要的是每秒计数值(counts/s)，所以，第二步骤 在于通过采样时间dT来划分原始检测器数据，以更新Y。如果其 值对于计算是有用的弦的数量小于预先确定的固定最小值，则函数 调用CEP，否则继续执行FUN_5并读取当前等离子体大半径和小 半径(分别为RP和a)以及等离子体垂直偏移ZP的值。在获得了 上述三个参数中每一个之后，检验它们的值是否在有效范围内；如 果不是在有效范围内，则执行CEP。在对ZP进行检验后，函数结 束。
在从输入缓冲器读取数据后的步骤是滤除检测器的原始数据， 以避免由于统计噪声引起的问题。这是FUN_6(见图18)的目的， 其中，解决方案由以下任何一个构成：a)求输入的平均值，或者b) 利用低通滤波器来滤除原始数据。两个操作都不是强制的，而是通 过在放电之前在上述参数设置阶段期间设置两个适当的标记来决 定它们的执行。
总之，FUN_6也必须检验当前使用的弦的数量是否大于允许的 最小值(弦的最小允许数量MNC是嵌入代码中的值)，并且这是在 FUN_6中执行的第一步操作。如果弦的数量不大于允许的最小值， 则函数执行CEP，否则，FUN_6检验用于求平均值阶段的标记。
如果该标记为1，则执行求平均值阶段；它在于对确定数量NIP 的输入图案求平均值(该数目嵌入代码中并且它主要取决于物理过 程的动力学和硬件采样时间)。如果已获取的图案数量等于或者大 于NIP，则计算平均值并将其用作原始数据；否则，在没有计算任 何平均值并使用当前采样的未平均原始数据，来进行其计算程序。
在第一阶段结束时或者如果第一标记为0，函数检测第二标记， 该标记应该使阶段b)得以执行。如果用于“滤除原始数据”阶段 的标记为1，则采用可以在参数设置阶段中改变其截止频率的低通 滤波器来执行该操作。
程序进行根据可用的弦来更新所使用的数据数组。对每秒计数 (counts/s)的当前最大值进行检验(用于可信赖的反演)：如果其 小于最小允许值(该值也在参数设置阶段设置)，则FUN_6调用 CEP，否则FUN_6在没有出错的情况下结束。
按照假设g)以及考虑用于确定SE局部发射率分布所采用的 程序，可以得到，Aρ相对于等离子体几何中心呈圆周对称。与之前 已提到的等离子体几何中心不必与真空管的中心相一致的情况一 样，这也需要采用中心为等离子体几何中心的坐标系。如图19中 所示，FUN_7的目的在于改变取决于系统坐标的算法量，以使它们 与新的坐标系一致，认为该坐标系在这点之上且其中心是(RP，0， ZP)。
实际上，用于以下计算以及需要重新计算的唯一量是弦与z轴 的交点。在完成上述步骤后，FUN_7检验之前计算的结果是否落入 数字有效范围内。如果是这种情况，则在没有出错的情况下函数终 止，否则，调用CEP。
在图20中示出了FUN_8的流程图，FUN_8用于确定弦与等离 子体边界的交点A和B。实际上，根据这些点的信息，可以计算量 ri(见图9B)，并且由此计算Jρi和F(公式(6)和(7))。对于每 条弦，计算交点A和B，求解以下方程组：
x2+z2＝a2    (19)
2＝Tx+I.     (20)
其中，T和I在之前已定义。
公式(19)限定了等离子体边界，而公式(20)限定了视线。 从以上方程组推导出二次方程式(关于x的函数)，如下：
(1+T2)x2+2TIx+I2-a2＝0.    (21)
通过下式给出x的解：
$x_{\mathrm{1,2}}=\frac{-\mathrm{TI}\stackrel{-}{+}\sqrt{\mathrm{Det}}}{1+T^2},---\left(22\right)$
其中，Det是行列式：
Det＝(TI)2-(1+T2)(I2-a2)    (23)
然后，z的解可以通过公式(20)来确定。
FUN_8对每条可用弦重复以下操作：通过公式(23)来计算行 列式；如果Det大于或者等于0，则意味着视线至少与等离子体边 界是相切的，因此，公式(22)至少有一个解；如果公式(22)有 至少一个解，则利用公式(22)和(20)来计算该解，否则该函数 检验下一条弦。在该循环结束后，FUN_8验证与等离子体边界相交 的弦的数量是否大于MNC：如果这是没有发出CEP的情况，则在 没有出错的情况下函数结束。
图21描述了用于计算矩阵F的函数FUN_9。基于在调用FUN_8 之后剩余的可用弦的数量，FUN_9执行循环。在每次迭代开始时， 将包括在等离子体中的视线的当前部分(通过FUN_8计算其极值) 再划分为NL个点Pi。然后，按照假设g)，计算从等离子体几何中 心O到每个点Pi的距离ri。如果这些点等间隔，实际上，由于对称 的原因，只用计算前NL/2个点的距离(例如，这些点属于图9B中 的线段AH)。得到ri后，可以确定Jρi(见图11及其说明)。如果该 计算出现错误，则调用CEP，循环异常中断并且函数终止，否则， FUN_9确定积分F。如果该计算出现错误，则调用CEP，循环异常 中断并且函数终止，否则，FUN_9跳转到下一迭代。
图22示出了FUN_10的步骤，该函数计算SE局部发射率分布 Aρ。基本上，该函数是简单的：首先，它确定矩阵M(见公式(13))； 然后，如果前一计算返回错误，则调用CEP并终止函数，否则，通 过公式(13)来计算分布Aρ。
在计算分布Aρ后，必须检验分布，以了解它是否是合适的分 布以及将其用于实时目的用途是否可靠。虽然在信号调节的初级阶 段，但是原始的Y测量值仍可能存在由于硬件或者软件故障而产生 的错误值，从而导致Aρ的错误计算。FUN_11的任务是检验分布Aρ， 在图23示出了其流程图。它由三段构成，并且如果检验到的特征 不令人满意，则其调用CEP。该检验过程如下(值得一提的是，最 后的检验结果是从实验分析中推导而来的)：
-1)边界上的Aρ的值不为0？(实际上，在公式(9)中所使 用的所有贝塞尔函数在等离子体边界上为0-同样见图10-因此， 这些函数的任何组合在该点上都必须为0)；
-2)Aρ的任一值小于0(Aρ＜0)或者太大？(SE局部发射率分 布是正物理量，因此，其不能小于0)；
-3)Aρ的局部最大值大于最大允许值NLM？(通常，NLM等 于4并且在参数设置阶段是可更改的)；
实际上，Aρ的在先估计只是交叉检验过程的第一步。第二步包 括估算Y和重构的线积分测量值YR之间的卡方x2(见图25)。因 此，在已检验了分布并且证实其是可信的之后，必须根据Aρ来计 算重构的线积分测量值YR。这是FUN_12的目的，如图24所示。 首先，确定矩阵N(见公式(2)和(11))；接着，通过公式(17) 来计算测量值YR。如果YR表示任何小于0的值，则函数发出CEP 并终止，否则它在没有出错的情况下终止。调用CEP的原因仅仅是 因为YR被再生为一个正的物理量，因此，它不可能小于0。
在图25示出了FUN_13的流程图。该函数计算之前提到的Y 和重构的线积分测量值YR之间的x2(NCF是用于确定Aρ的弦的最 终数量)。在此之后，FUN_13比较x2和阈值，在参数设置阶段可 更改阈值。如果x2大于阈值，则结果是为了实时的目的放弃当前的 Aρ，并且发出CEP。在将Y和YR都归一化为1后，不必计算x2。
随着调用FUN_13，FUN_2的主体结束(见图16)，此处，基 本上可以认为与用于确定Aρ的算法严格相关的所有计算都结束了。
在执行了FUN_2的主体后，得到且大量相关的数据必须写入 输出缓冲器中，以将数据存储在数据库和/或为实时的目的使用它 们。这通过FUN_3来实现(见图26)，该函数仅是FUN_2的写入 体的函数(同样见图16)。
由于归一化的SE局部发射率分布满足实时控制算法，因此， 第一步骤用于将Aρ归一化在0和1之间。接着，为了将分布Aρ记 录在数据库中以及为了将其用在实时控制方面，将分布Aρ写入输 出缓冲器中。值得注意的是，如果在调用FUN_10之前发出了CEP， 则写入到输出缓冲器中的Aρ值是从FUN_4中得到的。同样的程序 应用于重构的线积分测量值YR。最后，将其他相关参数(例如，标 记abort_cycle，卡方x2，以及CEP事件数)发送至输出缓冲器。在 这些参数之中，标记abort_cycle对于检验分布有效性尤其重要：该 标记值为0意味着当前的SE局部发射率分布可用于反馈控制，而 值1意味着当前的分布是不可用的。在后一种情况中，将在取决于 物理结果的固定时间间隔(≈100ms)内获得的最后的有效SE局 部发射率分布发送至控制系统。如果在该时间间隔内没有有效分布 可用，则执行用于脉冲终止的适当程序，其不在本专利的保护范围 内。
图27示出了根据在给定时间(t＝10s)内一次TS放电(放电 n°32570)的HXR数据、通过采用的反演法来计算的SE局部发射 率分布的实例。
引用的参考文献
[1]″Tomography of the fast electron bremsstrahlung emission during lower hybrid current drive on TORE SUPRA″，Y.Peysson and Frédéric Imbeaux， ，Rev.Sci.Instrum.70(10)，1999，pp.3987-4007.
[2 ]″Tokamaks″J.Wesson，Clarendon Press (Oxford)，1997，pp.108-121.