1 引言

自1970年世界上第一台电离层加热装置Platteville在美国Boulder（40.18°N，104.73°W）建成至今，世界各地又相继建成了众多的电离层高频人工变态装置，包括美国的HIPAS（以及后来的HAARP），俄罗斯的SURA，欧洲EISCAT计划中的Tromsø等等.利用这些大功率加热机得到了很多有价值的实验结果，包括加热产生的一级时空效应^{[1, 2]}和二级时空效应^[3~5].

在高频电磁波加热电离层中，泵波反射高度处折射率接近于0，使其电场增幅明显增大^[6]，该因素引起的偏移吸收导致泵波反射高度附近区域成为除低电离层（D、E层）外的另一个主加热区^[7].偏移吸收会引起该区域电离层内泵波反射高度处出现中大尺度的电子浓度和电子温度的变化.

关于加热中电子浓度、电子温度变化的理论以及数值模拟研究已有作者讨论过^[8~10]，其结果多为一维数值模拟，对于扩散系数和热传导系数都很小的D，E层是基本适用的，但在场向扩散系数比横向传导系数高数个量级的F层则存在明显的偏差，此时电离层被加热区域与泵波直接照射区域明显不同，加热导致的电离层空间形态也将有所不同，因此建立完整的高频加热电离层二维模型将更有助于对加热效应的深入理解和阐述.

对高频泵波加热电离层理论模型及数值模型的详细介绍和一般情况下的二维模拟结果已经在文献[11]中给出.模拟加热结果表明，当高频泵波入射到电离层中时，沿着穿过泵波反射区域的磁力线方向将会出现三个区域明显的电子浓度变化，如图 1.其中泵波反射高度附近电子浓度迅速下降，图中标为1区域，而在逆着穿过此区域的磁力线方向（2区域）与沿着穿过此区域的磁力线方向上（区域3），则会分别出现电子浓度稠团.（为了方便下文的表述，后文模拟结果中同样可以将加热区域分为与图 1类似的三个区域）.其中逆着磁力线方向上的稠团不甚明显，且会很快消失，较少被实验观测到，泵破反射区域的空洞与沿着穿过此区域的磁力线方向上的稠团则得到文献的证实^{[8, 12]}.加热中电子温度上升幅度随初始加热条件不同而介于2 5%~9 0%之间，最大电子浓度稠团和空洞的相对变化率分别为10%和-13%，其结果与文献[13]所获结果基本一致.此外，模拟结果中电离层电子温度、电子浓度的饱和与冷却时间尺度也与EISCAT雷达观测到的Tromsø加热实验结果基本相同^[12].

尽管前文介绍的模拟结果能够基本反映众多泵波加热电离层实验结果，然而在对高频加热中泵波反射区域进行诊断观测中，时常能够观测到一些不同寻常的电子浓度、电子温度变化^{[14, 15]}，如Arecibo就曾在高频加热中观测到大尺度区域内+300%电子温度、-25%电子浓度变化^[16]，曾有作者利用自聚焦不稳定性理论来解释这种反常的加热结果^[17].另外在一些高频泵波加热电离层实验中^{[18, 19]}，电子温度、电子浓度随时间的变化更加复杂，电子温度在加热机开关机一段时间后并没有达到稳定，以及电子浓度在加热中的不规则变化，本文介绍的不稳定机制或许可以解释以上实验结果.

前文在利用加热模型进行模拟加热中^[11]，虽然不同电离层电子浓度、电子温度的空间结构对模拟加热结果有一定影响，但是在大部分平静电离层背景下模拟结果基本上是一致的.然而当模拟加热的初始电子浓度、电子温度结构满足一定条件时，模拟加热电离层泵波反射区域出现一种正反馈机制，这一机制使得电子浓度、电子温度随加热机开关机的响应出现明显不同，且电子浓度、电子温度的变化也更加强烈.

本文首先利用文献[11]中的加热机制理论和数值模拟方法，以中纬度平静电离层为背景，进行地基高频电波加热的数值模拟.然后从电离层加热的物理机制的角度，结合具体模拟结果，逐项分析加热的物理模型，对地基高频电磁波加热的数值模拟结果正反馈机制进行分析，找出导致该特殊模拟加热结果的源项，并对模拟结果进行解释.

2 模拟结果

本文使用的模拟加热模型如下^[11]:动力学方程:

(1)

连续性方程：

(2)

电子能量守恒方程：

(3)

上式中k_B表示玻尔兹曼常数，D表示等离子体扩散系数，T_e、T_i分别为电子温度和离子温度，m_α、n_α、u_α（α可以为电子，NO⁺，O₂⁺，O⁺）分别为电子或离子的粒子质量、浓度和扩散速度，P_α，β_α分别为离子的产生率和消失率，k_e表示热传导系数，Q_HF表示泵波能量源项，L_e表示电子能量损失项，S_α，Q₀是为使各参量稳定于初始状态而引入的额外源项.

尤其值得注意的上式Q_HF，其大小直接影响模拟加热结果.此时Q_HF不但包括欧姆加热，还包括波-波、波-粒相互作用引起的反常吸收.其具体计算方法如下：首先做以下假设，电离层对泵波能量吸收在高度和水平方向上都满足以反射点为中心的高斯分布，前者可以由泵波天线方向图大致证明，后者可由实验观测知其大致满足高斯分布^[19].泵波反射点处源项Q_HFmax为

(4)

其中Q_L为欧姆吸收，反常吸收Q_anom利用饱和光谱方法^[20]计算而来，电场强度E₀利用经验公式E₀=求得，E_t为反常吸收激发阈值.ε₀为介电常数，ω，ω_p，ν_e分别为电波频率、等离子体频率，电子碰撞频率.

模拟加热地点仍为武汉（114.7°E，30.5°N），地磁纬度：18.8°，磁倾角：45°，加热时间2007年9月1日，加热有效辐射功率（ERP）100 MW，以上模拟加热条件与文献[11]完全相同，仅将加热时间选为LT16:00，高频泵波的频率为8.25 MHz，其时电离层临界频率为8.52 MHz，泵波反射高度同样在250 km附近.加热时间为300 s，并观察关机后600 s内T_e和n_e的变化.为了方便表述，规定位于加热机正上方泵波反射高度处的点称为加热中心（位于图 1中1区域）.模拟加热背景电离层模型选自IRI-2007，其电子浓度、电子温度二维剖面见图 2，大气模型选取NRLMSISE-00.

加热机开机不同时间后电子浓度电子温度的数值模拟结果如图 3所示.由图 3a，图 3b可知，+2 s后（+2 s表示加热开始时刻以后2s，后文与此相同，超过+300 s表示已经关机）电子温度上升幅度较小，加热中心处电子温度上升至2019 K，变化率为11.5%，比LT20:00^[11]加热2 s后电子温度42%的上升小的多，此时电子浓度没有明显的变化.造成以上差异是由于本文模拟中泵波反射高度处电子温度为2000 K，比夜间对应高度处1000 K高得多，其电子能量更加容易损失.+20 s时（图 3c，图 3d），电子温度变化率变为30 %，也比LT20:00加热时的78.5%要小.此时电子密度的场向分布也没有夜间明显.随着加热持续，电子温度继续上升，在+120 s后（图 3e，图 3f），电子温度上升速度减小.电子浓度、电子温度均出现明显的场向结构，场向不均匀体纵向尺度更大，几乎整个模拟平面分为左右两部分，尤其是电子温度1578K等高线（1662 K下面一条，相邻等高线间温度之差均为84 K），高度从270 km一直延续到385 km，加热使得原本由电子温度谷隔开的两个电子温度峰连接在一起.仔细对比此时加热中心以及电子温度谷值处电子温度的变化可知，它们分别从1850 K、1306 K升高到2362 K、1909 K，谷值处电子温度的相对变化量反而更高.对比图 3e、3f，图 3g、3h可知，此时电子温度、电子浓度仍然没有稳定（更加详细的电子温度、电子浓度随时间变化的分析将在第3节给出）.+600 s时（亦即加热机关机300 s）的电子浓度、电子温度空间剖面如图 3i、3j.此时电子温度有所降低，其最大相对变化率由+ 300 s（加热机关机时刻）时的52%变为28%，但没有如文献[11]模拟加热那样在加热结束后电子温度迅速下降到初始温度、电子浓度在加热结束后逐渐恢复到初始状态，其相对变化率反而由加热结束时的-7.5%降低为-10%.+900 s后，电离层电子温度比+600 s时进一步升高，而电子浓度则进一步降至-16%，另外电子温度最大变化率所在高度有所上升.

3 理论分析

在相同的物理模型和模拟方法下，使用基本相同的模拟条件，得到的模拟结果与文献[11]的结论差异明显.由于高频泵波加热电离层是一个非线性过程，涉及变化的参数众多，要找到导致上述结果的原因，完全通过公式解析推导非常困难.本文使用的方法是首先在模拟平面内选取部分有代表性的点，从模拟结果中导出部分参数值，结合加热模型，定量分析这些点上电子浓度、电子温度以及其他参数随时间的变化规律，并从中找出导致以上模拟结果的物理机制.

由于本文主要模拟区域为电离层F层，垂直于磁力线扩散系数相对很小，模型仅考虑了沿着磁力线方向上的扩散和热传导.为了更清晰地了解电子浓度、电子温度随时间变化，图 4，图 5分别给出经过加热中心所在磁力线上11个等间距点处（这11个点相对加热中心对称）电子温度和电子浓度相对加热前变化量随时间的变化（三维视图及侧视图）.

从图 4中可明显看出：加热开始后，所选各点的电子温度相继出现不同幅度的上升，其中加热中心顺磁力线往下区域内（类似图 1中3区域，为了表述简便，下文简称3区域）各点电子温度响应时间最短，此处电子温度的饱和时间约为40s，而加热中心及其逆磁力线往上区域（下文简称2区域）各点处电子温度并没有饱和，还存在持续的小幅上升，且离加热中心越远的点随着泵波能量传导时间的增加，对加热机开机响应也越慢.而文献[11]模拟加热中，各点电子温度虽然饱和时间长短有别，但最长没有超过70 s.本文模拟加热结果的另一个重要特性就是加热结束后，图 1中三个区域，电子温度达到饱和，各点处电子温度迅速恢复到加热前状态，而2区域内几点处电子温度经过最长120 s的下降后（离加热中心越远，下降幅度以及下降速度越小）开始上升，其中最上面几点处的电子温度经过足够长时间的缓慢上升后甚至超过了加热结束前的电子温度，而文献[11]加热中所有点电子温度冷却时间均介于50~100 s间.

随着高频泵波对电离层的加热引起电子温度上升，导致等离子体膨胀，此时从图 4看出加热中心（及其上两点）处的电子浓度下降，其沿磁力线上下区域电子浓度均开始上升，出现上下两个电子浓度稠团，上稠团（图 1中区域2，下文类似）由于受到重力的作用在50~100 s内消失.另外加热中心顺磁力线往下点处电子浓度虽然有稳定的趋势，但是直到加热结束时，电子浓度仍没有稳定.

加热结束后，关于电子浓度的变化仍然可以分为加热中心上下两个区域，2区域电子浓度在+900 s时可以认为稳定，而其以上区域内电子浓度则继续下降，且下降速度不断缓慢变大，这一下降持续到加热结束后很长时间.值得注意的是此时加热中心附近两点（+18 km、0 km）处电子浓度在经过300 s的向初始状态回归（电子浓度上升），随后却便变为与3区域内几点处电子浓度同样的下降趋势.这与文献[11]中电子浓度在关机后都会向加热前大小回归，并在经过约600 s后稳定明显不同.

综合以上模拟结果，特别是图 4中+18 km、0 km两点处电子浓度在加热结束后先上升后下降的异常现象可知，电子温度在加热机关机后没有恢复到初始状态，并在之后开始缓慢上升，这是引起异常变化的主要原因.要解释以上异常现象，必须从电子能量方程入手.公式（5）、（6）分别给出加热开始前和加热结束后电子温度的平衡方程，这两个时刻泵波源项Q_HF皆为零：

(5)

(6)

其中的下标I，H分别表示加热开始前和加热结束后各项.

将（6）式减去（5）式，消去Q₀后，得到加热结束后电子温度随时间变化率的方程，然后根据电离层电子温度、电子浓度等参数的空间剖面分布特点来分析模拟结果.为了使以后的计算尽量简便，考虑到电子能量公式（5）、（6）中的对流项、压缩项对电子温度的影响很小，忽略以上两项，此时公式变为：

(7)

其中.

从上式可知，加热结束后电子温度变化由右边R₁，R₂决定，他们分别对应于热传导项和能量损失项相对加热前的变化.下面来详细地分析右边两项的大小.

由于模拟附近区域内的中性气体对热传导的影响很小，忽略其影响后热传导系数公式变为

(8)

其中a为常数，可以看出k_e，T_e是正相关的.

对于第一项R₁，电子温度的变化能同时影响的大小，设：

将以上几式带入（7），得到

(9)

（9）式右边各项分别表示电子温度二阶梯度变化、扩散系数变化、电子温度二阶梯度及扩散系数变化协同项（下文将分别用D，E，F表示），由于上式每个时刻的值依赖于初始以及当前时刻电子温度、电子浓度的分布，另外式中参数均为对磁场方向求梯度.下面直接利用从原有的对应时间的模拟加热结果中抽取出一条磁力线上的各种参量求不同时刻的R₁以及其3分量D，E，F的大小，以解释模拟加热的异常现象.下文对公式（7）右边第二项R₂的分析使用相同方法.

选取得到的等间隔点同样位于经过加热中心同一条磁力线上，且以加热中心为对称中心，点间距为2 km，共131个点.图 6、图 7分别给出这131点处+320 s、+600 s时电子浓度、电子温度，其从另一个侧面反应了模拟加热结果.从图中可知，电子温度上升最大并非位于加热中心处，而是初始电子温度谷值处，此处电子浓度相对下降率也最大.

图 8和图 9分别是基于+320 s、+600 s时电子浓度电子温度计算得到的选取点上R₁和D、E、F的大小.对比图 8、9可知，随着加热机关机时间增加，R₁及其分量有所减少，这是由于此时电子温度开始向初始电子温度回归所致，可见直到此时，R₁及其各分量仍然是在向初始状态趋近的；其次D，F在正负交界点从关机20 s时的-10 km移动到了+10 km附近，这进一步导致了R₁最小值位置向右小幅移动.

综合上述分析有以下结论：在+320 s、+600 s，图 4模拟加热中电子温度出现异常的-30 km到+90 km以上区域内，R₁小于零，也就是该区域内电子温度上升引起了通过热传导损失的电子能量的增加.在模拟加热中，当电子温度被加热后R₁将抑制电子温度的继续上升，在当加热结束后，其会进一步引起温度的下降，可见R₁并不是导致电子温度异常的原因.

下面再来分析R₂的变化，R₂是加热前后由于电子与中性粒子，阳离子间复杂的相互作用引起的电子能量冷却率之差，其中L_e由众多部分组成，数值模型中考虑了以下几种较重要的能量冷却机制：（1）电子-离子（包括电子与NO⁺，O₂⁺，O⁺）弹性碰撞，（2）电子中性粒子（N₂，O₂，O，H_e）弹性碰撞，（3）中性分子的N₂，O₂旋转能级的激发，（4）N₂，O₂振动能级激发，（5）原子（O）电子能级激发，（6）原子（O）精细结构的激发等.但他们在高频泵波加热电离层中对电子温度的影响大小各不相同，经计算，其中1，4，6三项影响比其他项大得多，尤其在电离层较高高度处可以忽略其他因素的影响，这三项表达式分别如下^[18]：

L_e1电子与NO⁺，O₂⁺，O⁺碰撞：

(10)

L_e4电子与中性分子N₂，O₂的碰撞：

(11a)

(11b)

其中：

L_e6：O原子精细结构激发：

(12)

设L_E1=L_e1H-L_e1I，L_E4=L_e4H-L_e4I，L_E6=L_e6H-L_e6I，R₂=L_E1+L_E4+L_E6，下标I，H分别表示加热前后的值，图 10分别给出热机开机后不同时刻的131个点处L_E1、L_E4、L_E6以及他们之和R₂，以及对应时间计算得到的R₁的值.

从图 10a可知，+20 s时，随着电子温度的强烈上升，电子与各中性粒子、阳离子等相互作用频率增加，电子能量冷却率增加，所以R₂及其分量L_E1，L_E4，L_E6均大于零.另外由公式（11a）、（11b）、（12）可知L_E4，L_E6与N₂，O₂，O浓度以及电子浓度成正比，中性粒子浓度随高度呈指数衰减，而电离层最大电子浓度在+10 km附近，L_E4，L_E6最大值必将在+10 km以下某处.这与图 10a结果一致.而L_E1受到电子浓度与T_e-T_n共同的影响，其最大值位于+40km处.此时3区域内电子温度基本饱和，随着加热的继续，2区域内电子温度将继续上升，T_e-T_n的进一步增加使得L_E6变大，直到在一段时间后L_E6出现双驼峰.

从图 10b中可以看出，+300 s时，R₁在+30km以上区域内基本上与R₂相等，如果加热机此时不关机，加上泵波源项Q_HF，此时必然有

该区域内电子温度将继续上升，这就是为何图 4中电子温度直到关机前仍然在上升的原因.此时原来位于+40 km附近区域的L_E1峰值消失了，取而代之的是加热中心附近区域的凹陷，亦即该区域内与电子和阳离子弹性碰撞有关的能量冷却反而随着电子温度上升后减小.对比此时加热中心附近区域R₁，R₂大小可知，R₁-R₂明显小于零，当加热机在此时关机后，该区域内电子温度迅速下降.

+320 s时，随着电子温度下降，L_E4，L_E6逐渐减小，L_E1的作用开始凸现，R₂在-10 km以上区域内已经明显受到L_E1影响而减小，不过整体上R₂在所有点上仍大于零.几乎选取的所有131个点处R₁-R₂均小于零，这与图 4电子温度在关机20 s后下降相一致.+600 s时L_E1，L_E4，L_E6进一步减小，其中包括已经为负值的L_E1，导致此时所有选取点上R₁-R₂＞0，本还没有回复初始状态的电子温度开始上升.

下面结合加热中得到的各分量的值对下午加热中的异常现象进行解释.虽然本次模拟加热中泵波反射高度不高，但是由于此时电离层明显的电子温度双峰结构使得电子温度变化最大处位于320 km，随着高度的上升，N₂，O₂，O浓度呈指数衰减，与这三种浓度正相关的L_E4、L_E6开始变小，1区域内电子温度主要受到L_E1的影响，其开始主导电子温度的变化.

由上文可知，3区域内电子温度的所有异常变化都是由于L_E1没有随电子温度明显加热而上升甚至在某些区域内有所下降引起的.对比相同模拟条件下早上、下午、晚上不同时间模拟结果，引起下午L_E1特殊变化的因素如下：首先从初始电子温度图 2中可以看出，所选取的下午电离层内电子温度的高度剖面呈双驼峰结构，上下驼峰处电子温度均为1970℃，对应高度分别为450 km和230 km，谷值处电子温度相对峰值处比早晚更低，约为峰值温度的2/3，高度约为300 km，更低的谷值电子温度导致其在加热后更强的上升，+300 s时其电子温度最大变化率为46%，而LT10:00^[11]相应变化仅为22%.由图 6可知，虽然此时泵波反射点处电子温度仍为最高，但其绝对变化量反而比谷值处小.其次下午电离层反射区域内电子温度在加热前后都比文献[11]的更高，由图 11可知，温度越高，公式（10）中（T_e-T_i）·T_e^-1.5项随温度的变化越小，在电子温度超过2500℃时，该项随电子温度的上升将保持不变甚至下降，L_E1必然小于零.文献[19]就曾经提到过在高频加热中在电离层上F层内电子-离子的碰撞频率与电子温度之间存在复杂的依赖关系，当电子温度处于较高温度时，电子碰撞频率会随着电子温度的上升而下降.

由以上分析可知，下午特殊的电子温度空间结构才是模拟结果异常的根本原因，虽然并不是所有下午高频加热实验中都会出现类似不稳定性，但不可否认，此时电子浓度、电子温度特殊的空间结构更利于正反馈机制的形成.

4 结论与讨论

综合以上各模拟与分析结果可知，当电离层电子浓度、电子温度满足以下两个条件：第一，电离层电子温度随高度呈明显的双峰结构且电子温度谷值所在高度较高；第二，泵波的反射高度正好位于电子温度下峰值高度附近.一种正反馈机制形成，不稳定性被激发.此时电离层在高频泵波的加热下将会产生非常强烈的变化.通过对模拟结果的详细分析还可以得到以下结论：

（1）在加热泵波反射区域逆磁力线往上区域内，电子-离子的碰撞成为电子能量主要损失机制，当电子温度较高时，电子-离子碰撞引起的能量损失量基本不随电子温度的升高而升高，而是敏感的依赖于电子浓度的变化.

（2）由于已经处于加热机泵波照射区域以外，加热中心逆磁力线往上区域内电子温度、电子浓度受加热泵波的直接影响较小，而仅仅与电子热传导项和电子离子碰撞能量损失项有关，其中电子热传导项与电子温度的二阶梯度与热传导系数有关，与电子浓度变化无明显关联.

由于下午电离层背景电子浓度、电子温度都较高，且D、E层对泵波的吸收更强，电离层条件对高频加热电离层实验有一定负面影响，下午进行的加热实验较少，无法将其结果与本文模拟结果进行对比.但是文献[17]曾对在Arecibo进行的加热实验中加热中心逆磁力线往上区域内电子浓度、电子温度做过较详细的诊断观测，实验分析的结果与以上结论1、2基本相同.当电离层满足一定的条件时，在1、2结论的基础上正反馈机制引起的不稳定性是可以预期的，以上结论对今后的加热实验具有一定的指导作用，即在一些特殊的电离层状态下的电离层高频加热可能会得到非常强烈的加热效应.

为了简化模拟过程，本文没有考虑离子温度在加热过程中的变化，首先是因为其变化较小^[16]，另外从结论可知，离子温度的升高只会进一步的减弱电子-离子碰撞引起的电子能量损失.另外本文的模拟利用了文献[7]中对反射区域内电场估算方法，该方法没有考虑下午电离层D层，E层对泵波能量更加强的吸收，该吸收可能减小反射区域内的加热效应，这一因素将在以后的进一步研究中加以考虑.