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)中的相关元素。函数求解在满足约束条件的情况下使得目标函数最小的参数。

原创 鲁棒主成分分析（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、获取到研究区域边界后，即可参照下面程序进行区域格网图的绘制。出图效果为这样，需要人为手动修改标注条至合理的位置。

原创 matlab读取GLDAS数据，并绘制全球格网图

本文以GLDAS_Noah格网分辨率为1°×1°的数据为例。3、选择2002年4月至2009年12月的数据为例，get data生成数据下载链接，可以先选中“download instructions”中的“download the list of links .txt”，该文本文件中包含所要下载所有数据的网址链接，后期可以利用IDM软件快速下载。（这件简单解释一下为什么下载2002年4月至2009年12月数据，控制文档中的数据有缺失，因为这是为了和GRACE卫星的数据缺失月份对应上，而人为做的删减操作。

原创 利用Google Earth Pro软件结合Arcgis粗略导出研究区域的经纬度坐标

使用面转线、折点转点 等工具，在点矢量的属性表中添加经度（lon）纬度（lat）字段，并利用计算几何分别计算经纬度坐标，最后点击左上方里面的导出功能，导出为文本文档即可。打开桌面版google earth软件，旋转地球大致定位到自己所要研究的区域，选择软件左上方的“添加多边形”，在地图上依次点击研究区域边界形成一个封闭的多边形

原创 利用m_map代码包生成区域边界的经纬度坐标txt，并绘制图像

以南极洲为例，生成南极洲边界的经纬度坐标文本文档，即“antarctica_boundary.txt”

原创 关于使用word、mathtype和endnote时遇到的一些问题

4、文档中引用公式时，不建议使用Adobe PDF的高质量打印来生成图片高清的PDF，因为生成PDF中公式编号超链接将失效，可以使用Adobe Acrobat软件打开word直接另存为，生成的PDF图像高清而且超链接有用。②将光标置于绿色框内，右键点击"编辑域"，勾选“到段落的超链接”后确定，再切换域代码，之后将“（1）”设置为蓝色字体，转成PDF之后就可以进行超链接跳转啦。可以使用mathtype标签下的“转换公式”按钮，实现公式的统一（不过我发现的一个弊端是，公式中的括号也会因为公式的统一变成斜体）

原创 利用IDM直接抓取网页上的数据，人为选择进行下载

3、为加快下载速度，可以点击左上角此按钮（将选中的文件添加到IDM任务列表和下载列表中）1、在IDM中进行站点抓取，输入网址（需要登录的，则要输入用户名和密码）2、人为选中需要下载的数据，点击左上角三角形状的按钮进行下载。随后你将会看到，大批量选中的文件在下载中。

原创 WaterGAP Global Hydrology Model (WGHM)数据的下载、读取和全球格网图

找到缺失值NaN赋值为0。利用冯伟老师的“GRACE_Matlab_Toolbox-master”工具包绘制全球格网图。%截取1901年的数据。%将数据中的缺省值替换为0。下载对应水文分量的”nc4”文件，以TWS为例，利用如下matlab代码读取。

原创 利用Internet Download Manager（IDM）软件批量下载网站上已选择的数据

以GLDAS数据为例，展示利用IDM软件进行批量下载