【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}$时，朗谬尔波才能传播。

非线性最小二乘法拟合函数-3 在前面两个博客中，推导得到了Levenberg-Marquardt算法（简称LM算法）的迭代公式，这里将讲述如何使用Fortran编写一个简单的LM算法。 此处介绍待定系数只有一个的情况。 1 程序设计 1.1 Levenberg_Maquardt_Fit 简介：LM算法子例程，输入离散数据，迭代拟合待定系数 传递参数： n 离散数据的长度 x 离散数据的自变量 y 离散数据的因变量 a 待定系数 1.2 myfunc 简介：待拟合函数的数值计算子例程。为了具有普适性，假定待拟合函数的解析式未知，函数结果只能由此子例程数值计算得到。 传递参数： n 离散数据的长度 x 离散数据的自变量 a 待定系数的当前值 fx 待拟合函数的数值计算结果 1.3 Calculate_Jacobian 简介：计算待拟合函数的Jacobian矩阵的子例程。为了具有普适性，使用差分法计算一阶导数。 传递参数： n 离散数据的长度 x 离散数据的自变量 y 待拟合函数的数值计算结果 a 待定系数的当前值 J 待拟合函数的Jacobian矩阵 2 源代码 2.1 Levenberg_Maquardt_Fit 此处迭代次数上限设置为30次，可根据需求自行更改 对于自变量有多个的情况，可将输入x改为x1,x2,...，满足自变量个数即可 SUBROUTINE Levenberg_Maquardt_Fit(n, x, y, a) IMPLICIT NONE ! 输入参数 INTEGER, INTENT(IN) :: n ! 数据点个数 REAL, DIMENSION(n), INTENT(IN) :: x, y ! 数据点 ! REAL, DIMENSION(3), INTENT(INOUT) :: a ! 待求的拟合系数 REAL, INTENT(INOUT) :: a ! 待求的拟合系数 ! 定义常量 INTEGER, PARAMETER :: m = 1 ! 待求的系数个数 REAL, PARAMETER :: eps = 1.0E-6 ! 收敛阈值 ! 定义变量 REAL :: da(m), r(n), J(n,m), H(m,m) ! 拟合系数、残差、雅可比矩阵、Hessian矩阵 REAL :: lambda, alpha(m) ! 调节因子、步长 INTEGER :: i, iter ! 循环计数器 REAL :: fx(n) ! 初始化调节因子 lambda = 0.001 ! 开始迭代 iter = 0 DO WHILE(iter < 30) ! 迭代次数上限为10000 iter = iter + 1 ! 计算被拟合函数值 CALL myfunc(n,x,a,fx) ! 计算残差向量 r = y - fx ! 计算雅可比矩阵 CALL Calculate_Jacobian(n, x, fx, a, J) ! 计算Hessian矩阵 H = MATMUL(TRANSPOSE(J), J) ! 计算梯度向量 da = MATMUL(TRANSPOSE(J), r) ! 计算搜索方向 DO i = 1, m H(i,i) = H(i,i) + lambda END DO ! CALL SOLVE_LINEAR_SYSTEM(m, H, da, alpha) alpha = da(1)/H(1,1) ! 计算新的拟合系数 a = a + alpha(1) ! 更新调节因子 IF (NORM2(r) < eps) THEN EXIT ELSE IF (NORM2(r) < NORM2(r - MATMUL(J, alpha))/2) THEN lambda = lambda/10. ELSE lambda = lambda*10. END IF END DO END SUBROUTINE 2.2 myfunc 此处使用了一个非常简单的一次函数作为示例，使用时需要将fx = a*x修改为实际的待拟合函数的形式。若待拟合函数是一个子例程，则在此处使用CALL调用 若在Levenberg_Maquardt_Fit子例程中修改了自变量个数，此处需要对应修改 SUBROUTINE myfunc(n,x,a,fx) IMPLICIT NONE ! 输入参数 INTEGER, INTENT(IN) :: n REAL, INTENT(IN) :: x(n) ! 自变量 REAL, INTENT(IN) :: a ! 拟合系数 REAL, INTENT(OUT) :: fx(n) fx = a*x END SUBROUTINE 2.3 Calculate_Jacobian 本子例程使用了差分法计算一阶导数，因此几乎适用于所有的待拟合函数 若在前面修改了自变量的个数，此处需对应修改 SUBROUTINE Calculate_Jacobian(n, x, y, a, J) IMPLICIT NONE ! 输入参数 INTEGER, INTENT(IN) :: n ! 数据点个数 REAL, DIMENSION(n), INTENT(IN) :: x, y ! 数据点 REAL, INTENT(IN) :: a ! 拟合系数 ! 输出参数 REAL, DIMENSION(n,1), INTENT(OUT) :: J ! 雅可比矩阵 ! 定义变量 INTEGER :: i ! 循环计数器 REAL :: dela, y1(n) ! 步进 dela = 0.01 ! 计算 a+dela 对应的 y1 CALL myfunc(n,x,a+dela,y1) ! 计算雅可比 DO i = 1, n J(i,1) = (y1(i) - y(i))/dela END DO END SUBROUTINE 3 测试 编写程序主体，设定y约为x的两倍，调用LM子例程计算待定系数a program name implicit none INTEGER n real :: x(5),y(5),a x = (/1.,2.,3.,4.,5./) y = (/2.1,4.1,5.9,8.1,9.9/) ! y = 2*x**2 + n = 5 a = 1.0 call Levenberg_Maquardt_Fit(n, x, y, a) print *,a end program name 使用gfortran编译生成可执行文件，运行可执行文件后，在终端打印出 1.99818182 离散数据的因变量约为自变量的两倍，因此上述计算结果符合预期