对短波通信选频的思考 0 前言 最近几年接触并参与了几个短波通信选频的项目，从算法实现到服务界面设计，都有过深度参与的经历。但受限于各种因素，每次的设计在我看来都不是那么的完善和人性化。基于自己的经验和一些思考，简单谈下自己的一些认知和想法。 1 什么是短波通信选频 1.1 简要背景 1901年古列尔莫·马可尼（Guglielmo Marconi）首次实现跨洋通信的壮举，直接推动了无线通信技术的发展与实用化。随着现代对电离层反射的研究，短波（3–30 MHz）因F层反射效率高、损耗低，成为跨洋通信的主流频段。 1.2 简要原理 短波通信实际上是利用地面-电离层形成的类似波导管的结构，使得电磁波在电离层和地面之间来回反射，从而传播到几百甚至上千公里的距离。 1.3 选频场景 可以用以下两种场景覆盖选频的需求 当我们知道发射站、接收站位置以及电离层参数信息时，我们就可以通过理论或经验公式计算得到两地之间的可通频率范围，从而再根据理论或经验计算得到最佳的通信频率。 当我们知道发射站位置、通信频率以及电离层参数信息时，我们就可以通过理论或经验公式计算该频率可以传播到的空间范围，再根据损耗等参数信息，得到当前通信频率的最佳通信位置区间。 1.4 需要哪些技术和数据 电磁波传播路径计算：射线追踪技术、QPS模型下的传播路径公式、马丁等效原理、国际电联组织（ITU-R）提供的经验公式等 用途：根据站点位置、通信频率、电离层参数等信息，理论计算或经验估算电磁波的传播路径，从而判断当前频率是否有效 比较：从我的经验来看，以上这些技术所提供的结果不会有太大的区别，通常来说选择一种即可 成熟度：以上这些技术都比较成熟 电磁波损耗计算：这一块目前我只使用过ITU-R提供的经验公式进行过计算 用途：辅助判断电磁波的损耗是否严重，具体哪个可通频率的损耗最小 成熟度：在ITU-R官方文档中可查询，比较成熟 电离层参数：根据使用的技术不同，需求的具体参数也不同 射线追踪：通信时刻的三维电离层网格化数据（经度$\times$纬度$\times$高度） QPS模型下的传播路径公式：通信时刻的三维电离层网格化数据（经度$\times$纬度$\times$高度） 马丁等效原理：通信时刻反射点位置的电子密度剖面（若需要计算任意链路，约等于需要通信时刻的三维电离层网格化数据） ITU-R：通信时刻反射点位置的foF2、foF1、foE和M(3000)F2（若需要计算任意链路，约等于需要通信时刻前述各个参数的二维Map数据） 重要程度：由于选频技术已经非常成熟且区别不会很大，所以输入数据的好坏是影响选频结果好坏的关键所在。好的数据才能产出好的结果，差的数据什么技术也不好使。 分类：在应用场景，经常会出现现报和预报，选频的现报和预报则只与数据的现报和预报有关。现报：通常利用探测设备实时探测的电离层参数信息，实时同化处理成所需要的三维或二维电离层网格化数据，由于数据处理的复杂性，必然会存在一定的时间延迟；预报：通常利用已有的电离层模型，通过改变输入时间等参数，让模型输出所需要的三维或二维网格化数据，众所周知，目前的电离层模型准确度都不够高，所以预报结果只会是一个大致的估算结果，在电离层平静期可能不会有较大差距。 成熟度：无论哪种技术，对电离层参数的需求基本都是通信时刻全球（或重点区域）的三维或二维网格化数据。即便是现报，由于电离层探测设备不可能全空间覆盖，生成指定时刻准确的三维或二维网格化数据也是一件非常困难的事情。 2 工程化技术的选择 虽然前面提到了电离层三维网格化数据获取的困难性，但在一定误差范围内还是能做一些事情的。 2.1 选频技术的选择 首先，个人对纯经验模式是没那么推崇的，虽然可能效果短期内看起来还不错，但上限很有限且不够物理，所以排除ITU-R经验公式和马丁等效定理；其次，工程化过程中计算效率非常重要，过于复杂而导致时效性很差是得不偿失的，所以排除射线追踪技术（叠甲：对于优化很好、计算效率很高的射线追踪代码，肯定是优先考虑的）。那么就只剩下了QPS模型下的传播路径计算公式。 QPS模型是将电子密度剖面分为一段一段的抛物线，通过最小二乘拟合，得到一组最符合实际剖面的分段解析公式。将电子密度剖面解析的目的是对电磁波传播路径进行积分计算时，方便积分得到解析解。 QPS模型是上世纪提出的方法，从我目前的认知来看，或许完全不需要弄得这么复杂，当前数值方法很多且很成熟，直接使用数值解可能会更加的简单有效。 2.2 电离层参数的获取 现报 对于存在观测数据的情况，对观测数据进行同化处理，得到三维或二维网格化数据是一种很容易想到的方向。我对同化技术没有过研究，只了解之前项目用过的一种非常简单的方法。不过，同化技术已经发展了很多年，应该是有很多比较成熟的论文可以参考的。 预报 从我的经验来看，预报往往才是用户最想要的内容。很多时候，我们都是想提前知道未来某个时刻的通信状况如何。但这一点却是这个工程中最难的一点。 一种比较简单的处理方法是直接使用某个电离层模型未来时刻的输出作为对电离层参数的预报，这个方法在之前的项目中采用了很多次。 另外一种可能的方法是利用已有的电离层参数信息，采用深度学习等手段，训练生成未来时刻的电离层参数，就是不清楚对于三维空间网格，这种方法的难度和计算量会不会太大。 还有一种可能的方法是通过预报电离层模型输入参数，然后利用预报的模型输入参数驱动模型，从而生成预报的电离层参数，但我也从未对此进行过研究和调研。 结合现报的一种可能方法是通过预报观测台站位置的电离层数据，得到预报的台站数据，然后利用现报的同化模型，使用预报的台站数据同化得到预报的三维或二维网格化电离层参数。 2.3 选频结果 短波通信计算中有很多频率，如最大可用频率（MUF）、最低可用频率（LUF）、最高可能频率（HPF）、最佳工作频率（OWF）等，从我的经验来看，计算MUF和LUF就足够了，太多的频率参数反而会让用户很困惑。 除了频率外，大圆距离、传播模式、时延、损耗等也是比较重要的信息，能提供一定的辅助判断，最好也是能够输出的。 3 界面设计及展示思路 3.1 选频界面 对于1.3节中提到的第一个选频场景，可以设计如下界面： 通信时间：能够选择到某一天即可，不需要具体到时刻 发射站位置：可以直接输入经纬度坐标，也可以在地图中选点 接收站位置：可以直接输入经纬度坐标，也可以在地图中选点 开始计算：点击后，将上述表单信息传递给后端算法处理，自动计算通信时间那一天24h所有时刻的选频结果，时间分辨率最低1小时最高15分钟即可 折线图展示框：横坐标0至24h，纵坐标频率，折线图展示计算结果即可。如果能够计算损耗，采用二维彩图的形式将LUF到MUF之间的损耗或信噪比展示出来，效果会更好。 一次性计算24h的选频结果参考了VOACAP的设计思路，同时选频有明显的日变化，一天的数据也更有利于展示，使展示结果更加的美观。 对于现报和预报的区分，我觉得不需要在界面上显示出来，设定好时间后，让后台算法自动区分识别即可，这种处理可以避免前端过多冗余、复杂的选项。 3.2 指定频率传播范围 对于1.3节中提到的第二个选频场景，可以设计如下界面： 通信时间：需要具体到某个时刻，与电离层数据的时间分辨率有关（如30分钟） 通信频率：指定某个通信频率，3～30MHz之间 发射站位置：可以直接输入经纬度坐标，也可以在地图中选点 开始计算：点击后，将上述表单信息传递给后端算法处理，自动计算通信时间发射站向四周辐射的选频结果 二维地图展示：以发射站为中心，展示四周可通信范围，大致应该是个环状结构。如果可以计算损耗或信噪比，则对可通范围进行上色处理，效果会更好。 3.3 小结 个人觉得选频界面有上面两个就足够了，当然也可以根据需求自定义一些界面，如典型链路、电离层参数展示等界面。 4 后记 短波通信选频很多时候可能并不能够满足用户的所有需求，还需要一些对短波通信环境判断的辅助手段，如根据太阳活动、地磁活动等数据的反应，分析处理成一份有一定时效性和简单易懂结论的报告或指数。之前也参与了一部分相关内容的项目，也有一些思考和想法，以后有时间再写一写。 个人觉得研究应该是带着解决某个问题或疑问的目的去进行的，最终结果应该是实实在在的为某个问题的解决提供了帮助，或者解开了心中对未知领域的部分疑惑。过于追求新颖、独特，带来的未必就是有用的、有效的、实在的内容。 记得初中有一次数学课，老师讲解证明两个直角三角形全等的问题（已知斜边和直角边相等），给了两种证明方法。激发了不少同学提出新解法的兴趣，然后大家一起提出了好几种证明过程复杂且冗余的方法（我贡献了两三个= =），比如先用勾股定理计算另外一条直角边，再通过arctan计算斜边旁边的一个角，然后通过SAS证明全等。作为一时思考的乐趣固然很有意思，但从实际角度来说，HL定理已经是前人总结后最简单有效的方法了。

【GFZ】地磁指数与太阳指数 0 前言 介绍如何下载德国地学中心GFZ提供的地磁行星3小时Kp指数、相关地磁指数和太阳指数。 1 GFZ-FTP 1.1 Kp_ap_Ap_SN_F107 Python爬虫源码如下 import os import socket from ftplib import FTP from ftplib import error_perm ##----------------------------------------------------------------------## # INFO: ftp网站连接 ##----------------------------------------------------------------------## # Outputs: # ftp - ftp根目录 ##----------------------------------------------------------------------## # author: Washy # date: 2024/09/27 ##----------------------------------------------------------------------## def ftp_connect(): # FTP站点 host = 'ftp.gfz-potsdam.de' # 端口号 port = 21 # 文件夹路径 folderpath = '/pub/home/obs/Kp_ap_Ap_SN_F107' ftp = FTP() ftp.encoding = 'utf-8' try: ftp.connect(host, port) ftp.login() ftp.cwd(folderpath) except(socket.error, socket.gaierror): # ftp 连接错误 print("ERROR: cannot connect [{}:{}]".format(host, port)) return None except error_perm: # 用户登录认证错误 print("ERROR: user Authentication failed ") return None return ftp def is_ftp_file(ftp_conn, filename): try: if filename in ftp_conn.nlst(os.path.dirname(filename)): return True else: return False except error_perm: return False # 连接服务器 ftp = ftp_connect() # 需要下载的文件名 filename = "Kp_ap_Ap_SN_F107_2024.txt" # 下载文件 with open(filename, 'wb') as f: ftp.retrbinary('RETR ' + filename, f.write, 1024) # 断开服务器 ftp.close() 下载的数据格式示例如下 # PURPOSE: This file distributes the geomagnetic planetary three-hour index Kp and associated geomagnetic indices as well as relevant solar indices. # LICENSE: CC BY 4.0, except for the sunspot numbers contained in this file, which have the CC BY-NC 4.0 license # SOURCE: Geomagnetic Observatory Niemegk, GFZ German Research Centre for Geosciences # PLEASE CITE: Matzka, J., Stolle, C., Yamazaki, Y., Bronkalla, O. and Morschhauser, A., 2021. The geomagnetic Kp index # and derived indices of geomagnetic activity. Space Weather, https://doi.org/10.1029/2020SW002641 # # Kp, ap and Ap # The three-hourly equivalent planetary amplitude ap is derived from Kp and the daily equivalent planetary amplitude Ap is the daily mean of ap. # Kp is unitless. Ap and ap are unitless and can be multiplied by 2 nT to yield the average geomagnetic disturbance at 50 degree geomagnetic latitude. # Kp, ap and Ap were introduced by Bartels (1949, 1957) and are produced by Geomagnetic Observatory Niemegk, GFZ German Research Centre for Geosciences. # Described in: Matzka et al. (2021), see reference above. # Data publication: Matzka, J., Bronkalla, O., Tornow, K., Elger, K. and Stolle, C., 2021. Geomagnetic Kp index. V. 1.0. GFZ Data Services, # https://doi.org/10.5880/Kp.0001 # Note: the most recent values are nowcast values and will be replaced by definitive values as soon as they become available. # # International Sunspot Number SN # The international sunspot number SN (written with subscript N) is given as the daily total sunspot number version 2.0 introduced in 2015. # The sunspot data is available under the licence CC BY-NC 4.0 from WDC-SILSO, Royal Observatory of Belgium, Brussels. Described in: # Clette, F., Lefevre, L., 2016. The New Sunspot Number: assembling all corrections. Solar Physics, 291, https://doi.org/10.1007/s11207-016-1014-y # Note: the most recent values are preliminary and replaced by definitive values as soon as they become available. # # F10.7 Solar Radio Flux # Local noon-time observed (F10.7obs) and adjusted (F10.7adj) solar radio flux F10.7 in s.f.u. (10^-22 W m^-2 Hz^-1) is provided by # Dominion Radio Astrophysical Observatory and Natural Resources Canada. # Described in: Tapping, K.F., 2013. The 10.7 cm solar radio flux (F10.7). Space Weather, 11, 394-406, https://doi.org/10.1002/swe.20064 # Note: For ionospheric and atmospheric studies the use of F10.7obs is recommended. # # Short file description (for a detailed file description, see Kp_ap_Ap_SN_F107_format.txt): # 40 header lines, all starting with # # ASCII, blank separated and fixed length, missing data indicated by -1.000 for Kp, -1 for ap and SN, -1.0 for F10.7 # YYYY MM DD is date of UT day, days is days since 1932-01-01 00:00 UT to start of UT day, days_m is days since 1932-01-01 00:00 UT to midday of UT day # BSR is Bartels solar rotation number, dB is day within BSR # Kp1 to Kp8 (Kp for the eight eighth of the UT day), ap1 to ap8 (ap for the eight eighth of the UT day), Ap, SN, F10.7obs, F10.7adj # D indicates if the Kp and SN values are definitive or preliminary. D=0: Kp and SN preliminary; D=1: Kp definitive, SN preliminary; D=2 Kp and SN definitive # # # The format for each line is (i stands for integer, f for float): #iii ii ii iiiii fffff.f iiii ii ff.fff ff.fff ff.fff ff.fff ff.fff ff.fff ff.fff ff.fff iiii iiii iiii iiii iiii iiii iiii iiii iiii iii ffffff.f ffffff.f i # The parameters in each line are: #YYY MM DD days days_m Bsr dB Kp1 Kp2 Kp3 Kp4 Kp5 Kp6 Kp7 Kp8 ap1 ap2 ap3 ap4 ap5 ap6 ap7 ap8 Ap SN F10.7obs F10.7adj D 2024 01 01 33603 33603.5 2596 25 0.667 0.333 0.667 1.333 2.000 3.000 3.333 4.000 3 2 3 5 7 15 18 27 10 54 135.7 131.2 1 2024 01 02 33604 33604.5 2596 26 2.667 2.333 0.667 2.000 2.333 2.667 2.000 0.667 12 9 3 7 9 12 7 3 8 66 142.1 137.4 1 2024 01 03 33605 33605.5 2596 27 2.667 2.667 1.000 1.667 1.667 2.667 3.333 3.000 12 12 4 6 6 12 18 15 11 57 140.2 135.5 1 2024 01 04 33606 33606.5 2597 1 1.333 0.667 1.667 0.667 0.667 1.333 1.000 2.667 5 3 6 3 3 5 4 12 5 98 125.8 121.6 1 2024 01 05 33607 33607.5 2597 2 2.000 1.667 0.333 1.000 1.333 1.667 0.333 0.333 7 6 2 4 5 6 2 2 4 117 152.7 147.6 1

K-means聚类算法原理 0 前言 聚类是数据处理中常用的分析方法，此处简单介绍下K-means聚类算法原理。 1 K-means算法 1.1 算法简介 K-means标准算法是1957 年史都华·劳埃德(Stuart Lloyd)作为一种脉冲码调制的技术所提出，但直到1982年才被贝尔实验室公开出版在“ IEEE Transactions on Information Theory”中，原始文献为“Least square quantization in PCM”。术语“k-均值”于1967年才被詹姆斯·麦昆(James MacQueen) 在文献“Some methods for classification and analysis of multivariate observations”中提出，描述 K-means算法的完整理论并进行了详细的研究。 1.2 算法思想 K-means算法，又称K均值算法，其中K表示聚类簇数，means表示取每个簇所有数据的均值作为该簇的中心（或者称质心）。K-means算法的核心思想是将数据集中的n个对象划分为K个聚类，使得每个对象到其所属聚类的中心（或称为均值点、质心）的距离之和最小。这里所说的距离通常指的是欧氏距离，但也可以是其他类型的距离度量。 1.3 算法流程 对于待分类的数据X，随机选取K个簇中心 计算各个数据与各个簇中心的“距离”，并将各个数据划分至相距最近的簇 根据划分好的簇，计算每个簇的数据均值作为新的簇中心 重新计算所有数据与新的簇中心的“距离”，并重新分类 重复上述步骤，直至收敛 1.4 数学表达 假设存在一系列散点$X_i = (x_i,y_i)$，需要将其划分至$K$个簇 随机选择$K$个簇中心$(x_j,y_j)$，其中$j=1,2,\dots,K$ 计算$X_i$到各个簇中心的欧式距离（也可以是其他距离，此处以欧式距离为例），并将$X_i$划分至$r_{i,j}$最小的簇$S_k$ $$ r_{i,j} = \sqrt{(x_i-x_j)^2 + (y_i - y_j)^2}, \quad j=1,2,\dots,K. $$ 重新计算每个簇的均值作为新的簇中心 $$ x_k = \sum_{x_i \in S_k} x_i / N(S_k), \quad y_k = \sum_{y_i \in S_k} y_i / N(S_k), \quad k = 1,2,\dots,K. $$ 其中$S_k$表示第$k$个簇所包含的数据，$N(S_k)$表示第$k$个簇所包含的数据个数。 重复上述步骤，直至新的簇中心与旧的簇中心没有差异 1.5 实例演示 import numpy as np import matplotlib.pyplot as plt from sklearn.datasets import make_classification # 定义数据集 X, _ = make_classification(n_samples=1000, n_features=2, n_informative=1, n_redundant=0, n_clusters_per_class=1, random_state=1) # 根据欧氏距离将数据X分成两类 def rho_p1_p2(X,p1,p2): nums = X.shape[0] l1 = [] l2 = [] for idx in np.arange(nums): rho1 = np.linalg.norm(X[idx,:] - p1) rho2 = np.linalg.norm(X[idx,:] - p2) if rho1<rho2: l1.append(idx) else: l2.append(idx) return l1,l2 # 随机中心一 p1 = np.array([-2,-2]) # 随机中心二 p2 = np.array([2,2]) # 绘图 plt.figure(figsize=[15,8],dpi=300) # 第一幅图为原始数据 plt.subplot(2,3,1) plt.scatter(X[:,0],X[:,1]) for idx in np.arange(20): # 根据当前簇中心进行分类 l1,l2 = rho_p1_p2(X,p1,p2) # 计算新的簇中心 tmp1 = np.mean(X[l1,:],0) tmp2 = np.mean(X[l2,:],0) # 每进行一次分类绘制一次分类后的图像 plt.subplot(2,3,idx+2) plt.scatter(X[l1,0],X[l1,1]) plt.scatter(X[l2,0],X[l2,1]) # 判断中心是否发生变化 if np.linalg.norm(tmp1-p1)<1e-6 and np.linalg.norm(tmp2-p2)<1e-6: print(idx) break # 将簇中心替换为新的簇中心 p1 = tmp1 p2 = tmp2 plt.show() 可以看出随着迭代此次数的增加，数据逐渐被分为两类。 2 拓展 2.1 “距离”计算方法 曼哈顿距离（Manhattan Distance） $$ d(x,y) = \sum_{i=1}^n \vert x_i - y_i \vert $$ 欧几里得距离（Euclidean Distance） $$ d(x,y) = \sqrt{\sum_{i=1}^n (x_i - y_i)^2} $$ 切比雪夫距离（Chebyshev Distance） 切比雪夫距离起源于国际象棋中国王的走法，国际象棋中国王每次只能往周围的8格中走一步，那么如果要从棋盘中A格(x1,y1)走到B格(x2,y2)最少需要走几步？你会发现最少步数总是max(|x2-x1|,|y2-y1|)步。 $$ d(x,y) = \max \vert x_i - y_i \vert $$ 闵氏距离（Minkowski Distance） 对于点$x=(x_1,x_2,\dots,x_n)$ 与点$y=(y_1,y_2,\dots,y_n)$，闵氏距离可以用下式表示： $$ d(x,y) = \left( \sum_{i=1}^n \vert x_i - y_i \vert^p \right)^{1/p} $$ 闵氏距离是对多个距离度量公式的概括性的表述，p=1退化为曼哈顿距离；p=2退化为欧氏距离；切比雪夫距离是闵氏距离取极限的形式。 参考 K均值（K-means）聚类算法（Python3实现代码） python机器学习 | 聚类算法之K-Means算法介绍及实现 【机器学习-14】K-means聚类算法：原理、应用与优化 全面归纳距离和相似度方法(7种)

2024-05-27 X2.8级太阳耀斑 北京时间2024年5月27日14:53，太阳爆发了一个X2.8级X射线耀斑，于15:24结束，如下展示的是NASA 304A波段的太阳观测视频。 参考 王者归来 原活动区AR3664再次爆发X2.8级耀斑

Savitzky-Golay滤波器原理-01 0 前言 最近在看文章时，发现有人使用Savitzky-Golay滤波器对数据进行预处理，提取背景值，尝试之后发现效果显著。因此，简单研究下Savitzky-Golay滤波器的原理。 1 Savitzky-Golay滤波器原理 1.1 简介 Savitzky-Golay滤波器（通常简称为S-G滤波器）最初由Savitzky和Golay于1964年提出，发表于Analytical Chemistry 杂志。之后被广泛地运用于数据流平滑除噪，是一种在时域内基于局域多项式最小二乘法拟合的滤波方法。这种滤波器最大的特点在于在滤除噪声的同时可以确保信号的形状、宽度不变$^{[1]}$。 Savitzy-Golay 卷积平滑算法是移动平滑算法的改进。用 Savitzky-Golay 方法进行平滑滤波，可以提高光谱的平滑性，并降低噪音的干扰。Savitzy-Golay 平滑滤波的效果，随着选取窗宽不同而不同，可以满足不同场合的需求$^{[1]}$。 1.2 推导 思想：使用一定长度窗口的数据对窗口中心的数据进行平滑处理。 设窗口的宽度为$2m + 1$，窗口内的数据值为$y_i$， $i = -m,-m+1,\dots,-1,0,1,\dots,m-1,m$，使用$n$阶多项式对其进行拟合，有 $$ f_i = f(i) = \sum_{k=0}^{n} b_k i^k $$ 那么拟合数据与原始数据的残差平方和为 $$ R = \sum_{i=-m}^{m} (f_i - y_i)^2 = \sum_{i=-m}^{m} \left( \sum_{k=0}^{n} b_k i^k - y_i \right)^2 $$ 使用最小二乘的思想，为了使拟合效果最好，需要使残差的平方和最小，则需要使得所有拟合系数的一阶偏导为零。对于第$l$个系数$b_l$有 $$ \frac{\partial R}{\partial b_l} = 2 \sum_{i=-m}^{m} \left( \sum_{k=0}^{n} b_k i^k - y_i \right)i^l = 0 $$ 整理后可得 $$ \sum_{i=-m}^{m} y_i i^l = \sum_{k=0}^{n} b_k \sum_{i=-m}^{m} i^{k+l} $$ 对于$l \in [0,n]$可以得到$n+1$个方程，联立后可以解得所有系数$b_l$。写成矩阵形式，有 $$ \begin{bmatrix} (-m)^0 & \dots & 0^0 & \dots & m^0 \\ \vdots & & \vdots & & \vdots \\ (-m)^l & \dots & 0^l & \dots & m^l \\ \vdots & & \vdots & & \vdots \\ (-m)^n & \dots & 0^n & \dots & m^n \\ \end{bmatrix} \begin{bmatrix} (-m)^0 & \dots & (-m)^k & \dots & (-m)^n \\ \vdots & & \vdots & & \vdots \\ 0^0 & \dots & 0^k & \dots & m^n \\ \vdots & & \vdots & & \vdots \\ m^0 & \dots & m^k & \dots & m^n \\ \end{bmatrix} \begin{bmatrix} b_0 \\ \vdots \\ b_k \\ \vdots \\ b_n \\ \end{bmatrix} = \begin{bmatrix} (-m)^0 & \dots & 0^0 & \dots & m^0 \\ \vdots & & \vdots & & \vdots \\ (-m)^l & \dots & 0^l & \dots & m^l \\ \vdots & & \vdots & & \vdots \\ (-m)^n & \dots & 0^n & \dots & m^n \\ \end{bmatrix} \begin{bmatrix} y_{-m} \\ \vdots \\ y_0 \\ \vdots \\ y_m \\ \end{bmatrix} $$ 得到所有系数后，将其代回拟合方程，并计算中心点处的拟合值 $$ f_0 = f(0) = b_0 $$ 此时，我们就得到了Savitzky-Golay滤波器平滑后的中心点的数值大小。 1.3 一阶形式 当使用一阶多项式进行拟合时，代入上一节矩阵形式，很容易可以得到 $$ \begin{bmatrix} 1 & \dots & 1 & \dots & 1 \\ -m & \dots & 0 & \dots & m \\ \end{bmatrix} \begin{bmatrix} 1 & -m \\ \vdots & \vdots \\ 1 & 0 \\ \vdots & \vdots \\ 1 & m \\ \end{bmatrix} \begin{bmatrix} b_0 \\ b_1 \\ \end{bmatrix} = \begin{bmatrix} 1 & \dots & 1 & \dots & 1 \\ -m & \dots & 0 & \dots & m \\ \end{bmatrix} \begin{bmatrix} y_{-m} \\ \vdots \\ y_0 \\ \vdots \\ y_m \\ \end{bmatrix} $$ 整理后可得 $$ \left\{ \begin{split} &b_0 = \sum_{i=-m}^m y_i / \sum_{i=-m}^{m} 1 \\ &b_1 = \sum_{i=-m}^m y_i i / \sum_{i=-m}^{m} i^2 \\ \end{split} \right. $$ 从而可以得到中心点处的拟合值为 $$ f_0 = b_0 = \sum_{i=-m}^m y_i / (2m+1) $$ 此即为滑窗平均的表达式。 参考： 【UWB】Savitzky Golay filter SG滤波器原理讲解 Savitzky-Golay滤波器原理阐述

空间物理数据材料整理 0 前言 简单汇总下空间物理学领域部分数据、材料相关的FTP、网站等。 1 数据汇总 序号 简介 链接 备注 1 地磁与太阳活动指数 NOAA-FTP | GFZ-FTP | WDC 爬虫 2 ACE卫星数据 NOAA-FTP 3 GIM Map GIPP-FTP 4 GOLD数据，地球同步轨道远紫外成像仪 GOLD 5 ESA Swarm卫星数据 Swarm 6 NASA太阳观测图像 NASA 7 SEPC太阳H$\alpha$图像 SEPC 8 测高仪数据 GIRO 2 材料汇总 序号 简介 网址 备注 1 国际电联通信相关文档 ITU-R 参考材料 地磁活动指数与太阳活动指数 GFZ数据下载的一种方式分享

史瓦西黑洞最内稳定圆轨道计算 0 前言 最内稳定圆轨道（Innermost Stable Circular Orbit，ISCO）是测试粒子可以稳定地绕广义相对论中的大质量物体运行的最小边缘稳定圆轨道。ISCO的半径$r_{\rm isco}$取决于中心物体的质量和角动量（自旋），它标记了黑洞吸积盘的内边缘。 ISCO 不应与罗希极限混淆，罗希极限是物理物体在潮汐力将其破坏之前可以绕轨道运行的最内点。ISCO 关注的是理论测试粒子，而不是真实物体。一般来说，ISCO 将比罗氏极限更接近中心物体。 本文介绍下如何推导史瓦西黑洞的最内稳定圆轨道，是一篇学习笔记。 1 有效势 在自然坐标系下（$c=G=1$），史瓦西度规的球坐标形式如下 $$ d s^2 = -f dt^2 + \frac{d r^2}{f} + r^2 d \theta^2 + r^2 \sin^2 \theta d \phi \tag{1} $$ 其中$f = 1 - 2M/r$，$M$表示黑洞质量。因此度规张量形式为 $$ g_{\mu\nu} = \begin{bmatrix} -f & 0 & 0 & 0 \\ 0 & f^{-1} & 0 & 0 \\ 0 & 0 & r^2 & 0 \\ 0 & 0 & 0 & r^2 \sin^2 \theta \end{bmatrix} $$ 考虑黑洞周围一个质量为$m$的测试粒子，其拉格朗日作用量为 $$ \widetilde{L} = \frac{m}{2} g_{\mu \nu} \dot{ x }^\mu \dot{ x}^\nu = \frac{m}{2} \left( -f \dot{t}^2 + f^{-1} \dot{r}^2 + r^2 \dot{\theta}^2 + r^2 \sin^2 \theta \dot{\phi}^2 \right) $$ 为了计算方便，选取坐标系使得测试粒子的轨道面为赤道面，即$\theta = \pi/2$，则拉格朗日作用量简化为 $$ \widetilde{L} = \frac{m}{2} \left( -f \dot{t}^2 + f^{-1} \dot{r}^2 + r^2 \dot{\phi}^2 \right) \tag{2} $$ 根据四动量的定义$P_\mu = \partial \widetilde{L} / \partial \dot{ x}_\mu$，可知简化后的形式如下 $$ \begin{split} &P_t = \frac{\partial \widetilde{L}}{\partial \dot{t}} = -m f \dot{t} \equiv -E, \\ &P_r = \frac{\partial \widetilde{L}}{\partial \dot{r}} = mf^{-1} \dot{r}, \\ &P_\phi = \frac{\partial \widetilde{L}}{\partial \dot{\phi}} = m r^2 \dot{\phi} \equiv L. \end{split} \tag{3} $$ 其中$E$表示能量，$L$表示角动量，是两个守恒量。测试粒子的哈密顿量（Hamiltonian）为 $$ \begin{split} \widetilde{H} &= \dot{ x}^\mu P_\mu - \widetilde{L} = \frac{m}{2} g_{\mu \nu} \dot{ x}^\mu \dot{ x}^\nu \\ &= \frac{m}{2} \left( -f \dot{t}^2 + f^{-1} \dot{r}^2 + r^2 \dot{\phi}^2 \right) \end{split} \tag{4} $$ 由$\widetilde{H} = -mk/2$可知 $$ \dot{r}^2 = f(f\dot{t}^2 - r^2 \dot{\phi}^2 - k) = f \left( \frac{E^2}{m^2 f} - \frac{L^2}{m^2 r^2} - k \right) \tag{5} $$ 对上式进行整理，写为“动能”与“势能”之和的形式，即 $$ m^2 \dot{r}^2 + V(r) = E^2 \tag{6} $$ 其中$V(r)$为有效势，形式如下 $$ V(r) = f \left( \frac{L^2}{r^2} + m^2 k \right) \tag{7} $$ 2 圆轨道 测试粒子的运动是圆轨道的条件为 $$ \dot{r} = 0 \quad {\rm and} \quad \ddot{r} = 0 \tag{8} $$ 根据定义式$\ddot{r} \equiv d \dot{r} / dt$以及$\dot{r} \equiv dr/dt$可知 $$ \ddot{r} = \frac{d \dot{r}}{d t} = \frac{d \dot{r}}{dr} \frac{dr}{dt} = \dot{r} \frac{d \dot{r}}{dr} = \frac{1}{2} \frac{d \dot{r}^2}{dr} \tag{9} $$ 联立(5)式和(9)式可得 $$ \ddot{r} = 2ff'\dot{t}^2 - r(2f+rf') \dot{\phi}^2 - kf' \tag{10} $$ 将(5)式和(10)式代入(8)式可得 $$ \dot{t}^2 = \frac{2k}{2f - rf'}, \qquad \dot{\phi}^2 = \frac{kf'}{r(2f - rf')}. \tag{11} $$ 联立(3)式可知 $$ E^2 = \frac{2m^2kf^2}{2f - rf'}, \qquad L^2 = \frac{m^2 k r^3 f'}{2f-rf'}. \tag{12} $$ 3 最内稳定圆轨道 对于最内稳定圆轨道，除了圆轨道条件外，还需要满足有效势的二阶导数为零，即 $$ V''(r) = 0 \tag{13} $$ 首先，对$V(r)$求一阶导数，可得 $$ V'(r) = \frac{dV(r)}{dr} = f' \left( \frac{L^2}{r^2} + m^2k \right) + f \left( \frac{2LL'}{r^2} - \frac{2L^2}{r^3} \right) $$ 根据公式（3）的第三个公式可知$L' = 2L/r$，代入上式可得 $$ V'(r) = f' \left( \frac{L^2}{r^2} + m^2k \right) + \frac{2fL^2}{r^3} = \frac{2f + rf'}{r^3}L^2 + m^2kf' $$ 继续求导，有 $$ \begin{split} V''(r) &= -\frac{3(2f + rf')L^2}{r^4} + \frac{(2f' + f' + rf'')L^2 + (2f+rf')2LL'}{r^3} + m^2kf'' \\ &= (2f + 4rf' + r^2 f'') \frac{L^2}{r^4} + m^2kf'' \\ \end{split} \tag{14} $$ 将(12)式和(14)式代入(13)式，可得最内稳定圆轨道的大小为 $$ r_{\rm isco} = 6M $$ 参考 Innermost stable circular orbit https://ned.ipac.caltech.edu/level5/March01/Carroll3/Carroll7.html

高频传播特性的等效关系 0 前言 短波通信是一种常见的通信方式，其利用高频（High Frequency，HF）电磁波在电离层和地面之间的一次或多次反射，最远可将信号传至上万公里。在这个过程中，可以认为电离层和地面形成了类似“波导管”的结构。 在短波通信中，为了简化计算，通常会使用一些等效关系，这里将根据参考文献[1]中的内容，对这些等效关系进行介绍。 1 平面地球和平面电离层 1.1 “割线定律” - The "Secant Law" 在本节中，我们将考虑两个波的频率、虚拟路径和吸收之间存在的关系，一个波以斜入射反射，另一个波以垂直入射反射，二者的反射发生在相同的真实高度。为此，考虑以如下图所示的角度入射到平面电离层上的射线，其中电子密度随高度增加，从而发生全内反射。 Refs [1] 图4.1 平面地球与平面电离层的等效定理 在没有碰撞和外加磁场的情况下，在等离子体频率为$f_N$的水平高度上，频率为$f$的波的折射率$\mu$由下式给出： $$ \mu^2 = 1 - \left( \frac{f_N}{f} \right)^2 \tag{1} $$ 在反射层应用斯涅尔定律（Snell's law），即$\mu = \sin \phi_0$，可知 $$ f_N = f \cos \phi_0 $$ 如果$f_\nu$是与斜入射频率$f$在相同的真实高度（即相同的等离子体频率）上以垂直入射反射的频率，则$f_\nu =f_N$。因此 $$ f_\nu = f \cos \phi_0 \qquad {\rm or} \qquad f = f_\nu \sec \phi_0 \tag{2} $$ 该频率称为“等效垂直入射频率”，对应于$f$。上式中第二个式子就是所谓的割线定律。上述结果表明，在斜入射下，电离层反射的频率比正常入射时高得多。 1.2 Breit和Tuve定理 - Breit and Tuve's Theorem 另一个比较重要的关系称为Breit和Tuve定理。其中表示发射机T与接收机R之间传输的群（或等效）路径$P'$由等效三角形TAR的长度给出，即： $$ P' = TA + AR \tag{3} $$ 这可以通过以下论证来证明： $$ P' = \int_{TBR} \frac{ds}{\mu} = \int \frac{dx}{\mu \sin \phi} = \frac{1}{\sin \phi_0} \int dx = \frac{TR}{\sin \phi_0} = TA + AR $$ 请注意，反射的真实高度B总是小于A处的等效高度。 应该记住，只有当发送者和接收者位于电离层之外时（3）式才成立。如果发射机和接收机位于电离层内，即折射率为$\mu_1$的水平高度上，则（3）式的右侧必须除以$\mu_1$，这样$P'$仍然意味着群传播时间乘以自由空间速度。 1.3 马丁（等效路径）定理 - Martyn's (Equivalent Path) Theorem 如果$f$和$f_\nu$是从同一实际高度斜向和垂直反射的波的频率，则垂直反射信号的虚高等于斜向信号的等效三角路径的高度。 考虑等离子体频率为$f_N$的相同实际高度下，斜波和垂直波的折射率为$\mu_{ob}$和$\mu_\nu$，我们有 $$ \mu_{ob}^2 = 1 - \left( \frac{f_N}{f} \right)^2 \qquad {\rm and} \qquad \mu_\nu^2 = 1 - \left( \frac{f_N}{f \cos \phi_0} \right)^2 \tag{4} $$ 根据斯涅尔定律（Snell's law）可知$\mu_{ob} \sin \phi = \sin \phi_0$。联立这些方程可以得到 $$ \mu_{ob} \cos \phi = \mu_\nu \cos \phi_0 \tag{5} $$ 斜测信号的群路径为 $$ P' = \int_{TBR} \frac{ds}{\mu_{ob}} = 2 \int_0^{h_r} \frac{dh}{\mu_{ob} \cos \phi} = \frac{2}{\cos \phi_0} \int_0^{h_r} \frac{dh}{\mu_\nu} = \frac{2}{\cos \phi_0} h_\nu' $$ 其中 $$ h_\nu' = \frac{1}{2} P' \cos \phi_0 = \frac{1}{2} (TA + AR) \cos \phi_0 = AD $$ 因此 $$ P'(f) = h'(f_\nu) \sec \phi_0 \tag{6} $$ 这个定理表达了斜入射波的虚反射高度与等效垂直波的虚反射高度相等的重要关系。 1.4 马丁（吸收）定理 - Martyn's (Absorption) Theorem 略。 $$ \tag{7} $$ $$ \tag{8} $$ 2 电离层曲率的影响 Refs [1] 图4.2 斜入射射线几何图像 在弯曲的电离层中，斯涅耳定律的形式是 $$ \mu r \sin i = \mu_0 r_0 \sin i_0 \tag{9} $$ 式中$r$为从地心到折射率为$\mu$处的半径矢量的长度，半径矢量与光线夹角为$i$，如上图所示，图中$\mu_0$、$r_0$和$i_0$为任意参考值。以地面为参照，则有$\mu_0 = 1$、$r_0 = a$且$i_0 = (\pi/2) - \Delta$，其中$\Delta$为仰角。即 $$ \mu r \sin i = a \cos \Delta \tag{10} $$ 将式（1）代入式（10），使用$f_\nu$替换$f_N$，给出了频率$f$与相同实际高度$h_r$反射的等效垂直频率$f_\nu$之间的关系。即 $$ \left( \frac{f_\nu}{f} \right)^2 = 1 - \left( \frac{a \cos \Delta}{a + h_r} \right)^2 \tag{11} $$ 设$\phi_r'$为未折射光线的延线与半径矢量在$h_r$处的夹角，如上图所示。从几何关系中可以得到 $$ (a + h_r) \sin \phi_r = a \cos \Delta \tag{12} $$ 因此 $$ f_\nu = f \cos \phi_r \tag{13} $$ 由式（11）和式（13）可以看出，等效频率不仅与$\Delta$有关，而且与反射高度$h_r$有关。 $f_\nu$可以用（6）式或（7）式来定义，而不是用反射高度来定义，但与平面电离层的情况不同，$f_\nu$的值将取决于定义。 参考 [1] Ionospheric radio propagation

使用标势推导朗谬尔波 0 准备工作 引入标势$\varphi$和矢势$\mathbf{A}$，有 $$ \mathbf{E} = -\nabla \varphi, \qquad \mathbf{B} = \nabla \times \mathbf{A} $$ 则麦克斯韦方程组可以写为 $$ \nabla^2 \varphi = - \rho / \varepsilon_0, \qquad \nabla^2 \mathbf{A} - \nabla(\nabla \cdot \mathbf{A}) = - \mu_0 \mathbf{J} + \mu_0 \varepsilon_0 \nabla\frac{\partial \varphi}{\partial t} $$ ## 1 求解色散关系 1.1 双流程方程组 对于非磁化等离子体，双流体方程组的电子部分可以写为 $$ \left\{ \begin{split} &\nabla^2 \varphi = - \rho / \varepsilon_0 \\ &\frac{\partial n_e}{\partial t} + \nabla \cdot (n_e \mathbf{u_e}) = 0 \\ &n_e m_e \frac{\partial \mathbf{u_e}}{\partial t} = n_e e \nabla \varphi - \nabla p_e \end{split} \right. \tag{1} $$ 1.2 微扰法和平面波化 微扰法：离子作为正电荷背景，即$n_i = n_0$；设电子数密度由背景量和扰动量组成，即$n_e = n_0 + n_{e1}$。则电荷密度$\rho =-e n_{e1}$。 平面波化：设所有的扰动具有指数项$e^{i(\mathbf{k} \cdot \mathbf{r} - \omega t)}$，即具有平面波的形式，则有$\partial/\partial t = -i\omega$和$\nabla = i\mathbf{k}$。 因此双流体方程中电子部分可以改写为 $$ \left\{ \begin{split} &-k^2 \varphi = \frac{e}{\varepsilon_0} n_{e1} \\ &-i\omega n_{e1} + i n_0 \mathbf{k} \cdot \mathbf{u_e} = 0 \\ &-i\omega n_0 m_e \mathbf{u_e} = i \mathbf{k} n_0 e \varphi - i\mathbf{k} \gamma_e k_B T_e n_{e1} \end{split} \right. \tag{2} $$ 运动速度$\mathbf{u_e}$只考虑波矢$\mathbf{k}$方向的分量。则上述方程组进一步简化为 $$ \left\{ \begin{split} &-k^2 \varphi = \frac{e}{\varepsilon_0} n_{e1} \\ &-\omega n_{e1} + k n_0 u_{e\parallel} = 0 \\ &-\omega n_0 m_e u_{e \parallel} = k n_0 e \varphi - k \gamma_e k_B T_e n_{e1} \end{split} \right. \tag{3} $$ 联立上述方程组第二、第三两个方程，消去$u_{e\parallel}$可得 $$ n_{e1} = \frac{k^2 n_0 e}{k^2 \gamma_e k_B T_e - \omega^2 m_e} \varphi \tag{4} $$ 将上式代入方程组的第一个方程，消去$n_{e1}$和$\varphi$可得朗谬尔波的色散关系为 $$ \omega^2 = \frac{n_0 e^2}{\varepsilon_0 m_e} + \frac{\gamma_e k^2 k_B T_e}{m_e} = \omega_{pe}^2 + \frac{\gamma_e}{2} k^2 v_{the}^2 \tag{5} $$ 2 总结 从上述推导过程可以看出，使用标势推导色散关系相对而言更为简洁，且结论一致。标势的引入，也更有利于将结论推广至弯曲时空。此处，只是一个简单的尝试和复习以前的知识，并无独创之处。

朗谬尔振荡和朗谬尔波 1 朗谬尔（Langmiur）振荡 1.1 限制条件 考虑温度可以忽略不计的非磁化等离子体，对于与电子相关的现象（即高频部分），离子由于质量较大无法及时响应高频振荡，因此可看作不动的正电荷背景。当电子相对离子发生小扰动时，设扰动产生的扰动电场$\mathbf{E_1}$。 下图中红点表示离子，蓝点表示电子。 忽略磁场$\mathbf{B}$的等离子体称为非磁化等离子体，相对地，考虑磁场$\mathbf{B}$的等离子体称为磁化等离子体。 设$n_{e,i}$表示电子/离子数密度，$\mathbf{u_{e,i}}$表示电子/离子速度，$m_{e,i}$表示电子/离子质量，则电荷密度$\rho = e (n_i - n_e)$。 1.2 双流体方程组 由高斯定律可知 $$ \nabla \cdot \mathbf{E_1} = \frac{\rho}{\varepsilon_0} = \frac{e}{\varepsilon_0}(n_i - n_e) \tag{1} $$ 由法拉第电磁感应定律可知 $$ \nabla \times \mathbf{E_1} = - \frac{\partial \mathbf{B}}{\partial t} = 0 \tag{2} $$ 电子和离子的运动需要分开讨论，由于离子作为不动的正电荷背景，此处只考虑电子。由连续性方程可知 $$ \frac{\partial n_e}{\partial t} + \nabla \cdot (n_e \mathbf{u_e}) = 0 \tag{3} $$ 由运动方程可知 $$ n_e m_e \frac{\partial \mathbf{u_e}}{\partial t} = -n_e e(\mathbf{E_1} + \mathbf{u_e} \times \mathbf{B}) = - n_e e \mathbf{E_1} \tag{4} $$ 1.3 微扰法和平面波化 微扰法：离子作为正电荷背景，即$n_i = n_0$；设电子数密度由背景量和扰动量组成，即$n_e = n_0 + n_{e1}$。代入上述方程组，并忽略二阶小项，则方程组可改写为 $$ \left\{ \begin{split} &\nabla \cdot \mathbf{E_1} = -\frac{e}{\varepsilon_0} n_{e1} \\ &\nabla \times \mathbf{E_1} = 0 \\ &\frac{\partial n_{e1}}{\partial t} + n_0 \nabla \cdot \mathbf{u_e} = 0 \\ &m_e \frac{\partial \mathbf{u_e}}{\partial t} = -e\mathbf{E_1} \end{split} \right. \tag{5} $$ 平面波化：设所有的扰动具有指数项$e^{i(\mathbf{k} \cdot \mathbf{r} - \omega t)}$，即具有平面波的形式。利用$\partial/\partial t = -i\omega$和$\nabla = i\mathbf{k}$对上述方程组进行化简 $$ \left\{ \begin{split} &i\mathbf{k} \cdot \mathbf{E_1} = -\frac{e}{\varepsilon_0} n_{e1} \\ &i\mathbf{k} \times \mathbf{E_1} = 0 \\ &-i\omega n_{e1} + i n_0 \mathbf{k} \cdot \mathbf{u_e} = 0 \\ &-i\omega m_e \mathbf{u_e} = -e\mathbf{E_1} \end{split} \right. \tag{6} $$ 由$i\mathbf{k} \times \mathbf{E_1} = 0$可知$\mathbf{k} \parallel \mathbf{E_1}$，运动速度$\mathbf{u_e}$只考虑扰动电场$\mathbf{E_1}$方向的分量。则上述方程组进一步简化为 $$ \left\{ \begin{split} &i k E_1 = -\frac{e}{\varepsilon_0} n_{e1} \\ &i\omega n_{e1} - i n_0 k u_{e\parallel} = 0 \\ &i\omega m_e u_{e\parallel} = e E_{1} \end{split} \right. \tag{7} $$ 由上述方程组的第二和第三式消去$u_{e\parallel}$可得 $$ n_{e1} = - \frac{i k n_0 e}{\omega^2 m_e} \tag{8} $$ 将(8)式代入(7)式第一个方程，消去$E_1$、$n_{e1}$可得 $$ \omega^2 = \frac{n_0 e^2}{\varepsilon_0 m_e} \Rightarrow \omega_{pe} = \sqrt{\frac{n_0 e^2}{\varepsilon_0 m_e}} \tag{9} $$ 其中$\omega_{pe}$表示朗谬尔频率，也叫做（电子）等离子体频率。由群速度$v_g = d \omega/dk = 0$可知，该振荡只存在于振荡产生的位置，不会向外传播。 2 朗谬尔波 2.1 限制条件 当电子温度不为零且其他条件不变时，由于电子的热运动，可以将振荡区域的信息携带至邻近区域，从而使邻近区域也发生振荡。这样，发生在某处的振荡就能传播出去而形成波，这种波称为等离子体波或朗谬尔波$^{[1]}$。 [1] 《等离子体物理学》李定著，P92. 2.2 色散关系 此时，高斯定律、法拉点电磁感应定律和连续性方程的形式不发生改变，运动方程中需要引入热压梯度项$\nabla p_e$，如下 $$ n_e m_e \frac{\partial \mathbf{u_e}}{\partial t} = -n_e e(\mathbf{E_1} + \mathbf{u_e} \times \mathbf{B}) - \nabla p_e = -n_e e \mathbf{E_1} - \nabla p_e \tag{10} $$ 由状态方程$p_e \rho_e^{-\gamma_e} = {\rm consts}$和$p_e = n_e k_B T_e$以及$\rho_e = -e n_e$，可得 $$ \begin{split} &\nabla(p_e \rho_e^{-\gamma_e}) = \rho_e^{-\gamma_e} \nabla p_e + p_e \nabla \rho_e^{-\gamma_e} = \rho_e^{-\gamma_e} \nabla p_e - p_e \gamma_e \rho_e^{-\gamma_e - 1} \nabla \rho_e = 0 \\ \Rightarrow & \nabla p_e = p_e \gamma_e \rho_e^{-1} \nabla \rho_e = \gamma_e k_B T_e \nabla n_{e1} \end{split} \tag{11} $$ 其中$k_B$表示玻尔兹曼常数，$T_e$表示电子温度，$\gamma_e$表示电子比热比（在绝热假设中，$\gamma_e = 3$；在等温假设中，$\gamma_e = 1$）。将上式代入(10)式，则完整的方程组可以写为 $$ \left\{ \begin{split} &\nabla \cdot \mathbf{E_1} = -\frac{e}{\varepsilon_0} n_{e1} \\ &\nabla \times \mathbf{E_1} = 0 \\ &\frac{\partial n_{e1}}{\partial t} + n_0 \nabla \cdot \mathbf{u_e} = 0 \\ &n_e m_e \frac{\partial \mathbf{u_e}}{\partial t} = -n_e e \mathbf{E_1} - \gamma_e k_B T_e \nabla n_{e1} \end{split} \right. \tag{12} $$ 平面波化后，上述方程组改写为 $$ \left\{ \begin{split} &i\mathbf{k} \cdot \mathbf{E_1} = -\frac{e}{\varepsilon_0} n_{e1} \\ &i\mathbf{k} \times \mathbf{E_1} = 0 \\ &-i\omega n_{e1} + i n_0 \mathbf{k} \cdot \mathbf{u_e} = 0 \\ &-i\omega n_0 m_e \mathbf{u_e} = -n_0 e\mathbf{E_1} - i k \gamma_e k_B T_e n_{e1} \end{split} \right. \tag{13} $$ 同样地，运动速度$\mathbf{u_e}$只考虑扰动电场$\mathbf{E_1}$方向的分量，将上述方程组化简为 $$ \left\{ \begin{split} &i k E_1 = -\frac{e}{\varepsilon_0} n_{e1} \\ &i\omega n_{e1} - i n_0 k u_{e\parallel} = 0 \\ &i\omega n_0 m_e u_{e\parallel} = n_0 e E_{1} + i k \gamma_e k_B T_e n_{e1} \end{split} \right. \tag{14} $$ 由上述方程组的第二和第三式消去$u_{e\parallel}$可得 $$ n_{e1} = - \frac{ik n_0 e}{\omega^2 m_e - k^2 \gamma_e k_B T_e} E_1 \tag{15} $$ 将(15)式代入(14)式第一个方程，消去$E_1$、$n_{e1}$可得 $$ \omega^2 = \frac{n_0 e^2}{\varepsilon_0 m_e} + \frac{\gamma_e k^2 k_B T_e}{m_e} = \omega_{pe}^2 + \frac{\gamma_e}{2} k^2 v_{the}^2 \tag{16} $$ 其中$v_{the} = \sqrt{2 k_B T_e / m_e}$表示电子热速度。上式即为朗谬尔波的色散关系，从(16)式可以看出，只有当$\omega > \omega_{pe}$时，朗谬尔波才能传播。