WH_G_Y的博客

原创 GRACE数据的广义三角帽（TCH）计算不确定度

1、首先整理好需要计算不确定度的GRACE数据，这里以CSR level-2截断至60阶的球谐系数、DDK3和4滤波解、CSR mascon以及JPL mascon数据为例（以JPL mascon数据为参照）[2]姚朝龙, 李琼, 罗志才,等. 利用广义三角帽方法评估GRACE反演中国大陆地区水储量变化的不确定性[J]. 地球物理学报, 2019, 62(3):15.2、根据TCH中的相关公式，利用迭代初值计算噪声协方差矩阵R(n*n)中的相关元素。函数求解在满足约束条件的情况下使得目标函数最小的参数。

原创 所谓真正的告别

所谓真正的告别，没有长亭古道，更没有劝君更进一杯，就是在这样一个正午时分，有些人有些事就都永远停留在回忆里了

原创 实用小工具分享之“cftool”拟合

一般Xdata可以设置为时间变量t，Ydata即为我们要拟合的数据变量（下图红色框）c1、c2：周年振幅项；交流到这个时间序列的拟合问题，感觉这个用了好久的工具确实挺好用的，在这里简单分享出来。1、在matlab命令行窗口输入“cftool”，然后回车，等待片刻弹出工具框。matlab中的cftool小工具，用于序列的各种样式拟合，优点：操作简单、直观、方便、快捷。

原创 计算GRACE质量估计的精度，并绘制某月精度的全球格网图

到目前为止，大多数GRACE研究都集中在制作质量变化图上，很少讨论这些图中的误差。然而，如果GRACE要被用作评估和改进水文和海洋模型的工具，误差估计是必要的。此外，只有误差估计才能确定数据的某些特征是否真实，以及GRACE确定该特征的准确性。上来以示尊重，先给出经典参考文献：Wahr J, Swenson S, Velicogna I. Accuracy of GRACE mass estimates[J]. Geophysical Research Letters, 2006, 33(6).

原创 针对GRACE DDK3滤波解的主成分分析，并展示前四个主成分的时空模态

主成分分析是一种技术，广泛用于数据分析和机器学习中。它的主要目标是通过找到数据的主成分，将高维数据投影到一个低维空间，同时尽可能保留原数据的。PCA 的核心思想是将数据从原始坐标系转换到一个新的坐标系，这个新坐标系由数据的构成。

原创 2025新年的第一天，浅浅做个小结

23和24年共创造博文42篇，获赞279、收藏491和评论114次，博文被分享上千次。创造两年来承蒙大家的持续关注和厚爱，一次次和大家的私信交流，自己收获良多，非常感恩，同时也在一步步朝着自己的目标迈进。在这里立几条2025年flag：1、自己一个人的力量终究是有限的，本年度我想凭借大家的力量，只要大家遇到任何问题，亟待解决的都可以评论或者直接私信我，大家一起to be better；2、本年度创造博文不少于25篇，点赞、收藏和评论分别突破500、1000和300次；3、粉丝数量突破500；

原创 GRACE无约束解和约束解区域TWSA时间序列RMSE计算——以亚马逊流域为例

例如针对无约束解和约束解（选取DDK3滤波解），求亚马逊流域一段时间内水储量的变化，以官方机构JPL和CSR发布的mascon产品为参考，求解均方根误差（RMSE）。：明明从全球格网图上很直观的看出，无约束解具有很严重的条带误差，但就时间序列、相关系数和RMSE可以发现，无约束解和滤波解之间的差距并不大，甚至压根体现不出来差异。：当我们计算整个区域的一些评价指标后发现，滤不滤波的差异不明显的时候，不妨试一试逐个格网点进行求解，说不定会有符合预期的效果。

原创 运用棋盘检验展示GRACE Level-2球谐系数截断至不同阶的影响及其代码分享

本博文旨在针对仿真生成的棋盘格网，分别经过球谐分析得到截断至60、120、200阶数的结果，来清晰地展示截断对最后信号大小的影响。选择美国西部一片矩形区域（经度从西经106°-西经126°，北纬38°-北纬52°） 具体的边界经纬度坐标值存入文后的网盘链接中。是一种简单而有效的图像质量评估方法，它可以帮助用户了解图像处理和显示设备对图像质量的影响，从而选择最佳的图像处理方法和显示设备。最后，将格网转为球谐分别进行200、120和60阶截断，之后再转为区域格网绘制。

原创 鲁棒主成分分析（RPCA）matlab代码实现

主成分分析-维基百科，自由的百科全书 --- Robust principal component analysis - Wikipedia。[CC字幕]流体动力学机器学习：11.鲁棒主成分分析 (RPCA)_哔哩哔哩_bilibili。变量L为低秩矩阵，S为稀疏矩阵，count为迭代次数。

原创 我的创作纪念日

提示：你过去写得最好的一段代码是什么？提示：当前创作和你的工作、学习是什么样的关系。提示：可以和大家分享最初成为创作者的初心。提示：在创作的过程中都有哪些收获。提示：职业规划、创作规划等​​。

原创 在曲线图上最值和极值点位置进行适当标注

1、首先生成一组0-100的随机数，组内共有100个数据；3、在图上分别标注功率谱密度的最值和某些极值。2、求这组数据的功率谱密度，并绘图；%%选出一些较为明显的极值进行标注。

原创 利用matlab bar函数绘制较为复杂的柱状图，并在图中进行适当标注

如何自动选择合适的坐标位置对柱状图的数值大小进行标注？

原创 装载matlab一定要备个灭火器，以免一怒之下把电脑砸了引发火灾烧伤自己

装载matlab一定要备个灭火器，以免一怒之下把电脑砸了引发火灾烧伤自己

原创 Matlab中BaseZoom()函数实现曲线和图片的局部放大

可以先运行如下图片中的语句，看看是否报错；运行下载的BaseZoom工具包下的示例代码“demo_image.m”，操作很简单，鼠标左键选择需要局部放大的区域，蓝色状态时可以调整区域的位置，鼠标右键进行确定。发现生成的figure 1图片中鼠标变为十字丝，代码一直处于运行状态，且命令行窗口显示“运用鼠标左键选择需要放大的矩阵区域”；②在图上较为空旷的位置选中一片矩形区域，区域为蓝色状态时，可以任意调整区域的位置；③选择我们需要局部放在的曲线区域，同样处于蓝色状态时，可以调整位置，点击鼠标右键进行确定；

原创 球谐系数阶方差计算和图谱绘制

1、“GenerateC2030C2021TransIndex”函数 功能：调换球谐系数的排列顺序（例如按阶或者按次升序排列）2、“Tran_vector ”函数 功能：将方阵存储的球谐系数展开为列向量。绘制全球各大洲典型流域的时间序列图-CSDN博客。运用下面这期博文中的格网数据。

原创 ITSG、COST-G、Tongji和WHU Level-2数据产品读取绘图（Matlab）

Institute of Geodesy at Graz University of Technolog（格拉茨理工大学大地测量研究所） 2018版本，最高60阶球谐系数。代码中将所有球谐系数均人为截断至60阶，只对球谐系数做最基本的替换低阶项和扣除平均重力场处理，其余包括GIA改正、滤波处理和泄漏改正等均未涉及；：时变重力场解决方案组合服务解RL01版本，最高96阶球谐系数。：武汉大学 基于卫星间位势差的约束GRACE月重力场模型。：同济大学2022版本，最高96阶球谐系数。

原创 绘制全球各大洲典型流域的时间序列图

以CSR RL05球谐系数为例，生成无约束格网解（grid_CSR）、高斯300km滤波解（grid300_CSR）、去相关滤波解（gridP4M6_CSR）和组合滤波解（grid300P4M6_CSR），并以CSR mascon解（C_mas）作为参考。

原创 Matlab读取Swarm球谐系数，并绘制EWH全球格网图（存在疑问）

从结果可以看出40阶的Swarm球谐系数反演的噪声太大，即使经过强约束的高斯滤波后效果仍然不理想，不大清楚是Swarm数据本身所致还是代码实现有误？还请有懂的友友帮忙解答下，不胜感激！ICGEM官网下载 COST-G发布的40×40的球谐系数。

原创 matlab绘制时间序列图，横坐标轴如何标注为月-日

导入 matlab中，为数值矩阵；Excel表格中有类似于如下 年月日对应的数据。

原创 我的创作纪念日

原创 GRACE滤波数据处理之DDK系列滤波

Kusche 等（2007）从时变重力场模型球谐系数求解出发，对球谐系数的最小二乘解的法方程组应用贝叶斯估计，构造了 DDK 滤波。

原创 标准降水指数（SPI）、标准化降水蒸发蒸腾指数（SPEI）和帕尔默干旱指数（PDSI）的下载、读取和可视化

标准化降水指数通过计算给定时间内降雨量的累计概率,比较客观地表述了多时间尺度下的降水概率,反映降水因素对干旱的影响。Mckee et al.(1993)认为降水量的减少是导致干旱的主要因素,其他气候因子一般变化较小,对干旱产生的影响不大,因此SPI使用降水量作为唯一参数。SPI能较好表达因降水量的大小反映干旱状况,操作简单,所需资料仅为降水量,同时该干旱指数在各个地区和各个时段都具有良好的计算稳定性,是继PDSI之后另一种被广泛应用的干旱指数。值得说明滴是，“Time02_14”为以年为单位的时间变量；

原创 利用Wget批量下载网页链接中的数据

3、可以参照博文下载到指定文件夹；不指定的话，则下载到与文本文档一致的默认文件夹中。wget 的安装与使用（Windows）_wget windows-CSDN博客。2、以Argo温盐数据集Met office（EN4）为例，进行批量下载，wget下载文件到指定目录_wget指定下载目录-CSDN博客。在默认文件夹下放置需要下载的数据链接至txt中，出现此页面 即为下载进行中。

原创 利用Argo数据分别计算温度、盐度和温盐所造成的比容海平面变化

需下载计算海底压力、海水密度以及其他有关海水信息的“seawater”代码包，还需下载Yan-Ning Kuo编写的用于计算比容海平面变化的代码包（“steric_height_calculation.m”）按照上述盐容和热容海平面变化的计算公式，计算热容时，保持盐度为平均盐度不变；3、读取IPRC温盐数据，并绘制比容海平面变化全球趋势图与盐容、热容以及比容在太平洋某处的时间序列图。

原创 一个脚本文件读取CSR、JPL和GSFC的mascon产品，并绘制流域时序图

JPL：https://podaac.jpl.nasa.gov/dataset/TELLUS_GRAC-GRFO_MASCON_CRI_GRID_RL06.1_V3 (可能需要科学上网)由于下载的JPL mascon数据为“.nc4”文件，于是将其修改为“.nc”文件；并将文件名开头修改为“JPL_”，以与其它机构的mascon文件名对应起来。以2002年4月的全球EWH格网点为例，

原创 matlab中 for循环的循环间隔可以根据循环变动么

发现循环间隔并不能随着循环重新赋值，所以到底能不能变动循环间隔呢？欢迎知道的友友评论区或者私信交流。

原创 Windows 11 VS C++程序中安装和使用BLAS线性代数库

​本期是剪剪粘粘 (笑哭)，汇总整理下关于在VS中添加 基础线性代数子程序库（Basic Linear Algebra Subprograms，BLAS）来在window系统中利用相关函数实现线性运算算法库：基础线性代数子程序库（Basic Linear Algebra Subprograms，BLAS）介绍_以下哪个不是常用的基础线性代数子程序库blas库-CSDN博客例如实现如下C++代码中的“cblas_daxpy(i, gamma, &l2[1], -1, &l1[2], 1);”

原创 GRACE/GRACE-FO泄漏误差数据处理之尺度因子法

先用GLDAS水文模型的格网产品经过球谐分析得到球谐系数，再经过截断和高斯滤波或者去相关滤波处理得到滤波后的格网，逐格网点求解滤波前后的差异因子。发现经过尺度因子改正之后，泄漏至海洋的信号明显减弱，但是以长江流域为例，求EWH时间序列发现，乘以GRACE/GRACE-FO经过截断和滤波后的格网数据得到尺度因子改正后的解。3、 将GLDAS滤波前和滤波后的格网按时间序列排列进行最小二乘求尺度因子。2、格网数据转换为球谐（质量），再经过组合滤波得到GLDAS滤波解。尺度因子法的基本原理，直接上图。

原创 结合GLDAS模型估计GRACE/GRACE-FO全球的RMS比率（信噪比）

上述公式出自 Chen J L, Wilson C R, Seo K W. Optimized smoothing of Gravity Recovery and Climate Experiment (GRACE) time‐variable gravity observations[J]. Journal of Geophysical Research: Solid Earth, 2006, 111(B6). 中的公式（6）在JGR: Solid Earth本篇文章网站显示所有引用这篇文章的学者论文。

原创 整理分享GRACE领域可能常用的网站

未完待定，如果大家有更好更实用的网站也欢迎在评论区或者私信分享，后续准备再出一期关于国内外研究卫星重力的机构（大学和研究院等）。数据的获取固然很重要，但是像CSR、JPL和GFZ等网站关于数据的描述也至关重要，它能给我们从文献中无法获取的有关数据的处理和背景知识。近期各种杂事，自己也很久没更新了，想着整理一下GRACE领域相关的比较实用的网站分享给大家，主要还是年末了想着水一期（哈哈哈哈哈哈哈哈）GFZ 网站下的COST-G RL01 (Level-2B/3)(对GRACE的数据处理基本有一个完整的介绍)

原创 从水平衡方程角度验证WGHM模型中GWS分量的大小

最近在进行水文研究时发现，每个水文模型结合水平衡方程计算得到的地下水储量异常（GWSA）时各不相同的，归根结底还是因为每个模型输入水平衡方程的分量种类不一致导致的。简单点说即，A模型TWSA=GWSA+土壤水异常+雪水异常+树冠含水异常，B模型TWSA=GWSA+土壤水异常+雪水异常+树冠含水异常+水库储量异常。（1）利用TWSA -（绿色框内除黄色Groundwater以外的分量）， 绘制全球格网，并绘制区域GWSA时间序列（我们把水平衡方程求得GWSA命名为“GWSA_WB”）对比区域的时间序列图。

原创 利用Arcgis从面shp中得到需要区域的边界点经纬度坐标

1、首先需要选中需要的面区域（如上图），然后右击矢量面区域的文件名-数据-导出数据-修改导出位置和导出shp文件名。2、使用工具中的要素转线功能，将面转为线。

原创 利用GMT6.4及以上版本绘制DEM和地形起伏剖面图

1、利用下载好的等高线数据，调整合适的CPT文件绘制DEM、测点A与B 和测线AB。2、绘制AB线所切得剖面的高程变化情况。2、社区绘图实例 -

原创 VMware 16虚拟机下安装ubuntu 20.04，并调试好GMTSAR

③开机ubuntu的情况下，点击VMware中的虚拟机-继续安装VMware-Tools ，无论vm-tools是否安装成功，都可以尝试按照如下博文教程第一条，打开共享文件夹实现文件互传。①可以不在ubuntu中下载轨道数据ORBITS，网速够快的话，可以点击链接直接在windows中下载，之后在文件互传到虚拟机文件夹中。3、 镜像替换 （一个合适的镜像源能使得ubuntu下安装软件变得非常丝滑，试验很多镜像之后发现下面博文中的阿里云镜像是最适合我的）或者 也可以到github上下载许效华老师的代码包。

原创 matlab图片以合适的比例保存并放置在word中

1、假设绘制得到如下图窗，文件-导出设置-以600pi分辨率导出为jpg格式，图片两侧存在大量空白，需要到word中裁剪（如使用专门裁剪工具还会降低图片的清晰度）2、我们可以调整matlab图窗的大小（如下图所示，等待外侧图框的角点亮起，然后拖拽某一个角点调整图窗整体大小）4、最后再导出设置，可以发现，导出的图片大小与figure图窗大小一致，而且两侧没有空白区域。3、之后再选中内侧坐标区域，使得内部图框的角点亮起，同样的道理调整图框的大小。

原创 如何下载全球重点流域的边界经纬度vec文件

3、解压RBvect文件得到很多vec文件、全球gif图和每个vec名称索引的index文件。Total Runoff Integrating Pathways (TRIP) 网站（2、划到页面大致最下面的位置，找到whole files点击下载。

原创 GRACE泄露改正之全球正向建模恢复法（Forward-Modeling）

提前准备好：（1）CSR机构发布的Mascon产品，这里以CSR-Mascon RL06数据为例，准备了从2002年4月至2021年7月 共197个月的数据（实际试验的时候可以并不需要这么多数据）（2）全球陆地格网点的纬度加权值，即代码中的grid_weight。另一方面，南极区域外非冰盖质量变化信号泄漏至南极区域内，形成内泄漏，两者统称为泄漏误差。①M_obs为Mascon解经过球谐分析和500km高斯滤波后的结果，第一次循环开始前，假设M_tru=M_obs，每次循环的过程中，更新M_tru即可；

原创 GRACE数据的滤波方法实现（Han、Fan和Wiener）

首先以等效水柱高形式绘制测量信号的阶方差图像，再运用上述前两个公式进行最小二乘拟合，求得拟合参数a、b、c、d，最后利用第三公式计算滤波系数h（注意：对于最大阶数60的球谐产品，上述第一个公式适用于0-20阶，第二个公式适用于31-60阶数）示例代码中滤波半斤r0=200km、r1=300km、选定的次m1=15，绘制滤波系数图如下。

原创 绘制60阶无约束解和约束解的球谐系数格网图

①直接从机构数据读取，经过一阶项、二阶项和取平均等预处理操作后的61*121的球谐系数sc矩阵；②例如经过高斯300km之后的等效水柱高，再利用gmt_grid2cs函数转化为球谐系数。2、利用pcolor函数绘制球谐系数格网图。

原创 GRACE数据绘制区域格网图

1、首先下载区域的边界经纬度坐标（如下图展示的为全球主要流域的边界经纬度），解压RBvect文件可以找到index用记事本打开，以编号为0011的Changjiang Yangtze流域边界为例进行说明，2、获取到研究区域边界后，即可参照下面程序进行区域格网图的绘制。出图效果为这样，需要人为手动修改标注条至合理的位置。