小标题和正文1. m_map凭什么成了InSAR绘图的默认选项做InSAR时间序列的人应该都有同感数据解算完只是第一步真正让人头疼的是最后那几张图。形变速率图、累积形变图、时序曲线图每一张都要往论文或报告里放既要信息完整又要符合SCI期刊的审美底线。我最早用的是MATLAB自带的绘图函数结果画出来的图自带一股“学术丑图”的味道——海岸线粗糙坐标轴比例失衡研究区域周围的地理信息少得可怜放在论文里审稿人一眼就能看出来作者在可视化上没花心思。后来接触到m_map算是彻底把这个问题解决了。m_map是MATLAB环境下的一套开源地图绘制工具箱最早由海洋科学领域的研究人员开发维护初衷是给海洋、大气、地质学科提供一套不需要依赖商用GIS的地图底图解决方案。它内置了全球海岸线数据、国界数据、地形数据支持二十多种地图投影并且可以和MATLAB原生的绘图函数无缝混用。对于InSAR形变结果来说m_map能做的其实非常聚焦把LOS向或垂直向的形变速率值按照经纬度坐标直接铺在地图底图上同时控制好比例尺、方向指示、图例和裁剪区域。在InSAR圈子里m_map几乎成了默认选项。原因很简单InSAR数据的最终表达本质上是在地理坐标系下的栅格可视化而m_map恰好把“地图底图”“投影转换”“空间叠加”这三件事打包好了。相比GMTm_map对MATLAB用户友好得多不需要重新学一整套脚本语法相比ArcGISm_map又轻量得多不用频繁来回导入导出直接在MATLAB工作区里就能把形变场和地理要素画在一起。这一篇我就从零开始把m_map绘制InSAR形变速率图的完整流程拆开讲一遍。内容覆盖数据准备、绘图步骤、配色方案、常见报错和坑点基本照着走一遍就能画出能直接放进论文的图。适合作InSAR数据处理的研究生或者刚转入射干方向的科研人员参考。2. 画图前先把数据收拾利索从时序解算结果到绘图数组很多人拿到m_map就急着画图结果第一步就翻车。原因在于m_map本身只是一个绘图工具它不会帮你处理数据。你喂给它的必须是一组干净的、已经完成地理编码的经纬度网格和对应的形变速率值这一点和GMT绘图的逻辑是一样的。2.1 数据存储形式经纬度网格和速率矩阵是怎么来的以常规的SBAS-InSAR或PS-InSAR处理流程为例做完时序解算之后得到的形变速率通常是以GeoTIFF或者二进制栅格文件存储的。打开之后里面有四个核心数组经度数组单位是度纬度数组单位是度形变速率值单位通常是毫米/年质量图或相干性图用于掩膜。这里最容易犯的错误是直接读取GeoTIFF的像素矩阵而不关心地理参考信息。比如你在某软件里导出速率图时如果导出范围是一个矩形那么对应的经度和纬度实际上是按矩阵行列规则排列的网格而不是每个像素都单独带坐标。换句话说你需要先通过影像的仿射变换参数一般是六参数模型左上角坐标、像元尺寸、旋转项把行列号换算成经纬度。MATLAB里读取GeoTIFF通常用geotiffread或者readgeoraster后者是较新版本推荐的方式。以readgeoraster为例读取后会返回一个地理参考对象RR里有CoordinateGeographicCRS属性以及用于行列到经纬度换算的仿射变换矩阵。换算公式很简单lon R.RasterSize(2)T2? 不对列方向才是经度实际上更稳妥的做法是先看R的RasterSize再按下面的方式生成网格[rows, cols] size(velocity); lon0 R.XWorldLimits(1); % 左边界经度 lat0 R.YWorldLimits(2); % 上边界纬度 cellSize R.CellExtentInWorldX; % 像元尺寸 lon lon0 (0:cols-1) * cellSize; lat lat0 - (0:rows-1) * cellSize; [LON, LAT] meshgrid(lon, lat);这里唯一的坑是lat必须递减生成因为GeoTIFF的行方向是从北往南排的这一点跟一般的笛卡尔坐标系直觉相反。2.2 掩膜操作和不值剔除不处理干净画出来的图没法看InSAR解算结果里低相干区、失相干区、水体、植被覆盖重的地方速率值往往是NaN或者异常值。如果不做掩膜直接画会出现两种灾难性结果一种是图上大片空白或零值覆盖导致色标被压缩原本应该明显的形变信号在色带上只占一小段区间图看起来平平淡淡没有层次感另一种是残留的离群点把色标拉宽比如某个点速率标了-200毫米/年而研究区域实际形变只有±30毫米/年结果全部颜色都被这一个坏点带偏了真正的形变区反而显示不出差异。我常用的处理方式是分三步第一步把绝对值超过阈值的数据直接置NaN。比如库区形变一般在±50毫米/年以内超过±100毫米/年的基本可以判定为大气相位残余或解缠错误直接剔除。velocity(abs(velocity) 100) NaN;第二步用相干性或质量图做掩膜。相干性低于0.3的像素点相位解缠可信度不高对应速率不可信应当置NaN。coherence(coherence 0.3) NaN; velocity(isnan(coherence)) NaN;第三步检查是否存在经纬度坐标为NaN的边缘行或列。很多处理软件在裁剪过程中会在边缘留一圈NaN如果带着这一圈NaN去画图填色的时候会出现明显的边缘断层图面效果非常脏。处理完掩膜之后我一般会顺手做一个统计输出看看NaN占比是多少、速率最大值最小值、均值及标准差。这个输出既是质量控制的记录也能让你在论文方法部分写清楚用了什么筛选条件。fprintf(有效像元比例%.2f%%\n, (1 - sum(isnan(velocity(:)))/numel(velocity)) * 100); fprintf(速率范围%.2f %.2f mm/yr\n, min(velocity(:), [], omitnan), max(velocity(:), [], omitnan));不要小看这组数字很多时候你画出来的图色标很奇怪问题根源就是在数据处理阶段没有做好统计和质量筛选而不是绘图阶段出了问题。2.3 数据坐标参考系的检查WGS84和投影坐标别混着用再强调一个低级但高频的错误坐标参考系不统一。m_map默认使用的是经纬度坐标也就是地理坐标系基准面通常是WGS84。但有些InSAR处理流程输出的是投影坐标系比如UTM坐标坐标单位是米。直接把UTM坐标当成经纬度喂给m_map画出来的图坐标系就会完全错掉海岸线和形变位置对不上偏差通常大到几百公里。如果你拿到的是UTM坐标的栅格有两个选择用MATLAB的projinv函数从投影坐标反算经纬度回到InSAR处理软件里重新以WGS84坐标导出GeoTIFF。% 已知UTM坐标东向、北向和投影带号 [lat, lon] projinv(projcrs(32650), x, y); % 32650是UTM 50N的EPSG代码在我做过的几个项目里第二种方式更省时间。因为反算出来的经纬度网格往往不是规则的矩形网格后续在m_map里虽然也可以画但遇到数据范围较大的情况时容易出现边缘变形。直接在源头上用地理坐标导出网格整齐绘图时省心很多。3. m_map核心绘制流程从海岸线到速率填图的五步走数据准备好之后真正调用m_map的绘图流程反而不复杂。整个流程可以拆成五个核心步骤每一行代码都有明确目的。下面按我日常写脚本的顺序完整过一遍。3.1 初始化地图投影和绘图范围任何一幅m_map图都从m_proj开始。这一行代码决定了你整张图的地图投影类型和显示范围。InSAR形变图最常用的是lambert等角圆锥投影lambert和横轴墨卡托投影mercator。m_proj(lambert, long, [lon_min lon_max], lat, [lat_min lat_max]);这里lon和lat的范围不需要跟数据范围完全一致通常比数据范围略大一些周围留一点边距否则图面会显得很拥挤海岸线也贴到图框边上了。我一般留2%左右的外边距lon_min min(LON(:)) - 0.02 * (max(LON(:)) - min(LON(:))); lon_max max(LON(:)) 0.02 * (max(LON(:)) - min(LON(:))); lat_min min(LAT(:)) - 0.02 * (max(LAT(:)) - min(LAT(:))); lat_max max(LAT(:)) 0.02 * (max(LAT(:)) - min(LAT(:)));投影类型的选择逻辑是这样的如果研究区域纬度跨度大且形状呈南北走向lambert投影变形小如果区域面积较小且更关注局部细节mercator也能用但高纬地区面积被放大的问题要注意。还有一种情况是极地研究那就用polar stereographic不过InSAR常规应用里很少遇到。3.2 加载海岸线和地理要素地图底图是m_map的核心价值所在。m_coast和m_gshhs这两个函数提供了从粗糙到精细的全球海岸线数据。m_coast(patch, [0.85 0.85 0.85], edgecolor, none); m_gshhs(patch, [0.85 0.85 0.85], edgecolor, none);两者区别在于m_coast自带的是相对低分辨率的数据全球范围内使用没有问题但如果你放大了看某个河口或半岛区域岸线形状可能跟实际偏差几百米m_gshhs可以加载更精细的分辨率等级最细的是full但数据量大绘制速度慢。我的习惯是研究区域范围小于1°×1°时用m_gshhs(fine)或full保证岸线和真实海岸尽可能贴合范围大于几个度时m_coast就足够了精细数据反而会因为岸线太密把视线干扰掉。关于岸线的陆地填充颜色我通常用[0.85 0.85 0.85]这种浅灰色。原因是形变速率图的色标通常是从蓝色到红色如果陆地填充太重会和冷色系形变区产生视觉冲突浅灰色既能区分陆地和海洋又不会抢占形变信号本身的视觉焦点。3.3 绘制形变速率填色图并控制颜色映射核心绘图命令是m_pcolor它是m_map对MATLAB原生pcolor的封装。用法上跟pcolor几乎一样但有两点必须注意第一pcolor默认会把单元格索引对应的四个角点连线并填色如果你直接传LON、LAT、velocity得到的结果会有一个格子偏移图形看起来像是“错位”的。解决办法是使用m_pcolor直接传网格坐标即可m_map内部会做处理或者用imagesc配合m_plot来打点。第二MATLAB的pcolor存在“栅格线”问题即图上会显示单元格之间的细线。解决这个问题最常用的办法还是用m_pcolor它内部已经把EdgeColor设置成了none。如果你用的是原生pcolor记得手动设set(h, EdgeColor, none);我的标准绘图代码是h m_pcolor(LON, LAT, velocity); set(h, EdgeColor, none); colormap(jet); % 或者自定义色带 caxis([-vmax vmax]);其中caxis的设定非常关键。InSAR形变速率图通常采用对称色标即正负速率使用对称的颜色区间这样形变符号对比一目了然。vmax一般选择速率标准差的2~3倍或者取绝对值最大值向上取整。3.4 添加色标、比例尺、指北针和地形阴影形变图不只是填色还需要一系列地图要素让读者读懂空间信息。m_map提供了一套专门的函数m_rbcolor % 添加垂直色标 m_scale % 添加比例尺不过实际用下来我对m_rbcolor的排版效果一直不太满意它的位置控制相对固定很难精确调整到和图文对齐。所以我一般采用colorbar原生函数来代替c colorbar(southoutside); c.Label.String LOS Velocity (mm/yr); c.TickLabels strsplit(num2str(ticks));这样调整起来更灵活尤其是需要统一多幅子图色标的时候可以用set(c, Position, ...) 精确控制位置。比例尺用m_scale指定一个位置坐标它会在图上绘制一段带刻度的线段m_scale([lon_corner lat_corner], 10, k); % 10表示比例尺代表的实际距离公里指北针我建议画在图框外的角落避免遮挡数据区域。m_map没有内置指北针函数但可以通过m_plot加annotation的方式实现网上有成熟的northarrow脚本可以参考。3.5 叠加研究区边界、断层线和GPS测站这是InSAR论文里非常常见的组合图底图是形变速率上面叠加断层、边界或者离散点位。m_map的好处是同一套坐标体系里继续用m_plot、m_line、m_patch绘制地理要素坐标仍然传经纬度即可。m_plot(fault_lon, fault_lat, r-, linewidth, 1.5); m_plot(gps_lon, gps_lat, kv, markersize, 6, markerfacecolor, k);这一步看似简单但有一个隐藏要求断层线和边界的顺序要在填色图绘制完之后再叠加否则会被覆盖。填色图在最底层然后叠加岸线最后才是断层、点位、标注之类的矢量要素。4. 把图调好看且不误导读者色彩设计、NaN显示与裁剪细节形变速率图的颜色处理是最容易被忽视也最值得花时间的部分。你花的精力不会直接写进论文的文字里但审稿人和读者对图的第一印象往往决定了他们对数据质量的预判。4.1 为什么不要直接用默认的jet以及怎么自定义连续色带MATLAB默认的jet色带是从蓝到红视觉冲击力很强但它有个非常严重的感知缺陷颜色过渡不均匀。在jet里青色和黄色区域的亮度变化非常急促而蓝色和红色区域几乎不变这会导致读者对形变幅值的感知出现偏差。比如±5毫米/年的微小形变在jet色带上几乎难以区分因为它们在色标上对应的颜色跳变太微弱。更专业的选择是MATLAB内置的parula或者各类科学色带比如cmocean系列中的balance色带。balance是蓝白红的对称色带非常适合InSAR正负形变的表达而且是白底零值居中视觉上非常干净。cmocean并不是MATLAB官方内置的需要单独下载但我强烈建议纳入你的工具库。colormap(cmocean(balance));如果没有安装cmocean也可以用MATLAB推荐的线性色带构建方法cmap [0.00 0.00 0.55; 0.10 0.15 0.70; 0.25 0.45 0.90; 0.60 0.80 1.00; 1.00 1.00 1.00; 1.00 0.85 0.60; 1.00 0.60 0.30; 0.90 0.30 0.10; 0.55 0.00 0.00]; colormap(cmap);这种手写色带的逻辑是从深蓝到浅蓝到白色再到浅红到深红中间白表示零形变两侧对称。颜色数量和每段的RGB值可以按需增加但核心原则是——色带两端要深、中间要浅、深浅过渡要均匀避免某一段颜色跳变特别陡峭。4.2 NaN值的显示策略透明、补零还是插值很多人在caxis和NaN之间纠结。实际上m_pcolor配合pcolor时NaN值默认是不绘制的也就是像素点保持透明背后显示的是白纸或岸线的底图颜色。这个行为大多数情况下是理想的但如果你的研究区域内失相干区域占比很大图面上就会出现大片的空白反而分散读者注意力。处理策略要看具体场景如果是小范围的局部失相干河流、植被带保持透明即可读者能自然理解那是相干性差的区域如果失相干面积太大导致形变场看起来“稀疏”一般建议在图中用半透明的灰色多边形遮盖这些区域而不是把速率补成0。补0是最大的错误因为0在色标上对应的颜色是有物理意义“零形变”的会把失相干区域和真实零形变区域混淆。如果想在图上明确标记失相干区可以先把NaN区域提取出来再叠加一层灰色mask isnan(velocity); hold on; m_pcolor(LON, LAT, double(mask)); set(gca, Color, [0.7 0.7 0.7]); % 或者给mask矩阵设置特殊的colormap不过这个方案实施起来要注意colormap冲突同一张图只能有一套colormap灰色mask会和速率色带发生冲突。所以最终的折中方案是保持NaN透明通过论文文字说明“剔除相干性低于0.3的像素”这个方法简单直接审稿人也都接受。4.3 裁剪范围的经验法则让图面焦点集中但保留地理参考研究区域的裁剪范围会直接影响图的质量。范围太小周围缺乏地理参照物读者很难定位范围太大形变信号只占图面一小块有效信息密度低。我的经验法则是以形变区为中心外扩至包含至少两个以上有辨识度的地理要素。比如研究区附近有河流、湖泊、主要道路、大型城镇裁剪范围至少要把这些要素框进来。这样读者看图时可以通过这些熟悉的地理位置快速定位而不是盯着经纬度标注算坐标。具体实现时直接在m_proj里调整边界lonlim [lon_min - 0.3, lon_max 0.3]; latlim [lat_min - 0.3, lat_max 0.3]; m_proj(lambert, long, lonlim, lat, latlim);扩展多少取决于研究区尺度。一个经验参数是当研究区宽度约50公里时四周各外扩10~15公里宽度100公里以上时外扩20~30公里。外扩尺寸以经纬度表述时需要根据当地纬度做换算用deg2km函数辅助判断。5. 实测中躲不开的几个坑以及怎么绕过去这一部分写的是我这几年来用m_map画InSAR图时真实踩过的坑有一个算一个每个都让人花过不少冤枉时间。5.1 报错“index exceeds matrix dimensions”问题出在数据维度而不是代码遇到过好几次网上搜了一圈大部分答案都是让检查数组长度。真正的原因在InSAR场景里往往是这样GeoTIFF读取时velocity的维度是[rows, cols]但你生成的LON和LAT网格如果用了meshgrid输出的维度是[cols, rows]两者维度一致但顺序反了传进m_pcolor就直接报错或者画出旋转了90度的图。这个问题的排查方法是在绘图前用size检查所有变量的维度并强制统一[rows, cols] size(velocity); if size(LON, 1) cols size(LON, 2) rows LON LON; LAT LAT; end这行代码我几乎写进了每一个绘图脚本里。在没有严格约定时ROW变COL, COL变ROW是InSAR数据从外部软件导入MATLAB最常见的问题。5.2 图形显示正常但海岸线位置和形变偏移必须检查m_proj范围另一个让人困惑的问题是形变数据画出来位置正确但加上海岸线后海岸线和形变区明显错位几百米甚至几公里。这个问题的根源往往是数据坐标参考系与m_map默认坐标系不一致。举个我实际遇到的例子某次处理某区域数据处理软件导出的GeoTIFF虽然写的是WGS84但实际上是经过特定椭球拟合的局部坐标底图海岸线来自全球数据集直接用会有约500米的系统性偏移。解决方案是在GIS软件里检查数据是否正确或者用几个已知地面控制点道路交叉口、桥梁对比形变图与光学影像确认偏移量后手动平移数据。简单说在m_map里海岸线偏移不是绘图函数的问题而是数据坐标参考系没校准的问题。这一点在开始处理大批量数据前一定要确认否则后面对比验证会很难进行。5.3 色标上限、caxis、colormap三者关系混乱导致图面颜色失真画出来的图颜色和caxis设置对不上是新手最常见的问题。核心原因是MATLAB绘图的颜色映射逻辑每个像素点的颜色不是由它的数值直接决定的而是由“数值在caxis范围里的相对位置”与colormap的映射关系决定的。举一个典型错误colormap(jet); caxis([-30 30]); h m_pcolor(LON, LAT, velocity); % 先设置caxis再绘图如果先设置了caxis再绘图某些情况下图形对象不会立即更新导致画面颜色仍然按默认范围映射。正确顺序是先绘图再设置caxish m_pcolor(LON, LAT, velocity); set(h, EdgeColor, none); caxis([-vmax vmax]); colormap(cmap);另一个问题是同时设置多个子图时MATLAB会把所有子图共享同一套caxis。要避免这种情况每一张子图都要单独绘图并设置自己的caxis或者用caxis([-统一值 统一值])强制所有子图使用同一套颜色区间这在多图对比里反而是推荐做法因为色标不一致的对比图会误导读者。5.4 保存图片时文字变成方块或乱码MATLAB绘图完保存成PDF或EPS时中文标注乱码是老问题。InSAR图里的中文标注主要是研究区名称、图例文字虽然SCI期刊大多要求英文标注但很多中文论文或者阶段汇报仍然需要中文。我的做法是两点第一尽量用英文或拼音标注。这不是崇洋媚外而是从实际成图效果出发英文标注在任何格式下都不会有字体问题。第二如果一定要用中文保存使用高分辨率PNG不要用PDF/EPS矢量格式。因为MATLAB矢量输出时对中文字体的嵌入支持非常差换成PNG格式后文字是按像素输出的基本不会乱码。print(gcf, velocity_map, -dpng, -r600);600dpi的PNG用于论文完全够用于PPT汇报更是绰绰有余。5.5 多子图共用色标导致颜色错误叠加如果一张图里有多个子图比如LOS速率图、累积形变图、相干图一排放三张最容易出的问题是绘图顺序和colormap相互污染。正确做法是逐个绘图每一个子图进行完整的初始化、绘图、设置caxis和colormap流程。MATLAB的colormap是图形级别的属性一个Figure里有多个axes时colormap是全局共享的除非用colormap(ax, cmap)指定ax。subplot(1,3,1); m_pcolor(LON, LAT, velocity); caxis([-vmax vmax]); colormap(cmocean(balance)); colorbar; subplot(1,3,2); m_pcolor(LON, LAT, velocity_accum); caxis([-amax amax]); colormap(cmocean(balance)); colorbar;注意这里每次subplot后都要重新设置colormap不然第二张图会沿用上一张图的色带导致颜色错乱。6. 把m_map工具链融入日常InSAR出图流程的几点心得m_map并不是一个功能单一的小函数它本质上是一个画地图的框架。当你把它纳入自己的出图流程后效率提升是很明显的。我现在的标准流程是写完一套绘图脚本后后续每个研究区的出图只需要替换数据和边界坐标其他代码基本不动。脚本里包含了完整的投影初始化、岸线绘制、数据掩膜处理、填色绘图、要素叠加、色标调整和保存导出一个项目从头到尾的形变图质量都保持一致。另外m_map还支持将多个不同的数据集叠加到同一张底图上这就为InSAR结果和地质背景的综合分析提供了很大的便利。比如在研究某个滑坡体时我可以在同一张图上叠加速率场、断裂带位置、已有滑坡边界和GPS位移矢量信息层次分明审稿人看图就能直接抓住要点。如果研究区域比较大或者你需要更精细的底图数据也可以考虑和其他工具配合使用。m_map能加载用户自定义的高分辨率海岸线数据或DEM晕渲图配合上自己做的坡度、地形阴影出图效果会更接近专业制图软件的成品。有一点我特别想提醒m_map用的海岸线数据本身是全球范围内的中等精度数据在一些近岸精细化研究里岸线偏差可能达到数百米。如果论文中有严格的地表要素对比分析建议用更高精度的岸线或直接用光学影像叠加验证不要完全依赖m_map自带的海岸线。还有一个小技巧就是给m_map绘制的图设置统一的输出尺寸。很多人在MATLAB的figure窗口里直接保存结果拉出来的图比例时宽时窄放在论文里插图大小不一后期还得在排版软件里反复调整。我一般在绘图开始时用set(gcf, PaperPositionMode, manual); set(gcf, PaperUnits, centimeters); set(gcf, PaperPosition, [0 0 18 15]);这样可以保证每次输出的图尺寸一致后续排版省很多事。最后再说一点个人体会InSAR数据处理到最后图件往往是决定成果表达质量的临门一脚。数据解算得再好图件表达不清楚会大大降低成果的说服力。m_map虽然不是最现代的绘图方案但在MATLAB社区里成熟稳定、教程多、排错容易对InSAR研究者来说仍然是最稳妥的选择。把这套流程练熟你会发现出图这件事本身也可以变成一件让人有成就感的事。 
长河君 · 2026/10/13 20:24:19 · 阅读 1.2k m_map绘制InSAR形变速率图：从数据准备到出图全流程
小标题和正文1. m_map凭什么成了InSAR绘图的默认选项做InSAR时间序列的人应该都有同感数据解算完只是第一步真正让人头疼的是最后那几张图。形变速率图、累积形变图、时序曲线图每一张都要往论文或报告里放既要信息完整又要符合SCI期刊的审美底线。我最早用的是MATLAB自带的绘图函数结果画出来的图自带一股“学术丑图”的味道——海岸线粗糙坐标轴比例失衡研究区域周围的地理信息少得可怜放在论文里审稿人一眼就能看出来作者在可视化上没花心思。后来接触到m_map算是彻底把这个问题解决了。m_map是MATLAB环境下的一套开源地图绘制工具箱最早由海洋科学领域的研究人员开发维护初衷是给海洋、大气、地质学科提供一套不需要依赖商用GIS的地图底图解决方案。它内置了全球海岸线数据、国界数据、地形数据支持二十多种地图投影并且可以和MATLAB原生的绘图函数无缝混用。对于InSAR形变结果来说m_map能做的其实非常聚焦把LOS向或垂直向的形变速率值按照经纬度坐标直接铺在地图底图上同时控制好比例尺、方向指示、图例和裁剪区域。在InSAR圈子里m_map几乎成了默认选项。原因很简单InSAR数据的最终表达本质上是在地理坐标系下的栅格可视化而m_map恰好把“地图底图”“投影转换”“空间叠加”这三件事打包好了。相比GMTm_map对MATLAB用户友好得多不需要重新学一整套脚本语法相比ArcGISm_map又轻量得多不用频繁来回导入导出直接在MATLAB工作区里就能把形变场和地理要素画在一起。这一篇我就从零开始把m_map绘制InSAR形变速率图的完整流程拆开讲一遍。内容覆盖数据准备、绘图步骤、配色方案、常见报错和坑点基本照着走一遍就能画出能直接放进论文的图。适合作InSAR数据处理的研究生或者刚转入射干方向的科研人员参考。2. 画图前先把数据收拾利索从时序解算结果到绘图数组很多人拿到m_map就急着画图结果第一步就翻车。原因在于m_map本身只是一个绘图工具它不会帮你处理数据。你喂给它的必须是一组干净的、已经完成地理编码的经纬度网格和对应的形变速率值这一点和GMT绘图的逻辑是一样的。2.1 数据存储形式经纬度网格和速率矩阵是怎么来的以常规的SBAS-InSAR或PS-InSAR处理流程为例做完时序解算之后得到的形变速率通常是以GeoTIFF或者二进制栅格文件存储的。打开之后里面有四个核心数组经度数组单位是度纬度数组单位是度形变速率值单位通常是毫米/年质量图或相干性图用于掩膜。这里最容易犯的错误是直接读取GeoTIFF的像素矩阵而不关心地理参考信息。比如你在某软件里导出速率图时如果导出范围是一个矩形那么对应的经度和纬度实际上是按矩阵行列规则排列的网格而不是每个像素都单独带坐标。换句话说你需要先通过影像的仿射变换参数一般是六参数模型左上角坐标、像元尺寸、旋转项把行列号换算成经纬度。MATLAB里读取GeoTIFF通常用geotiffread或者readgeoraster后者是较新版本推荐的方式。以readgeoraster为例读取后会返回一个地理参考对象RR里有CoordinateGeographicCRS属性以及用于行列到经纬度换算的仿射变换矩阵。换算公式很简单lon R.RasterSize(2)T2? 不对列方向才是经度实际上更稳妥的做法是先看R的RasterSize再按下面的方式生成网格[rows, cols] size(velocity); lon0 R.XWorldLimits(1); % 左边界经度 lat0 R.YWorldLimits(2); % 上边界纬度 cellSize R.CellExtentInWorldX; % 像元尺寸 lon lon0 (0:cols-1) * cellSize; lat lat0 - (0:rows-1) * cellSize; [LON, LAT] meshgrid(lon, lat);这里唯一的坑是lat必须递减生成因为GeoTIFF的行方向是从北往南排的这一点跟一般的笛卡尔坐标系直觉相反。2.2 掩膜操作和不值剔除不处理干净画出来的图没法看InSAR解算结果里低相干区、失相干区、水体、植被覆盖重的地方速率值往往是NaN或者异常值。如果不做掩膜直接画会出现两种灾难性结果一种是图上大片空白或零值覆盖导致色标被压缩原本应该明显的形变信号在色带上只占一小段区间图看起来平平淡淡没有层次感另一种是残留的离群点把色标拉宽比如某个点速率标了-200毫米/年而研究区域实际形变只有±30毫米/年结果全部颜色都被这一个坏点带偏了真正的形变区反而显示不出差异。我常用的处理方式是分三步第一步把绝对值超过阈值的数据直接置NaN。比如库区形变一般在±50毫米/年以内超过±100毫米/年的基本可以判定为大气相位残余或解缠错误直接剔除。velocity(abs(velocity) 100) NaN;第二步用相干性或质量图做掩膜。相干性低于0.3的像素点相位解缠可信度不高对应速率不可信应当置NaN。coherence(coherence 0.3) NaN; velocity(isnan(coherence)) NaN;第三步检查是否存在经纬度坐标为NaN的边缘行或列。很多处理软件在裁剪过程中会在边缘留一圈NaN如果带着这一圈NaN去画图填色的时候会出现明显的边缘断层图面效果非常脏。处理完掩膜之后我一般会顺手做一个统计输出看看NaN占比是多少、速率最大值最小值、均值及标准差。这个输出既是质量控制的记录也能让你在论文方法部分写清楚用了什么筛选条件。fprintf(有效像元比例%.2f%%\n, (1 - sum(isnan(velocity(:)))/numel(velocity)) * 100); fprintf(速率范围%.2f %.2f mm/yr\n, min(velocity(:), [], omitnan), max(velocity(:), [], omitnan));不要小看这组数字很多时候你画出来的图色标很奇怪问题根源就是在数据处理阶段没有做好统计和质量筛选而不是绘图阶段出了问题。2.3 数据坐标参考系的检查WGS84和投影坐标别混着用再强调一个低级但高频的错误坐标参考系不统一。m_map默认使用的是经纬度坐标也就是地理坐标系基准面通常是WGS84。但有些InSAR处理流程输出的是投影坐标系比如UTM坐标坐标单位是米。直接把UTM坐标当成经纬度喂给m_map画出来的图坐标系就会完全错掉海岸线和形变位置对不上偏差通常大到几百公里。如果你拿到的是UTM坐标的栅格有两个选择用MATLAB的projinv函数从投影坐标反算经纬度回到InSAR处理软件里重新以WGS84坐标导出GeoTIFF。% 已知UTM坐标东向、北向和投影带号 [lat, lon] projinv(projcrs(32650), x, y); % 32650是UTM 50N的EPSG代码在我做过的几个项目里第二种方式更省时间。因为反算出来的经纬度网格往往不是规则的矩形网格后续在m_map里虽然也可以画但遇到数据范围较大的情况时容易出现边缘变形。直接在源头上用地理坐标导出网格整齐绘图时省心很多。3. m_map核心绘制流程从海岸线到速率填图的五步走数据准备好之后真正调用m_map的绘图流程反而不复杂。整个流程可以拆成五个核心步骤每一行代码都有明确目的。下面按我日常写脚本的顺序完整过一遍。3.1 初始化地图投影和绘图范围任何一幅m_map图都从m_proj开始。这一行代码决定了你整张图的地图投影类型和显示范围。InSAR形变图最常用的是lambert等角圆锥投影lambert和横轴墨卡托投影mercator。m_proj(lambert, long, [lon_min lon_max], lat, [lat_min lat_max]);这里lon和lat的范围不需要跟数据范围完全一致通常比数据范围略大一些周围留一点边距否则图面会显得很拥挤海岸线也贴到图框边上了。我一般留2%左右的外边距lon_min min(LON(:)) - 0.02 * (max(LON(:)) - min(LON(:))); lon_max max(LON(:)) 0.02 * (max(LON(:)) - min(LON(:))); lat_min min(LAT(:)) - 0.02 * (max(LAT(:)) - min(LAT(:))); lat_max max(LAT(:)) 0.02 * (max(LAT(:)) - min(LAT(:)));投影类型的选择逻辑是这样的如果研究区域纬度跨度大且形状呈南北走向lambert投影变形小如果区域面积较小且更关注局部细节mercator也能用但高纬地区面积被放大的问题要注意。还有一种情况是极地研究那就用polar stereographic不过InSAR常规应用里很少遇到。3.2 加载海岸线和地理要素地图底图是m_map的核心价值所在。m_coast和m_gshhs这两个函数提供了从粗糙到精细的全球海岸线数据。m_coast(patch, [0.85 0.85 0.85], edgecolor, none); m_gshhs(patch, [0.85 0.85 0.85], edgecolor, none);两者区别在于m_coast自带的是相对低分辨率的数据全球范围内使用没有问题但如果你放大了看某个河口或半岛区域岸线形状可能跟实际偏差几百米m_gshhs可以加载更精细的分辨率等级最细的是full但数据量大绘制速度慢。我的习惯是研究区域范围小于1°×1°时用m_gshhs(fine)或full保证岸线和真实海岸尽可能贴合范围大于几个度时m_coast就足够了精细数据反而会因为岸线太密把视线干扰掉。关于岸线的陆地填充颜色我通常用[0.85 0.85 0.85]这种浅灰色。原因是形变速率图的色标通常是从蓝色到红色如果陆地填充太重会和冷色系形变区产生视觉冲突浅灰色既能区分陆地和海洋又不会抢占形变信号本身的视觉焦点。3.3 绘制形变速率填色图并控制颜色映射核心绘图命令是m_pcolor它是m_map对MATLAB原生pcolor的封装。用法上跟pcolor几乎一样但有两点必须注意第一pcolor默认会把单元格索引对应的四个角点连线并填色如果你直接传LON、LAT、velocity得到的结果会有一个格子偏移图形看起来像是“错位”的。解决办法是使用m_pcolor直接传网格坐标即可m_map内部会做处理或者用imagesc配合m_plot来打点。第二MATLAB的pcolor存在“栅格线”问题即图上会显示单元格之间的细线。解决这个问题最常用的办法还是用m_pcolor它内部已经把EdgeColor设置成了none。如果你用的是原生pcolor记得手动设set(h, EdgeColor, none);我的标准绘图代码是h m_pcolor(LON, LAT, velocity); set(h, EdgeColor, none); colormap(jet); % 或者自定义色带 caxis([-vmax vmax]);其中caxis的设定非常关键。InSAR形变速率图通常采用对称色标即正负速率使用对称的颜色区间这样形变符号对比一目了然。vmax一般选择速率标准差的2~3倍或者取绝对值最大值向上取整。3.4 添加色标、比例尺、指北针和地形阴影形变图不只是填色还需要一系列地图要素让读者读懂空间信息。m_map提供了一套专门的函数m_rbcolor % 添加垂直色标 m_scale % 添加比例尺不过实际用下来我对m_rbcolor的排版效果一直不太满意它的位置控制相对固定很难精确调整到和图文对齐。所以我一般采用colorbar原生函数来代替c colorbar(southoutside); c.Label.String LOS Velocity (mm/yr); c.TickLabels strsplit(num2str(ticks));这样调整起来更灵活尤其是需要统一多幅子图色标的时候可以用set(c, Position, ...) 精确控制位置。比例尺用m_scale指定一个位置坐标它会在图上绘制一段带刻度的线段m_scale([lon_corner lat_corner], 10, k); % 10表示比例尺代表的实际距离公里指北针我建议画在图框外的角落避免遮挡数据区域。m_map没有内置指北针函数但可以通过m_plot加annotation的方式实现网上有成熟的northarrow脚本可以参考。3.5 叠加研究区边界、断层线和GPS测站这是InSAR论文里非常常见的组合图底图是形变速率上面叠加断层、边界或者离散点位。m_map的好处是同一套坐标体系里继续用m_plot、m_line、m_patch绘制地理要素坐标仍然传经纬度即可。m_plot(fault_lon, fault_lat, r-, linewidth, 1.5); m_plot(gps_lon, gps_lat, kv, markersize, 6, markerfacecolor, k);这一步看似简单但有一个隐藏要求断层线和边界的顺序要在填色图绘制完之后再叠加否则会被覆盖。填色图在最底层然后叠加岸线最后才是断层、点位、标注之类的矢量要素。4. 把图调好看且不误导读者色彩设计、NaN显示与裁剪细节形变速率图的颜色处理是最容易被忽视也最值得花时间的部分。你花的精力不会直接写进论文的文字里但审稿人和读者对图的第一印象往往决定了他们对数据质量的预判。4.1 为什么不要直接用默认的jet以及怎么自定义连续色带MATLAB默认的jet色带是从蓝到红视觉冲击力很强但它有个非常严重的感知缺陷颜色过渡不均匀。在jet里青色和黄色区域的亮度变化非常急促而蓝色和红色区域几乎不变这会导致读者对形变幅值的感知出现偏差。比如±5毫米/年的微小形变在jet色带上几乎难以区分因为它们在色标上对应的颜色跳变太微弱。更专业的选择是MATLAB内置的parula或者各类科学色带比如cmocean系列中的balance色带。balance是蓝白红的对称色带非常适合InSAR正负形变的表达而且是白底零值居中视觉上非常干净。cmocean并不是MATLAB官方内置的需要单独下载但我强烈建议纳入你的工具库。colormap(cmocean(balance));如果没有安装cmocean也可以用MATLAB推荐的线性色带构建方法cmap [0.00 0.00 0.55; 0.10 0.15 0.70; 0.25 0.45 0.90; 0.60 0.80 1.00; 1.00 1.00 1.00; 1.00 0.85 0.60; 1.00 0.60 0.30; 0.90 0.30 0.10; 0.55 0.00 0.00]; colormap(cmap);这种手写色带的逻辑是从深蓝到浅蓝到白色再到浅红到深红中间白表示零形变两侧对称。颜色数量和每段的RGB值可以按需增加但核心原则是——色带两端要深、中间要浅、深浅过渡要均匀避免某一段颜色跳变特别陡峭。4.2 NaN值的显示策略透明、补零还是插值很多人在caxis和NaN之间纠结。实际上m_pcolor配合pcolor时NaN值默认是不绘制的也就是像素点保持透明背后显示的是白纸或岸线的底图颜色。这个行为大多数情况下是理想的但如果你的研究区域内失相干区域占比很大图面上就会出现大片的空白反而分散读者注意力。处理策略要看具体场景如果是小范围的局部失相干河流、植被带保持透明即可读者能自然理解那是相干性差的区域如果失相干面积太大导致形变场看起来“稀疏”一般建议在图中用半透明的灰色多边形遮盖这些区域而不是把速率补成0。补0是最大的错误因为0在色标上对应的颜色是有物理意义“零形变”的会把失相干区域和真实零形变区域混淆。如果想在图上明确标记失相干区可以先把NaN区域提取出来再叠加一层灰色mask isnan(velocity); hold on; m_pcolor(LON, LAT, double(mask)); set(gca, Color, [0.7 0.7 0.7]); % 或者给mask矩阵设置特殊的colormap不过这个方案实施起来要注意colormap冲突同一张图只能有一套colormap灰色mask会和速率色带发生冲突。所以最终的折中方案是保持NaN透明通过论文文字说明“剔除相干性低于0.3的像素”这个方法简单直接审稿人也都接受。4.3 裁剪范围的经验法则让图面焦点集中但保留地理参考研究区域的裁剪范围会直接影响图的质量。范围太小周围缺乏地理参照物读者很难定位范围太大形变信号只占图面一小块有效信息密度低。我的经验法则是以形变区为中心外扩至包含至少两个以上有辨识度的地理要素。比如研究区附近有河流、湖泊、主要道路、大型城镇裁剪范围至少要把这些要素框进来。这样读者看图时可以通过这些熟悉的地理位置快速定位而不是盯着经纬度标注算坐标。具体实现时直接在m_proj里调整边界lonlim [lon_min - 0.3, lon_max 0.3]; latlim [lat_min - 0.3, lat_max 0.3]; m_proj(lambert, long, lonlim, lat, latlim);扩展多少取决于研究区尺度。一个经验参数是当研究区宽度约50公里时四周各外扩10~15公里宽度100公里以上时外扩20~30公里。外扩尺寸以经纬度表述时需要根据当地纬度做换算用deg2km函数辅助判断。5. 实测中躲不开的几个坑以及怎么绕过去这一部分写的是我这几年来用m_map画InSAR图时真实踩过的坑有一个算一个每个都让人花过不少冤枉时间。5.1 报错“index exceeds matrix dimensions”问题出在数据维度而不是代码遇到过好几次网上搜了一圈大部分答案都是让检查数组长度。真正的原因在InSAR场景里往往是这样GeoTIFF读取时velocity的维度是[rows, cols]但你生成的LON和LAT网格如果用了meshgrid输出的维度是[cols, rows]两者维度一致但顺序反了传进m_pcolor就直接报错或者画出旋转了90度的图。这个问题的排查方法是在绘图前用size检查所有变量的维度并强制统一[rows, cols] size(velocity); if size(LON, 1) cols size(LON, 2) rows LON LON; LAT LAT; end这行代码我几乎写进了每一个绘图脚本里。在没有严格约定时ROW变COL, COL变ROW是InSAR数据从外部软件导入MATLAB最常见的问题。5.2 图形显示正常但海岸线位置和形变偏移必须检查m_proj范围另一个让人困惑的问题是形变数据画出来位置正确但加上海岸线后海岸线和形变区明显错位几百米甚至几公里。这个问题的根源往往是数据坐标参考系与m_map默认坐标系不一致。举个我实际遇到的例子某次处理某区域数据处理软件导出的GeoTIFF虽然写的是WGS84但实际上是经过特定椭球拟合的局部坐标底图海岸线来自全球数据集直接用会有约500米的系统性偏移。解决方案是在GIS软件里检查数据是否正确或者用几个已知地面控制点道路交叉口、桥梁对比形变图与光学影像确认偏移量后手动平移数据。简单说在m_map里海岸线偏移不是绘图函数的问题而是数据坐标参考系没校准的问题。这一点在开始处理大批量数据前一定要确认否则后面对比验证会很难进行。5.3 色标上限、caxis、colormap三者关系混乱导致图面颜色失真画出来的图颜色和caxis设置对不上是新手最常见的问题。核心原因是MATLAB绘图的颜色映射逻辑每个像素点的颜色不是由它的数值直接决定的而是由“数值在caxis范围里的相对位置”与colormap的映射关系决定的。举一个典型错误colormap(jet); caxis([-30 30]); h m_pcolor(LON, LAT, velocity); % 先设置caxis再绘图如果先设置了caxis再绘图某些情况下图形对象不会立即更新导致画面颜色仍然按默认范围映射。正确顺序是先绘图再设置caxish m_pcolor(LON, LAT, velocity); set(h, EdgeColor, none); caxis([-vmax vmax]); colormap(cmap);另一个问题是同时设置多个子图时MATLAB会把所有子图共享同一套caxis。要避免这种情况每一张子图都要单独绘图并设置自己的caxis或者用caxis([-统一值 统一值])强制所有子图使用同一套颜色区间这在多图对比里反而是推荐做法因为色标不一致的对比图会误导读者。5.4 保存图片时文字变成方块或乱码MATLAB绘图完保存成PDF或EPS时中文标注乱码是老问题。InSAR图里的中文标注主要是研究区名称、图例文字虽然SCI期刊大多要求英文标注但很多中文论文或者阶段汇报仍然需要中文。我的做法是两点第一尽量用英文或拼音标注。这不是崇洋媚外而是从实际成图效果出发英文标注在任何格式下都不会有字体问题。第二如果一定要用中文保存使用高分辨率PNG不要用PDF/EPS矢量格式。因为MATLAB矢量输出时对中文字体的嵌入支持非常差换成PNG格式后文字是按像素输出的基本不会乱码。print(gcf, velocity_map, -dpng, -r600);600dpi的PNG用于论文完全够用于PPT汇报更是绰绰有余。5.5 多子图共用色标导致颜色错误叠加如果一张图里有多个子图比如LOS速率图、累积形变图、相干图一排放三张最容易出的问题是绘图顺序和colormap相互污染。正确做法是逐个绘图每一个子图进行完整的初始化、绘图、设置caxis和colormap流程。MATLAB的colormap是图形级别的属性一个Figure里有多个axes时colormap是全局共享的除非用colormap(ax, cmap)指定ax。subplot(1,3,1); m_pcolor(LON, LAT, velocity); caxis([-vmax vmax]); colormap(cmocean(balance)); colorbar; subplot(1,3,2); m_pcolor(LON, LAT, velocity_accum); caxis([-amax amax]); colormap(cmocean(balance)); colorbar;注意这里每次subplot后都要重新设置colormap不然第二张图会沿用上一张图的色带导致颜色错乱。6. 把m_map工具链融入日常InSAR出图流程的几点心得m_map并不是一个功能单一的小函数它本质上是一个画地图的框架。当你把它纳入自己的出图流程后效率提升是很明显的。我现在的标准流程是写完一套绘图脚本后后续每个研究区的出图只需要替换数据和边界坐标其他代码基本不动。脚本里包含了完整的投影初始化、岸线绘制、数据掩膜处理、填色绘图、要素叠加、色标调整和保存导出一个项目从头到尾的形变图质量都保持一致。另外m_map还支持将多个不同的数据集叠加到同一张底图上这就为InSAR结果和地质背景的综合分析提供了很大的便利。比如在研究某个滑坡体时我可以在同一张图上叠加速率场、断裂带位置、已有滑坡边界和GPS位移矢量信息层次分明审稿人看图就能直接抓住要点。如果研究区域比较大或者你需要更精细的底图数据也可以考虑和其他工具配合使用。m_map能加载用户自定义的高分辨率海岸线数据或DEM晕渲图配合上自己做的坡度、地形阴影出图效果会更接近专业制图软件的成品。有一点我特别想提醒m_map用的海岸线数据本身是全球范围内的中等精度数据在一些近岸精细化研究里岸线偏差可能达到数百米。如果论文中有严格的地表要素对比分析建议用更高精度的岸线或直接用光学影像叠加验证不要完全依赖m_map自带的海岸线。还有一个小技巧就是给m_map绘制的图设置统一的输出尺寸。很多人在MATLAB的figure窗口里直接保存结果拉出来的图比例时宽时窄放在论文里插图大小不一后期还得在排版软件里反复调整。我一般在绘图开始时用set(gcf, PaperPositionMode, manual); set(gcf, PaperUnits, centimeters); set(gcf, PaperPosition, [0 0 18 15]);这样可以保证每次输出的图尺寸一致后续排版省很多事。最后再说一点个人体会InSAR数据处理到最后图件往往是决定成果表达质量的临门一脚。数据解算得再好图件表达不清楚会大大降低成果的说服力。m_map虽然不是最现代的绘图方案但在MATLAB社区里成熟稳定、教程多、排错容易对InSAR研究者来说仍然是最稳妥的选择。把这套流程练熟你会发现出图这件事本身也可以变成一件让人有成就感的事。