InSAR数据处理与绘图自动化：GMT+Shell+Matlab工具链实战指南

1. 项目概述：InSAR数据处理与绘图的“瑞士军刀”集

如果你正在处理合成孔径雷达干涉测量（InSAR）数据，无论是做形变监测、沉降分析还是地质灾害评估，那你一定对从原始数据到最终出版级图件这个漫长流程中的“工具切换”深有体会。SARscape、GMTSAR、ISCE这些专业软件包固然强大，但它们往往在数据格式转换、批量处理、以及最终成果图的精细美化环节留下空白。这时，一套得心应手的命令行工具和脚本就成了连接各个孤岛、提升效率的关键。这个项目，或者说这份经验总结，就是关于如何将GMT（Generic 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在这个链条中扮演两个角色：

1. 数据处理器：对于相位解缠后需要进行的时空滤波、相位到形变的转换、时间序列分析、模型拟合（如线性速率、季节性信号提取）等涉及矩阵运算和复杂算法的步骤，Matlab脚本比Shell更合适。
2. 补充绘图器：当需要绘制非地图类图表，如某个点的时间序列图、形变剖线图、统计直方图、相关性散点图时，用Matlab可以快速实现并保持与数据处理部分一致的代码环境。

2.4 协同工作流示例一个典型的时序InSAR成果图生成流程可能是：

1. Bash脚本启动：遍历所有解缠后的相位文件（.unw）。
2. 格式转换：在Bash中调用GDAL，将.unw文件（通常是二进制+头文件）转换为GMT的grd格式。
3. GMT绘图：Bash脚本调用gmt grdimage绘制每一幅干涉图的相位图，并统一添加比例尺、颜色条。
4. Matlab介入：Bash脚本将一系列grd文件路径传递给Matlab脚本。Matlab读取这些栅格，进行时序反演（如SBAS），计算平均形变速率和时序位移。
5. 结果回传：Matlab将计算得到的速率图（grd格式）和指定点的时序数据（文本格式）输出。
6. 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 -I+a15+nt0.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/22+c22+w50k+f+u -O -K >> output.ps

  • -Df：使用全分辨率海岸线数据（需提前下载GMT的GSHHG数据）。
  • -W：绘制海岸线。
  • -G：陆地填充色。
  • -N1：绘制国界线（注意数据源的敏感性和绘图用途，学术出版需谨慎）。
  • -L：添加比例尺。这个参数非常实用，能自动计算并标注比例尺长度和单位。

• gmt psscale：添加颜色条。颜色条是科学图件的灵魂，必须清晰准确。

  gmt psscale -Cphase.cpt -Dx15c/5c+w12c/0.5c+h -Bxaf+l"Phase (rad)" -By+l"Cycles" -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 REGION="113.5/114.5/22.0/23.0" PROJECTION="M15c" CPT_FILE="my_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" end

3.4 Matlab数据处理与衔接Matlab在这里主要做两件事：读GMT的grd文件进行分析，以及输出GMT可读的文本或grd文件。

• 读取GMT grd文件：GMT的grd（NetCDF格式）可以用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);

  这个过程较为繁琐。更简单的方法是使用第三方工具箱，如gmtmex（GMT官方提供的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.8

4.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 -I+a15+nt0.2 # 3. 叠加高分辨率海岸线 gmt coast -R -J -Df -W0.8p,black -G240/240/240 -N1/0.5p,50/50/50 # 4. 添加颜色条 gmt colorbar -Cviridis -Dx15c/-1c+w12c/0.5c+h -Bxa0.1f0.02+l"LOS Displacement (m)" -By+l"" # 5. 绘制剖面线位置 gmt plot profile_coords.txt -R -J -W2p,red,solid -l"Profile 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.1+c22+w20k+f+u -Tdg114.3/22.9+w1c+f2+l 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 -Bxaf+l"Distance along profile (km)" -Byafg+l"LOS Displacement (m)" -BWSen # 可选：填充曲线与零线之间的区域 gmt plot profile_data.txt -i2,5 -R -J -Glightblue@50 -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 -l"Original" gmt plot profile_processed.txt -i0,2 -R... -J... -W2p,blue -l"Smoothed" gmt plot profile_processed.txt -i0,3 -R... -J... -W2p,red,- -l"Linear Fit" gmt legend -DjTR+o0.2c -F+gwhite+p0.5p gmt end

5. 常见问题、调试技巧与避坑指南

在实际操作中，你会遇到各种报错和意外情况。这里记录了一些典型问题及其解决方法。

5.1 GMT相关报错与解决

• 错误：grdimage: Warning: 1 (of 1) grid file [xxx.grd] doesn't have a recognized grid format

  • 原因：GMT无法识别网格文件格式。最常见原因是grd文件不是真正的NetCDF格式，或者内部维度、变量名不符合GMT预期。
  • 排查：
    1. 用ncdump -h xxx.grd查看文件头信息。检查是否存在x,y,z（或lon,lat,z）变量。
    2. 用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]refpoint+wwidth[/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绘制的图像，上下或左右颠倒。
• 诊断与解决：
  1. 在Matlab中读取grd后，同时显示size(data)和[length(lon), length(lat)]，看维度是否匹配（应该是[length(lat), length(lon)]？）。
  2. 尝试不同的组合：data',flipud(data),fliplr(data),flipud(data')。并与GMT原图对比。
  3. 最可靠的方法：在Matlab中，用ncdisp('file.grd')查看变量详情，注意是否有direction或order相关的属性。有时，数据是按“行优先”存储的，而Matlab是“列优先”。
  4. 建立一个已知的小型测试网格（例如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数据处理与可视化需求，从重复劳动中解放出来，更专注于科学问题本身。我的经验是，从一个具体的小目标开始（比如“自动画出我这幅干涉图”），边做边学，积累自己的代码片段库，逐渐就能搭建起强大的个人分析流水线。