找到
5
篇与
空间物理学
相关的结果
-
对短波通信选频的思考 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
-
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滤波器原理阐述
