1. 项目概述InSAR数据处理与绘图的“瑞士军刀”集如果你正在处理合成孔径雷达干涉测量InSAR数据无论是做形变监测、沉降分析还是地质灾害评估那你一定对从原始数据到最终出版级图件这个漫长流程中的“工具切换”深有体会。SARscape、GMTSAR、ISCE这些专业软件包固然强大但它们往往在数据格式转换、批量处理、以及最终成果图的精细美化环节留下空白。这时一套得心应手的命令行工具和脚本就成了连接各个孤岛、提升效率的关键。这个项目或者说这份经验总结就是关于如何将GMTGeneric Mapping Tools、bash/csh脚本以及Matlab组合起来构建一个高效、可复现的InSAR数据处理与绘图流水线。简单来说它解决的核心痛点是自动化与灵活性。GMT负责产出地理投影正确、美观专业的矢量地图底图bash或csh脚本取决于你的系统偏好像胶水一样把数据预处理、格式转换、调用GMT命令、调用Matlab进行数值分析等步骤串联起来而Matlab则擅长处理矩阵运算、相位解缠结果的后续分析、或者生成一些GMT不那么擅长的复杂统计图表。这套组合拳打下来你就能从一个只会点击图形界面的操作员进化成能精准控制每一个处理步骤、一键生成从原始干涉图到最终形变序列图所有中间产品的“流水线工程师”。无论是处理单对影像还是时序InSAR如SBAS、PSI的大量数据栈这套方法都能显著提升你的工作效率和结果的可重复性。2. 核心工具链选型与协同逻辑为什么是GMT Shell Matlab这个组合这背后是基于每个工具的核心优势和在InSAR流程中的自然分工。2.1 GMT地图绘制的“定海神针”GMT并非一个GIS软件而是一个由近百个命令行工具组成的集合专门用于处理地理数据并生成高质量的PostScript或现代格式如PDF, PNG图件。在InSAR绘图中它的不可替代性体现在投影精确InSAR结果本质是地理编码后的栅格数据如geotiff。GMT原生支持UTM、地理坐标等多种投影能确保你的形变图与行政边界、地形底图完美套合这是很多科学绘图软件或Matlab自带绘图函数难以媲美的。出版级质量GMT生成的矢量图如PS、PDF线条平滑、字体清晰完全满足学术期刊对图件分辨率的要求。你可以精细控制颜色条cpt、图例、比例尺、指北针等所有地图元素。批处理友好所有绘图参数都通过命令行选项指定极易嵌入脚本实现自动化成图。2.2 Bash/Csh流程自动化的“中枢神经”Shell脚本是整个流水线的大脑负责调度和文件管理。Bash在Linux/macOS或Windows的Git Bash、WSL中几乎是标配语法功能强大社区资源丰富是大多数人的首选。Csh/Tcsh在某些历史较久的地球物理或超算环境中仍有使用其语法如set变量、foreach循环对某些用户来说更直观。选择哪一个取决于你的工作环境和团队习惯。本项目会兼顾两者给出关键语法对照。核心任务脚本负责遍历数据文件夹、批量转换格式例如用gdal_translate将h5或二进制文件转为GMT可读的netCDF或grd、构建并执行复杂的GMT绘图命令链、调用Matlab处理中间数据并管理临时文件。2.3 Matlab数值分析与特色图件的“专业顾问”Matlab在这个链条中扮演两个角色数据处理器对于相位解缠后需要进行的时空滤波、相位到形变的转换、时间序列分析、模型拟合如线性速率、季节性信号提取等涉及矩阵运算和复杂算法的步骤Matlab脚本比Shell更合适。补充绘图器当需要绘制非地图类图表如某个点的时间序列图、形变剖线图、统计直方图、相关性散点图时用Matlab可以快速实现并保持与数据处理部分一致的代码环境。2.4 协同工作流示例一个典型的时序InSAR成果图生成流程可能是Bash脚本启动遍历所有解缠后的相位文件.unw。格式转换在Bash中调用GDAL将.unw文件通常是二进制头文件转换为GMT的grd格式。GMT绘图Bash脚本调用gmt grdimage绘制每一幅干涉图的相位图并统一添加比例尺、颜色条。Matlab介入Bash脚本将一系列grd文件路径传递给Matlab脚本。Matlab读取这些栅格进行时序反演如SBAS计算平均形变速率和时序位移。结果回传Matlab将计算得到的速率图grd格式和指定点的时序数据文本格式输出。GMT最终成图Bash脚本再次调用GMT用速率图grd绘制主图并用gmt plot将Matlab生成的时序数据文本绘制成小图插入主图角落。 整个流程通过一个主控脚本可能是Bash来调度实现从原始数据到包含形变速率和典型点时间序列的复合出版图件的全自动生成。3. 关键命令与语法实例详解下面我们进入实战环节拆解每个工具在InSAR流程中的关键命令和脚本写法。3.1 GMT常用命令模块GMT命令繁多但用于InSAR绘图的核心模块集中在以下几个gmt grdconvert/gdal_translate数据输入桥梁。虽然GMT有grdconvert但处理复杂的SAR数据格式如ISCE输出的Erdas .img格式、ROI_PAC的.rsc.unw格式时GDAL库的gdal_translate命令往往更可靠。例如将ISCE生成的.geo.unw.geotiff转换为GMT的grd# Bash示例 gdal_translate -of NetCDF input.geo.unw.geo.tif phase.grd注意确保GMT编译时支持NetCDF并且GDAL版本与数据格式兼容。转换后务必用gmt grdinfo phase.grd检查网格范围、像素尺寸和单位是否正确。gmt grdimage绘制干涉相位或形变栅格图的核心命令。关键参数包括gmt grdimage phase.grd -R113.5/114.5/22.0/23.0 -JM15c -Cphase.cpt -Baf -BWSen -Ia15nt0.5 -P output.ps-R指定区域经度/纬度。务必与你的数据区域严格一致可以从.grd文件的头信息中获取。-JM设置墨卡托投影和地图宽度。-C指定颜色表.cpt文件。InSAR相位图常用循环色系如polar形变图常用线性色系如vik。可以使用gmt makecpt自定义。-I添加光照效果山体阴影-a15是方位角nt0.5是透明度能让地形起伏感更强突出干涉条纹。-B绘制地图边框、刻度及注释。-Baf是自动添加主要和次要刻度-BWSen表示在西边和南边绘制边框和刻度。gmt pscoast叠加海岸线、国界、河流等地理要素。这是让InSAR图具有地理参考意义的关键一步。gmt pscoast -R -J -Df -W0.5p,black -Glightgray -N1/0.5p,red -Lg115/22c22w50kfu -O -K output.ps-Df使用全分辨率海岸线数据需提前下载GMT的GSHHG数据。-W绘制海岸线。-G陆地填充色。-N1绘制国界线注意数据源的敏感性和绘图用途学术出版需谨慎。-L添加比例尺。这个参数非常实用能自动计算并标注比例尺长度和单位。gmt psscale添加颜色条。颜色条是科学图件的灵魂必须清晰准确。gmt psscale -Cphase.cpt -Dx15c/5cw12c/0.5ch -BxaflPhase (rad) -BylCycles -O output.ps-D精确定位颜色条。x15c/5c表示颜色条左下角位于页面坐标(15cm, 5cm)处w12c是宽度h表示水平放置默认为垂直v。-B设置颜色条刻度注释。l参数用于添加标签。3.2 Bash脚本编程要点Bash脚本用于串联上述GMT命令并处理文件。循环处理批量文件#!/bin/bash # 遍历当前目录下所有 .unw.grd 文件 for grd_file in *.unw.grd; do # 提取文件名不含后缀 base_name$(basename $grd_file .unw.grd) # 构建输出文件名 ps_file${base_name}.ps # 执行GMT绘图命令 gmt begin $base_name ps gmt grdimage $grd_file -R... -J... -C... gmt pscoast -R -J -Df -W... gmt psscale -C... -D... gmt end # 将PS转换为PDF gmt psconvert $ps_file -A -Tf echo 已完成: $base_name done实操心得在循环体内使用变量替换命令参数时务必用双引号包裹变量如$grd_file以防止文件名中含有空格时脚本报错。这是新手常踩的坑。参数化与配置文件对于固定的研究区可以将-R, -J等参数定义为变量甚至写入一个单独的配置文件config.sh用source config.sh引入提高脚本的可维护性。# config.sh REGION113.5/114.5/22.0/23.0 PROJECTIONM15c CPT_FILEmy_deformation.cpt# main.sh source config.sh gmt grdimage data.grd -R$REGION -J$PROJECTION -C$CPT_FILE ...3.3 Csh脚本语法对照Csh的语法与Bash差异较大主要注意变量设置和循环。变量设置与引用#!/bin/csh set region 113.5/114.5/22.0/23.0 set projection M15c gmt grdimage input.grd -R$region -J$projection ... # 注意变量赋值用 set引用时直接使用 $变量名循环处理foreach grd_file (*.unw.grd) set base_name basename $grd_file .unw.grd # 使用反引号执行命令并赋值 gmt begin $base_name ps # ... GMT命令 gmt end gmt psconvert $base_name.ps -A -Tf echo 已完成: $base_name end3.4 Matlab数据处理与衔接Matlab在这里主要做两件事读GMT的grd文件进行分析以及输出GMT可读的文本或grd文件。读取GMT grd文件GMT的grdNetCDF格式可以用Matlab的ncread函数直接读取。% 读取形变速率grd文件 filename velocity.grd; % 注意GMT grd文件可能使用‘x’, ‘y’, ‘z’作为变量名也可能用‘lon’, ‘lat’, ‘z’ try x ncread(filename, x); % 或 lon y ncread(filename, y); % 或 lat vel ncread(filename, z); catch % 如果变量名不对尝试读取所有变量信息 info ncinfo(filename); disp(info.Variables); end % 将x, y网格化为矩阵格式便于绘图和分析 [X, Y] meshgrid(x, y); % 注意vel矩阵的方向可能与Matlab的meshgrid预期不一致有时需要转置(transpose)或翻转(flipud)注意事项GMT和Matlab对矩阵的行列存储顺序行优先 vs 列优先和坐标系左上角原点 vs 左下角原点定义可能不同。这会导致读入的矩阵图像“倒置”或“镜像”。一个常见的解决方法是在Matlab中读取后使用vel flipud(vel);或类似操作进行校正。务必用imagesc(x, y, vel)初步显示与GMT原图对比验证。输出文本供GMT绘图将某个点的时序形变数据输出为GMT的plot命令可读的文本。% 假设有时间和形变数据 time [2018.0, 2018.5, 2019.0, ...]; % 十进制年 deformation [0, 5.2, -3.1, ...]; % 毫米 % 保存为两列文本 data_out [time(:), deformation(:)]; save(time_series.txt, data_out, -ascii); % 在Bash脚本中后续可以用 gmt plot time_series.txt -R... -B... -W2p,red 来绘制曲线输出grd文件将Matlab计算出的新栅格如滤波后的形变场写回为GMT可读的grd。% 假设有新的网格数据 new_vel以及对应的xvec, yvec向量 % 创建NetCDF文件 ncid netcdf.create(filtered_velocity.grd, CLOBBER); % 定义维度 dimid_x netcdf.defDim(ncid, x, length(xvec)); dimid_y netcdf.defDim(ncid, y, length(yvec)); % 定义变量 varid_x netcdf.defVar(ncid, x, double, dimid_x); varid_y netcdf.defVar(ncid, y, double, dimid_y); varid_z netcdf.defVar(ncid, z, double, [dimid_x, dimid_y]); netcdf.endDef(ncid); % 写入数据 netcdf.putVar(ncid, varid_x, xvec); netcdf.putVar(ncid, varid_y, yvec); netcdf.putVar(ncid, varid_z, new_vel); % 注意转置Matlab是列优先NetCDF通常是行优先。 netcdf.close(ncid);这个过程较为繁琐。更简单的方法是使用第三方工具箱如gmtmexGMT官方提供的Matlab接口或export_fig社区中的一些辅助函数。但掌握原生NetCDF写入有助于理解数据交换的本质。4. 完整实操案例从干涉图到形变速率剖面图我们通过一个完整案例将上述所有知识点串联起来。目标处理一个干涉对生成的形变栅格los_disp.grd单位米绘制带有地理背景的形变填色图并在图上画一条剖面线AB提取并绘制该剖面的形变曲线。4.1 步骤一准备环境与数据假设我们已在Bash环境下拥有以下文件los_disp.grd视线向形变栅格文件。config.sh配置文件定义了区域、投影等。profile_coords.txt文本文件包含剖面线起点A和终点B的经纬度每行一个点经度 纬度。113.6 22.2 114.2 22.84.2 步骤二主绘图Bash脚本plot_deformation.sh#!/bin/bash # 加载配置 source config.sh # 1. 启动GMT现代模式会话直接生成PDF gmt begin deformation_map pdf # 2. 绘制形变栅格图使用viridis色系 gmt grdimage los_disp.grd -R$REGION -J$PROJECTION -Cviridis -Ia15nt0.2 # 3. 叠加高分辨率海岸线 gmt coast -R -J -Df -W0.8p,black -G240/240/240 -N1/0.5p,50/50/50 # 4. 添加颜色条 gmt colorbar -Cviridis -Dx15c/-1cw12c/0.5ch -Bxa0.1f0.02lLOS Displacement (m) -Byl # 5. 绘制剖面线位置 gmt plot profile_coords.txt -R -J -W2p,red,solid -lProfile A-B # 6. 在起点和终点添加标记 gmt plot -R -J -Sc0.3c -Gred -W0.5p,black EOF 113.6 22.2 114.2 22.8 EOF # 7. 添加比例尺和指北针 gmt basemap -R -J -Lg113.7/22.1c22w20kfu -Tdg114.3/22.9w1cf2l gmt end echo 主形变图绘制完成deformation_map.pdf # 8. 提取剖面数据 # 使用gmt grdtrack沿剖面线采样 gmt grdtrack profile_coords.txt -Glos_disp.grd profile_data.txt # profile_data.txt 格式经度 纬度 距离(从起点算起,km) 形变量(m) # 9. 绘制剖面图 gmt begin profile pdf # 设置绘图区域X轴为距离(0-最大距离)Y轴为形变值(自动调整) # 先获取形变值的范围用于设置Y轴范围 min_max$(gmt info profile_data.txt -C -o5,6) # 获取第5列(距离)和6列(形变)的min/max # 拆分为变量 (假设info输出为xmin xmax ymin ymax) read xmin xmax ymin ymax $(echo $min_max) # 绘制剖面曲线距离单位转换为km形变单位转换为mm gmt plot profile_data.txt -i2,5 -R0/$xmax/$ymin/$ymax -JX15c/8c -W2p,blue -BxaflDistance along profile (km) -ByafglLOS Displacement (m) -BWSen # 可选填充曲线与零线之间的区域 gmt plot profile_data.txt -i2,5 -R -J -Glightblue50 -t50 gmt end echo 剖面图绘制完成profile.pdf实操心得gmt grdtrack是提取剖面数据的利器。-i选项在gmt plot中用于指定输入数据的列索引从0开始。在脚本中通过gmt info和命令替换$()动态获取数据范围来设置-R参数能使脚本适应不同的数据更加通用。4.3 步骤三使用Matlab进行剖面数据的平滑与拟合有时直接从栅格中提取的剖面数据噪声较大我们需要在Matlab中进行平滑或多项式拟合。% profile_analysis.m data load(profile_data.txt); distance_km data(:, 3); % 第3列是距离(km) disp_m data(:, 4); % 第4列是形变(m) % 1. 移动平均平滑 window_size 5; % 5个点的窗口 smoothed_disp movmean(disp_m, window_size); % 2. 线性拟合假设形变是距离的线性函数 p polyfit(distance_km, disp_m, 1); fit_disp polyval(p, distance_km); % 3. 将平滑和拟合后的数据保存为新文件供GMT绘制 output_data [distance_km, disp_m*1000, smoothed_disp*1000, fit_disp*1000]; % 转换为毫米 header Distance(km) Original(mm) Smoothed(mm) LinearFit(mm); fid fopen(profile_processed.txt, w); fprintf(fid, %s\n, header); fclose(fid); dlmwrite(profile_processed.txt, output_data, -append, delimiter, \t, precision, %.4f); % 4. 也可以在Matlab中直接绘图对比 figure; plot(distance_km, disp_m*1000, k., MarkerSize, 8); hold on; plot(distance_km, smoothed_disp*1000, b-, LineWidth, 2); plot(distance_km, fit_disp*1000, r--, LineWidth, 2); xlabel(Distance along profile (km)); ylabel(LOS Displacement (mm)); legend(Original, [Smoothed (win, num2str(window_size), )], Linear Fit); grid on;然后可以在Bash脚本中调用Matlab处理数据再使用GMT绘制更精美的剖面图。# 在plot_deformation.sh末尾添加 matlab -batch profile_analysis -nosplash -nodesktop # 使用GMT绘制处理后的剖面 gmt begin enhanced_profile pdf gmt plot profile_processed.txt -i0,1 -R... -J... -W1p,gray -lOriginal gmt plot profile_processed.txt -i0,2 -R... -J... -W2p,blue -lSmoothed gmt plot profile_processed.txt -i0,3 -R... -J... -W2p,red,- -lLinear Fit gmt legend -DjTRo0.2c -Fgwhitep0.5p gmt end5. 常见问题、调试技巧与避坑指南在实际操作中你会遇到各种报错和意外情况。这里记录了一些典型问题及其解决方法。5.1 GMT相关报错与解决错误grdimage: Warning: 1 (of 1) grid file [xxx.grd] doesnt have a recognized grid format原因GMT无法识别网格文件格式。最常见原因是grd文件不是真正的NetCDF格式或者内部维度、变量名不符合GMT预期。排查用ncdump -h xxx.grd查看文件头信息。检查是否存在x,y,z或lon,lat,z变量。用gdalinfo xxx.grd查看是否能被GDAL识别。如果不能说明文件可能已损坏或格式特殊。解决使用gdal_translate进行格式转换是更稳妥的入口。确保输出格式为-of NetCDF。错误绘图区域-R与数据区域不匹配导致空白图或部分显示。原因-R参数设置错误或者数据本身的坐标范围用gmt grdinfo查看与预期不符。解决在脚本中使用命令替换自动获取数据范围region$(gmt grdinfo input.grd -I-) gmt grdimage input.grd -R$region -J...-I-选项会输出-R所需的min/max格式字符串。问题生成的PS/PDF文件颜色条或图例位置不理想。解决-D参数用于精确定位。理解其语法-D[g|j|J|n|x]refpointwwidth[/height][jjustify][odx[/dy]]是关键。g使用地图坐标定位需在-R -J之后。j/J使用相对定位如JMR表示地图内右下角。x使用页面坐标单位cm/inch。多调试几次找到最适合你图件布局的位置。可以先画一个简单的图确定参考点坐标。5.2 Shell脚本调试技巧脚本执行权限bash: ./script.sh: Permission denied解决chmod x script.sh变量未定义或命令未找到在脚本开头添加set -euxo pipefail。-e有错误立即退出-u使用未定义变量时报错-x打印执行的命令便于追踪-o pipefail管道中任何命令失败则整个管道失败。对于命令未找到检查命令是否在PATH中或使用绝对路径。路径中包含空格这是Shell脚本的经典陷阱。始终用双引号包裹变量。# 错误 gmt grdimage $input_file ... # 正确 gmt grdimage $input_file ...5.3 Matlab与GMT数据交换的“方向”陷阱这是最隐蔽的问题之一。Matlab的meshgrid生成的X, Y矩阵与GMT保存的grd数据在内存中的排列方式可能正好转置或翻转。症状在Matlab中imagesc(lon, lat, data)显示的图像与GMT用grdimage绘制的图像上下或左右颠倒。诊断与解决在Matlab中读取grd后同时显示size(data)和[length(lon), length(lat)]看维度是否匹配应该是[length(lat), length(lon)]。尝试不同的组合data,flipud(data),fliplr(data),flipud(data)。并与GMT原图对比。最可靠的方法在Matlab中用ncdisp(file.grd)查看变量详情注意是否有direction或order相关的属性。有时数据是按“行优先”存储的而Matlab是“列优先”。建立一个已知的小型测试网格例如5x5分别用GMT和Matlab生成、读取、显示来摸清转换规律。5.4 性能优化建议对于大批量绘图避免在循环中反复启动和关闭GMT会话。使用GMT现代模式的gmt begin和gmt end将一系列绘图命令包裹起来效率更高。减少文件I/O如果Matlab只是进行简单的矩阵运算考虑使用GMT自带的gmt grdmath进行网格计算如加减乘除、滤波这比在Matlab和GMT之间来回读写文件要快得多。并行处理如果处理成百上千个干涉图可以利用Shell的并行工具如GNU parallel或xargs -P来并行运行多个GMT绘图进程充分利用多核CPU。这套工具链的学习曲线初期可能有些陡峭尤其是需要同时熟悉GMT的数百个命令选项、Shell脚本的语法以及Matlab与外部数据的交互。但一旦掌握你将获得无与伦比的灵活性和自动化能力能够应对各种复杂的InSAR数据处理与可视化需求从重复劳动中解放出来更专注于科学问题本身。我的经验是从一个具体的小目标开始比如“自动画出我这幅干涉图”边做边学积累自己的代码片段库逐渐就能搭建起强大的个人分析流水线。
