【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种)

Ubuntu安装gcc及gsl库 1 安装gcc 使用如下命令安装gcc，等待安装成功即可 sudo apt install gcc 若提示报错，大概率是没有更改镜像源，可参考博客Ubuntu修改源进行修改 2 安装gsl 前往GSL官网下载GSL-latest最新版本，可从下面任意链接进入FTP网站 最近的GNU镜像 GNU FTP主站点 将压缩包解压至任意路径 打开终端进入解压路径，使用如下命令进行安装，等待安装成功即可 ./configure sudo make sudo make install 如果./configure步骤提示找不到命令，则是因为权限不够，使用sudo chmod +x configure命令添加权限即可 参考 linux命令行为什么输入sudo ./configure提示找不到命令

Ubuntu修改源 0 前言 记录下Ubuntu使用的一些基本操作。 1 修改源 使用默认源会因为网络问题无法正常更新、安装包等，因此需要修改为国内镜像，国内镜像链接如下（任选一个即可）： # 清华 http://mirrors.tuna.tsinghua.edu.cn/ubuntu/ # 中科大 http://mirrors.ustc.edu.cn/ubuntu/ # 阿里云 http://mirrors.aliyun.com/ubuntu/ # 网易 http://mirrors.163.com/ubuntu/ 1.1 直接修改 使用如下命令打开文件，并将其中所有的链接修改为国内源链接 sudo vi /etc/apt/sources.list 1.2 代码修改 使用如下命令，将【默认源】修改为【清华源】 sudo sed -i 's/archive.ubuntu.com/mirrors.tuna.tsinghua.edu.cn/g' /etc/apt/sources.list 如果源已被修改，只需对应替换链接，如将【清华源】修改为【中科大源】 sudo sed -i 's/mirrors.tuna.tsinghua.edu.cn/mirrors.ustc.edu.cn/g' /etc/apt/sources.list 2 更新包管理器 修改为国内源之后，使用如下命令更新包管理器 sudo apt-get update sudo apt-get upgrade 参考 Ubuntu修改源镜像方法（22.04也能用）附带常用源镜像地址 终端一行命令更换ubuntu国内镜像源

Mac安装gsl库及配置 0 前言 前段时间在Mac上运行C程序，需要调用gsl库，使用过程中遇到一些问题，在网上找了不少博客才解决，在此记录下。 1 安装gsl库 首先需要安装Homebrew和gcc，可参考Mac安装Homebrew和M1芯片Mac安装gcc 使用brew命令安装gsl（2024-06版本为2.7.1） brew install gsl 等待安装完成即可 默认安装路径为/opt/homebrew/Cellar/gsl/2.7.1 2 配置gsl库 进入终端，打开根目录下的.zprofile文件 vim .zprofile 在文件末尾添加头文件路径至C检索目录 export C_INCLUDE_PATH=$C_INCLUDE_PATH:/opt/homebrew/Cellar/gsl/2.7.1/include 在文件末尾添加链接库至检索C/C++目录 export LIBRARY_PATH=$LIBRARY_PATH:/opt/homebrew/Cellar/gsl/2.7.1/lib 保存并关闭文件，然后重新加载 source .zprofile 在终端输入如下命令（cpp-13是因为我安装的gcc版本是13.x），如果能看到/opt/homebrew/Cellar/gsl/2.7.1/include路径，则说明头文件路径已经能够被检索 cpp-13 -v 在终端输入如下命令，返回值一致则说明链接库已经配置成功 # 命令 gsl-config --libs # 返回值 -L/opt/homebrew/Cellar/gsl/2.7.1/lib -lgsl -lgslcblas 3 补充 3.1 clang与GNU GCC M1 Mac默认安装有clang用于编译C程序，因此直接使用gcc调用的是clang，测试如下 # 命令 gcc -v # 返回值 Apple clang version 15.0.0 (clang-1500.3.9.4) Target: arm64-apple-darwin23.4.0 Thread model: posix InstalledDir: /Library/Developer/CommandLineTools/usr/bin 如果想调用GNU GCC，则需要加上版本号，如13.x版本加上-13，测试如下 # 命令 gcc-13 --version # 返回值 gcc-13 (Homebrew GCC 13.2.0) 13.2.0 Copyright (C) 2023 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 确定GNU GCC的版本号，可以用gfortran的版本查看，返回值中13.2.0即为当前安装的GNU GCC版本 # 命令 gfortran --version # 返回值 GNU Fortran (Homebrew GCC 13.2.0) 13.2.0 Copyright (C) 2023 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. clang与GNU GCC在某些命令上存在差异，因此是不能直接互相替代的，比如-fopenmp命令在前者会报错 3.2 其他调用方式 如果不使用第2节的配置，可以使用-I和-L直接指定头文件和链接库路径，只需在编译时加入如下命令即可 -I/opt/homebrew/Cellar/gsl/2.7.1/include -L/opt/homebrew/Cellar/gsl/2.7.1/lib 参考 cmake 添加头文件目录，链接动态、静态库 LINUX中编译C/C++指定头文件和链接库的搜索路径 【C++编译】gcc的-l参数和-L参数 GCC -l选项：手动添加链接库

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

Win10使用msys2安装gsl库 0 前言 最近需要在电脑上跑别人的C程序，发现其中用到了gsl库，记录下如何配置相关环境及安装该库。 1 安装MSYS2 前往MSYS2官网下载msys2-x86_64-20240507.exe软件 按照引导进行安装 安装完成后，在弹出的命令行窗口输入pacman -Syu更新包数据库和基础包 如果不小心关闭了窗口，可双击安装目录下的msys2.exe打开窗口 继续输入pacman -Su更新其余基本软件包 2 配置环境 进入MSYS2安装目录C:\msys64，双击msys2.exe打开命令行窗口 输入pacman -S --needed base-devel mingw-w64-x86_64-toolchain安装mingw-w64 gcc工具 输入pacman -S mingw-w64-x86_64-gcc-fortran安装gfortran工具 将C:\msys64\mingw64\bin添加至系统环境变量 进入【此电脑】，空白处右键，选择【属性】，点击左侧【高级系统设置】，点击右下角【环境变量】 双击系统变量中的【Path】，点击【新建】，写入C:\msys64\mingw64\bin，然后确定即可 3 安装GSL 前往GSL官网下载GSL-latest最新版本，可从下面任意链接进入FTP网站 最近的GNU镜像 GNU FTP主站点 将压缩包解压至任意路径（路径中最好不包含中文） 进入MSYS2安装目录，双击mingw64.exe打开命令行窗口 在命令行中进入上述GSL的解压路径 输入./configure && make && make install，等待安装完毕，此过程较长 安装完成后，可在MSYS2文件夹下的mingw64/bin文件夹中看到libgsl-27.dll（本文安装的GSL为2.7版本） 参考 Windows10安装VScode + mingw64 + GSL 用 VS Code + MSYS 搞定 Windows 上的 Fortran 开发 超详细教程：windows安装MSYS2（mingw && gcc）——更新于2021.11

Win10配置Python环境及安装apexpy包 0 前言 之前一直使用Mac办公，今天因为工作需求需要在Windows上运行之前写的代码，在安装apexpy包时，遇到了不少问题，经过一番折腾，发现是gfortran和Python版本问题导致的，在此记录下。 1 安装Python 目前（2024-05）最新的Python版本为3.12.3，由于apexpy包还没有完全支持最新版本，所以需要使用相对落后的版本，经测试发现3.10.11版本可正常使用 前往Python官网下载3.10.11版本： 64位：Windows installer (64-bit) 32位：Windows installer (32 -bit) 下载完成后，双击安装包，勾选安装界面的两个复选框，点击安装（建议安装在默认路径） 使用win+R快捷键，输入cmd，点击【确定】进入终端 输入python --version，打印Python 3.10.11则安装成功。如果版本号不对，可能是之前安装过python，可选择卸载其他版本或创建新的运行环境（网上自行搜索教程，建议使用必应搜索，非常不推荐使用百度） 输入pip list查看当前已安装的包 2 安装apexpy包 apexpy包需要有gcc和gfortran的运行环境，可参考Win10安装mingw64配置最新版gcc与gfortran环境博客 使用win+R快捷键，输入cmd，点击【确定】进入终端 输入pip install apexpy，等待安装结束即可 该包安装完成后，会顺带安装numpy包 现在的程序还需要使用netCDF4、scipy、matplotlib 同样在终端中使用pip install <package>进行安装，将<package>改为对应包的名称即可 3 安装编辑器 此处推荐使用Visual Studio Code，可前往官网进行下载安装 安装完成后，点击左侧【扩展】图标 搜索Chinese (Simplified) (简体中文)，安装第一个搜索结果 搜索python，安装第一个搜索结果 搜索Jupyter，安装第一个搜索结果 一般来说安装上面三个扩展就能够正常运行代码了，如有其他需求也可自行搜索安装 4 后记 涉及到fortran的东西，总是存在一堆奇怪的问题，害我折腾了这么久。

Win10安装mingw64配置最新版gcc与gfortran环境 0 前言 最近因为多件事情的需求，需要在Windows电脑上配置Fortran环境，由于网上大多数的博客介绍的方法安装的Fortran版本较低，使用过程中会出现各种问题，最终找到了解决办法，在此处记录下。 1 安装MinGW-w64 在线安装版本无法获取最新版本，只能到8.x，最新版本需要从GitHub上获取，链接如下： https://github.com/niXman/mingw-builds-binaries/releases 当前（2024-05）最新版本为13.2，如下图所示。其中i686为32位，x86_64为64位，对于win10电脑，下载图中红框压缩包 将下载好的x86_64-13.2.0-release-win32-seh-ucrt-rt_v11-rev1.7z压缩包解压至C盘根目录，此时路径为C:\mingw64 将C:\mingw64\bin添加至系统环境变量 进入【此电脑】，空白处右键，选择【属性】，点击左侧【高级系统设置】，点击右下角【环境变量】 双击系统变量中的【Path】，点击【新建】，写入C:\mingw64\bin，然后确定即可 使用win+R快捷键，输入cmd，点击【确定】进入终端 在终端中输入gcc --version，正常打印版本号则安装成功 在终端中输入gfortran --version，正常打印版本号则安装成功 如果以上操作都没问题，但版本号不对，可进行如下排查 在终端输入where gfortran，查看打印出来的路径是否存在其他的路径 若存在，可选择删除其他版本的gcc/gfortran，或者将环境变量中其他版本的路径删除，或者将环境变量中新安装的版本路径移动至最上方 注意：gcc与gfortran的版本号必须一致 参考： MinGW-W64 下载、安装与配置（支持最新版的GCC，目前 GCC 13.2.0） 2 后记 这件事也是反复折腾很多次了，每次都会遇到不同的问题。不过这次算是彻底搞明白了，以后直接参考这次结果就行。

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滤波器原理阐述