开头不用我多解释Anusplin这个名字搞过气象插值的人应该都有所耳闻。我最早是在简书上看到有人一步步演示Anusplin插值的完整流程那篇文章帮我打开了门但等我自己动手把“全省几百个自动气象站的月平均气温变成一张连续分布图”这件事真正跑通时才发现中间有大量细节是原帖一笔带过的。作为一个GIS基础并不扎实、Fortran和命令行都只停留在“听说过”层面的小白整个重现过程可以说是一路踩坑一路补课。这篇文章就是把我从零开始逐步实现Anusplin气象插值的全过程、所有配置、所有坑都记录下来给同样被老派工具折磨的人做个参考也欢迎大家指正。1. 为什么是Anusplin薄盘样条与IDW、克里金的本质差别1.1 薄盘光滑样条在气象插值里为什么能打Anusplin是澳大利亚国立大学的M. F. Hutchinson团队开发的薄盘光滑样条插值程序集最早可以追溯到上世纪80年代。它广泛应用于气象、水文、生态等领域的气候要素空间化尤其是气温、降水、辐射这类受地形影响明显的要素。它的核心思路可以拆成两部分一部分是“常规趋势”另一部分是“空间光滑残差”。通俗点说Anusplin认为某个站点的气温大致等于“高程的函数比如海拔每升高100米气温降多少加上一个随空间位置变化的光滑调整项”。第一项描述了气象要素随地形变化的宏观规律第二项负责吸收局地小气候、观测误差等偏离趋势的部分。这个设计正好击中了气象插值的痛点。气温随海拔递减是一个普遍规律降水和地形的关系虽然更复杂但也有迹可循。假如只用纯空间插值方法比如IDW你插出来的是“水平距离上的加权平均”地形信息完全用不上而Anusplin把DEM作为协变量放进去相当于在插值模型里引入了垂直方向的解释力。为什么叫“薄盘光滑样条”可以这样理解想象一块很薄的金属板要穿过这些站点观测值并且板面尽量平滑。样条的意义在于它不强求函数刚好穿过每个观测点而是允许有一定偏差偏差大小由“光滑参数”控制。这个光滑参数怎么选最合适Anusplin用广义交叉验证GCV自动确定这也是它比手动调参的插值方法更让人省心的地方。1.2 同一组数据三种插值方法差异直观对比我用同一份某省2019年7月平均气温数据大概400多个站点分别做了IDW、克里金和Anusplin插值差异非常直观IDW图上全是“牛眼”一个站点就是一个圈山区等温线像被刀切过一样毫无地形逻辑。普通克里金斑块感比IDW好一些但同样没有把高程信息引进去山区和平原的过渡很生硬。Anusplin加入DEM协变量之后山区的温度明显随海拔降低而递减等温线走向和山脉走向大体一致视觉上“顺眼”太多了。这种差异在平原地区可能不明显但只要有山地、高原Anusplin的优势就会非常突出。很多公开发表的气候区划、生态模型研究里气温和降水栅格数据就是用Anusplin生成的这不是没有原因的。对于GIS小白来说我的建议很直接如果你手头有DEM要插值的气象要素又明显受地形影响那Anusplin值得花一个下午去折腾如果只是做平原地区、要素变化很平缓的简单展示用IDW或克里金就够没必要上这个老派工具。2. 首次安装Windows与Linux两条路线怎么选2.1 官网下载、注册与许可协议Anusplin的官方网站是澳大利亚国立大学Fenner School环境与社会学院维护的页面直接搜“Anusplin”基本能第一个找到。网站要求先填写一份申请/注册信息包括姓名、单位、邮箱、用途等提交后会收到下载链接和许可协议。有一点需要特别注意Anusplin是非商业用途免费的软件它的许可协议明确要求在使用时引用相关文献并且不能把程序嵌入商业系统牟利。如果你是在学校或者科研院所工作正常提交申请即可通常一两个工作日就能收到下载链接。下载下来之后文件包里一般包含程序文件、示例数据、一个很老的PDF格式用户手册以及一些示例命令文件。这个手册虽然排版简陋但信息量非常大后面所有配置参数的说明都以它为准。建议下载后第一时间解压把手册打印一份或者在另一个屏幕上常开随时翻。2.2 Linux下源码编译的关键步骤我自己的主力工作机是Linux所以先选了源码编译路线。Anusplin是用Fortran写的编译依赖gfortran。环境准备好之后进入源码目录执行make命令即可。这里只说我实测有效的流程sudo apt install gfortran tar -zxvf anusplin.tar.gz cd anusplin make编译完成之后目录下会生成splina、splinb、splinc等可执行文件。整个过程如果不出意外几分钟就能完成。如果make报错大多数情况是缺少某些系统库比如libgfortran装一下就行sudo apt install libgfortran5需要注意不同小版本的源代码目录结构和makefile可能略有差异遇到报错先看官方手册的“Installation”章节里面写得很清楚。不要一报错就怀疑软件有问题先把依赖装齐成功率很高。2.3 Windows下的运行思路如果你只有Windows机器也不用慌。官网同时提供Windows版本的预编译程序下载解压之后目录里就是可以直接运行的splina.exe等文件。不过Anusplin是纯命令行程序需要在CMD里切到对应目录执行操作方式对所有用惯ArcGIS图形界面的GIS用户来说第一次可能会觉得“怎么这么原始”。另一个更顺畅的玩法是在Windows下用WSLWindows Subsystem for Linux直接装一个Ubuntu子系统然后在WSL里按Linux方式编译和运行。这样既能保留Windows里的ArcGIS/QGIS做后处理又能有一个干净的Linux环境跑Anusplin。我的经验是不要试图把Anusplin“安装”到系统路径里直接把它放在一个固定文件夹比如D:\anusplin\或者~/anusplin/然后每次在这个目录下操作。它不需要注册表、不需要环境变量就是个古老的可执行文件反而简单。3. 站点数据与DEM协变量格式对了才谈得上插值3.1 站点文本表第一行到底写什么Anusplin对输入数据格式的执着是我在整个过程中印象最深的部分。它读取的不是Excel也不是Shapefile而是特定排列的纯文本文件。站点数据文件的基本格式是第一行写站点数量从第二行开始每一行代表一个站点依次是经度、纬度、高程、观测值。以我用的七月平均气温数据为例文件名temp_july_station.txt内容大致长这样412 118.12 32.48 45.2 26.8 118.38 32.62 88.6 25.9 119.01 32.74 62.4 27.1 ...这里有几个细节非常容易出错经度和纬度的单位必须是十进制度不能用度分秒。如果原始数据是度分秒需要先统一转换。高程和观测值之间用空格分隔不要用逗号也不要用制表符混排。Anusplin的读取逻辑很死板老老实实用空格隔开最稳妥。站点数量这一行不能省略也不能写成其他说明文字。我一开始想当然地写成了“站点数412”程序直接报错连解析阶段都过不去。观测值如果有缺测不要留空也不要填-9999Anusplin对缺测有自己的处理字段。小白阶段最稳妥的办法是直接用完整站点数据不缺测的月份先行筛选。如果你的数据文件里不止一个气象要素而是同时有气温、降水、辐射等多列也是支持的只要在后续配置里把变量数量写对。但第一次跑通流程强烈建议只放一个要素减少变量。3.2 DEM协变量范围、分辨率与ASCII导出Anusplin的另一个核心输入是协变量栅格文件通常就是DEM。它要求的是ASCII Grid格式也就是和ArcGIS里“栅格转ASCII”工具导出的.asc文件一致。在ArcGIS里操作路径是加载DEM图层右键导出选择“栅格转ASCII”。导出前一定要确认DEM的坐标系是GCS_WGS_1984经纬度坐标。如果原始DEM是其他坐标系比如UTM投影需要先用“投影栅格”工具转成WGS84再导出否则经纬度和高程混在一起插值结果会乱套。关于分辨率的设置我踩过一个大坑。最开始我直接用原始30米分辨率的DEM导出文件巨大Anusplin跑了一天一夜都没出结果。后来查手册才知道Anusplin虽然是栅格插值工具但它的网格点数是有上限的过高的分辨率除了拖慢速度对插值结果并没有本质提升。一般省级尺度的气象插值用0.05度约5公里分辨率的DEM就足够了如果是小范围研究区可以用0.01度甚至更高。重采样方法选双线性或最邻近均可我一般选双线性地形过渡更自然。还有非常重要的一条原则DEM的覆盖范围最好比研究区边界往外扩一圈至少要包含所有参与插值的站点并且给出一定余量。原因后面会详细说这里先记住一句话Anusplin是在“整个DEM覆盖范围”上生成连续曲面不是按你的研究区边界去裁剪的。3.3 关于投影Anusplin只认经纬度这一点值得单独拿出来强调。Anusplin计算距离和空间关系的坐标系默认是经纬度。如果你把站点坐标转成了投影坐标比如高斯-克吕格又把DEM也转成了投影坐标理论上Anusplin也能算但输出的网格坐标就是投影坐标后续和通用GIS数据叠加容易出问题。更麻烦的是投影坐标系的单位是米经纬度坐标的单位是度二者的“距离”语义完全不同。Anusplin在计算样条光滑程度时默认的误差方差比参数是依据经纬度坐标设计的如果输入坐标单位变了这个参数也需要重新调整小白阶段完全没必要给自己加这个难度。所以我的建议非常明确站点坐标统一用十进制度经纬度DEM统一用GCS_WGS_1984导出ASCIIAnusplin计算完之后再在ArcGIS或者QGIS里把结果投影到目标坐标系。输出可以最后转输入不要乱转。4. 配置文件与命令行splina每行参数到底什么意思4.1 splina.cmd配置逐行解析Anusplin不是双击就运行的程序它需要你写一个文本配置文件把各种参数告诉它。这个配置文件没有固定扩展名官方示例里经常叫splina.cmd我沿用这个命名。以下是我跑通七月平均气温插值时的配置模板逐行说明它的含义# splina.cmd 示例配置 temp_july_station.txt dem_0.05.asc temp_july_out 1 2 1 1.0 0逐行解析第一行站点数据文件名。注意这个路径是“相对于运行命令时所在目录”的最好用相对路径减少出错。第二行协变量栅格文件名DEM的ASCII文件。第三行输出文件的前缀。程序会生成temp_july_out.sur、temp_july_out.dat、temp_july_out.log等一系列文件。第四行变量数量。指站点文件里需要插值的气象要素个数填1。第五行样条次数。填2代表三次样条这是最常用的设置适合大多数气象要素。第六行协变量数量。填1表示只有DEM一个协变量。第七行误差方差比。填1.0是官方推荐的默认值表示先验误差方差和样条光滑方差相等。如果你对数据的测量误差有更精确的认识可以改但小白阶段不要动。第八行是否输出诊断信息。填1会输出更多交叉验证细节填0则只保留精简日志。我建议第一次跑填1可以多看到一些内部信息。需要注意不同版本的Anusplin命令行文件里的参数数量和顺序可能会有一点点出入所以拿到软件包之后先打开自带的示例文件对照一下再填千万不要凭记忆硬套。我吃过的亏就在这里网上找了一个老版本的配置模板直接套到新版本上结果程序识别不了。4.2 变量数、协变量数与样条次数怎么定这三个参数是配置文件里最容易让人懵的我详细说一下我的理解。变量数量NUMBER OF VARIABLES指的是你希望同时插值多少个气象要素。如果站点文件里只有一列气温数据那就填1。如果你把气温和降水两列都放在了站点文件里想一次插值两个要素就填2。多变量模式在模型结构上更复杂输出结果里会包含两套曲面。第一次跑通流程建议只填1。协变量数量NUMBER OF COVARIATES是模型里除了经纬度位置之外参与解释空间趋势的自变量个数。填1代表只用高程一个协变量DEM。Anusplin支持多个协变量比如高程加纬度、高程加距海岸距离这属于进阶玩法后面可以扩展。对于小白1就够用这也是Anusplin相对克里金的优势所在——它把“加协变量”这件事简化成了一个数字。样条次数ORDER OF SPLINE填2是三次样条这是气象插值中最常用的选择。填3、4理论上可以有更复杂的光滑曲面但对站点数据量的要求也会更高容易过度拟合。老老实实填2不要作。4.3 运行命令与日志输出怎么看在Linux下运行命令非常简单./splina splina.cmd在Windows下则是splina.exe splina.cmd运行之后终端会刷出一堆调试信息同时工作目录下会生成前缀为temp_july_out的一系列文件。程序运行时间取决于DEM网格点数和站点数量我的0.05度分辨率、全省范围的数据几分钟就能跑完如果用高分辨率DEM可能跑几小时甚至更久。跑完之后最重要的事情是打开.log文件重点看两个指标GCV分数和误差方差。GCV分数越小说明模型的广义交叉验证误差越小拟合与泛化的平衡越好。误差方差如果出现异常值比如比观测值本身的方差还大说明模型可能没收敛或者输入数据有问题。我在第一次跑通时GCV数值大概在个位数误差方差也在合理范围心里才踏实。如果看到日志里出现“FAILED TO CONVERGE”或者“NEGATIVE VARIANCE”之类的字样不要慌绝大多数情况不是参数问题而是数据问题比如站点数量太少、坐标越界、DEM范围没覆盖站点。先检查输入数据而不是去调参数。5. 结果文件解读与GIS出图从.dat到栅格图5.1 输出文件里哪些才是真正需要的Anusplin运行结束后会生成一大堆后缀不同的文件刚看到的时候很容易懵。我实际用下来的文件清单大致如下文件后缀内容用途.log运行日志含GCV、误差方差、迭代过程判断模型是否收敛、质量如何.sur拟合曲面系数文件程序内部使用一般不需要手动处理.dat插值结果的ASCII格网文件核心结果用于转栅格.res残差与站点诊断文件交叉验证、站点残差分析.cov协变量系数文件查看高程协变量的系数估计.lis列表输出备查其中最重要的是.dat文件它就是你插值结果的可视化前身。这个文件的格式和DEM的ASCII文件类似有行列数、坐标范围、像元大小等头部信息后面跟着一排排数值。ArcGIS可以直接识别这种格式。.sur和.cov这些文件我基本不打开但建议保留有时候需要复核协变量系数比如想确认“气温随高程递减率到底是多少度每千米”就可以去.cov文件里找到这个系数。对于纯做图的人来说有了.dat就够了。5.2 ASCII转栅格与投影转换把.dat转成真正能显示的栅格这个步骤和转DEM完全一样。在ArcGIS里用“ASCII转栅格”工具输入文件选temp_july_out.dat输出栅格选temp_july_out_grd数据类型选FLOAT观测值通常是浮点数不要选INTEGER。转换完成之后栅格是没有定义坐标系的状态需要手动给它指定坐标系。因为在Anusplin里输入的就是GCS_WGS_1984的经纬度所以这里直接给栅格定义成GCS_WGS_1984即可。操作路径是“定义投影”工具选Geographic Coordinate Systems里的WGS 1984。如果你习惯用QGIS也可以用GDAL命令行一步完成gdal_translate -of GTiff temp_july_out.dat temp_july_out.tifGDAL会自动识别ASCII Grid格式转出来的TIFF默认就是经纬度坐标后面再用gdalwarp做投影转换。整个过程灵活度比ArcGIS更高但小白用ArcGIS的图形界面更直观。投影转换这一步我会在出图前做。如果我的研究区用阿尔伯斯等积投影或者UTM就在“投影栅格”工具里设置目标坐标系。转换时注意选择重采样方法为双线性因为气温场是连续变量最邻近重采样会产生块状瑕疵。5.3 用研究区边界掩膜提取与验证投影完之后整张图的范围还是DEM的范围通常比研究区大一圈。这时候用研究区边界去做“按掩膜提取”把研究区外的区域去掉出图就会干净很多。操作路径是ArcToolbox里的“按掩膜提取”工具输入栅格选投影后的气温栅格掩膜数据选研究区边界面文件。这里有一个常见问题如果掩膜面太碎、包含很多空洞提取结果可能会出现一些细小的无值缝隙。解决办法是把行政区边界先做一次消除去掉小的飞地再用于掩膜。做完掩膜之后不要急着出图先做一次简单验证。把原始站点加载进来用“多值提取至点”工具提取每个站点位置的插值结果然后打开属性表新建一个字段计算“插值减实测”的差值。看一下这个差值的平均值和范围如果在1到2摄氏度以内说明插值结果整体可接受如果出现十几度的偏差就要回头检查站点数据或者配置参数了。我用这个方法验证过七月平均气温大部分站点残差在正负1度左右山区个别站点偏差稍大整体是可以放入报告的水平。6. 我踩过的坑从深山到沟里一步一步找原因6.1 坑1DEM范围比研究区小边界插值全是“布丁”这是我最开始犯的错。我用ArcGIS把DEM裁剪到了研究区边界认为这样更“精确”结果Anusplin跑完之后研究区边界内侧出现了一大圈颜色特别突兀的斑块像巧克力蛋糕上的布丁怎么都解释不通。后来查手册、看日志才明白Anusplin是在整个DEM覆盖范围上生成曲面的边界外的区域没有协变量信息样条就会“放飞自我”产生剧烈的震荡。这种震荡会从边界向内扩散污染到研究区边缘几十公里的范围。解决办法就是前面说的DEM范围一定要大于研究区边界至少外扩半个到一个网格单元的跨度最好直接下载整个DEM数据块覆盖研究区及其周边所有站点。最后再在GIS里用研究区边界掩膜提取这样边界效应就会落在研究区之外出来的图干干净净。6.2 坑2站点表第一行写成列名程序直接读挂这个错误说起来很丢人但可能很多小白都会犯。我一开始用Excel整理数据习惯了第一行写表头导出成文本之后就变成了lon lat elev temp 118.12 32.48 45.2 26.8 ...Anusplin把第一行当成站点数量来读读出一个字符串程序直接报错。那时候我还以为是软件装错了重装了两遍才意识到是数据格式问题。正确做法是用文本编辑器Notepad、VS Code都行打开站点数据文件手动检查第一行是否就是纯数字的站点数量如果不是就改掉。处理Excel导出的数据时建议导出为“空格分隔的文本文件.prn或者.txt”而不是CSVCSV里的逗号会让坐标和观测值挤成一个大字段。6.3 坑3ASCII转栅格后和DEM对不齐有一段时间我出的图气温栅格和DEM栅格叠在一起总是错开半个像元边界像锯齿一样对不上。原因在于我输出栅格时没有设置和DEM一致的“捕捉栅格”环境。在ArcGIS里“ASCII转栅格”工具的默认输出范围就是文件自带的坐标范围理论上应该和DEM一致。但如果你之前对DEM做过重采样、裁剪、投影转换它的范围和像元对齐状态可能已经和原始ASCII文件不完全相同了。解决办法是在ArcToolbox环境设置里找到“栅格分析”把“捕捉栅格”设成DEM图层再把“像元大小”设成和DEM一致然后再执行ASCII转栅格。这样输出栅格就会自动对齐到DEM的像元网格上后续任何叠加分析都不会出现错位。6.4 坑4Windows记事本中文路径命令文件找不到文件这个坑主要发生在Windows平台上。我用记事本创建splina.cmd文件里面填了输出文件前缀并且把整个工程放在D:\气象插值\2024年7月\这种带中文和空格路径的文件夹里。结果运行splina.exe的时候程序各种找不到文件日志里的路径还是乱码。原因有两层一是记事本保存文本文件默认用的是ANSI编码或者带BOM的UTF-8Anusplin这个老程序对带BOM的文件头非常敏感它会把BOM当成第一个字符读进去导致文件名解析失败二是程序内部对非ASCII字符路径支持很差中文路径直接乱码。解决办法很粗暴所有Anusplin相关文件统一放在纯英文路径下比如D:\anusplin\july\文件名也全部用英文小写。配置文件用Notepad或者VS Code保存为UTF-8无BOM格式。这个准备工作只要一次做对能省好几个小时。6.5 坑5降水插值出现负值要不要截断用Anusplin插值气温很顺利但换到降水数据时就出现了问题插值结果在局部区域出现负值这在气象上完全不合理。我最初的想法是直接把这个负值重分类成0但后来觉得这样做太粗暴容易在研究区留下一个突然的“断崖”。查了文献才发现Anusplin的样条模型默认假设数据是正态分布的而降水数据通常偏态分布大量站点为0或很小值直接插值就会在0值区域产生负值。官方手册里给的建议是对数据进行Box-Cox变换或者开方变换让数据分布更接近正态再插值最后把结果变换回来。我实际操作下来用四次方根变换效果比较好插值前对每个站点的降水值开四次方插值完成后对栅格做四次方运算还原。这样得到的降水分布不会出现负值而且0值区域的过渡也自然很多。如果你的降水数据里有大量0值这是一个非常值得记住的技巧。写到这里Anusplin从下载、安装、数据准备、配置运行到GIS出图的完整链路基本就闭环了。这套流程我前前后后跑了十几轮每次被某个细节卡住的时候都会想起简书那篇带我入门的大神文章——但也正因如此我才更想把那些原帖里没展开讲的“坑”和“为什么”补全。如果你也是GIS小白刚开始碰Anusplin我唯一强烈的建议是先拿一个范围不大、站点不多、只有单个要素的假数据把整条流程跑通再上真实数据每次改配置之前把上一个能跑通的splina.cmd备份一份。老派工具不复杂但它对你的耐心要求比现代软件高得多。