一、栅格数据基本概念
1.1 什么是栅格数据
栅格数据(Raster Data)是GIS中两大空间数据模型之一，它以规则排列的像元(Pixel/Cell)阵列来表达空间信息。与矢量数据使用离散的点、线、面不同，栅格数据将整个研究区域划分为均匀的网格单元，每个网格单元存储一个或多个属性值。这种数据结构天然适合表达连续变化的地理现象，如地表温度、高程、降雨量、植被覆盖度等。
从直觉上看，栅格数据就像一张数字化的照片 —— 照片由无数像素组成，每个像素有其颜色值；同理，栅格数据由无数像元组成，每个像元记录着对应地表位置的某种属性值。卫星遥感影像、数字高程模型(DEM)、扫描地图、土地利用分类图等都属于典型的栅格数据。
1.2 像元(Pixel/Cell)与行列结构
像元是栅格数据的最小组成单位。在GIS中，"Pixel"和"Cell"这两个术语经常混用，但严格来说存在细微区别：Pixel多用于遥感影像语境，指传感器采集的一个采样单元；Cell则更常见于GIS分析语境，指栅格网格中的一个单元格。
栅格数据的行列结构(Row-Column Structure)遵循矩阵的组织方式。整个栅格被组织为M行 × N列的二维矩阵。行号从顶部开始向下递增(从0或1开始)，列号从左侧开始向右递增。每个像元可以通过其行列号 (row, col) 唯一确定位置。例如，一幅1000×1000的遥感影像包含100万个像元。
像元属性的数据类型
像元值可以是多种数据类型：整型(Integer) — 用于分类数据如土地利用类型(1=农田, 2=森林, 3=水体)；浮点型(Float) — 用于连续数据如高程(1234.56m)或温度(23.7°C)；布尔型(Boolean) — 用于二值数据如洪水淹没区(0/1)。单波段栅格每个像元存一个值，多波段(Multi-band)栅格每个像元可同时存储多个波段值。
1.3 栅格坐标系与空间参考
栅格数据具有双重坐标系统。图像坐标系(Image Coordinate System)以像元的行列号定位，原点通常位于左上角，行号向下增长，列号向右增长。地理坐标系/投影坐标系(Geographic/Projected Coordinate System)则将每个像元映射到真实的地球表面位置。
这种映射关系通过仿射变换(Affine Transformation)来描述，包含六个参数：
X_geo = X_origin + col × PixelWidth + row × RotationX
Y_geo = Y_origin + col × RotationY + row × PixelHeight
其中：
X_origin, Y_origin = 栅格左上角的地理坐标
PixelWidth, PixelHeight = 像元在X和Y方向上的地面尺寸(PixelHeight通常为负值)
RotationX, RotationY = 旋转参数(无旋转时为0)
在GeoTIFF文件中，这些参数存储在TIFF的GeoKey标签中；在ESRI格式中，则保存在同名的.tfw(World File)或.prj(投影文件)中。
1.4 NoData概念
NoData是栅格数据中一个关键概念，表示"此像元没有有效数据"。它不同于数值0 — 0是一个合法的数据值(例如高程为0表示海平面)，而NoData意味着数据缺失或不适用。常见的NoData设定值包括 -9999、-3.4e38、255(8位无符号整型)等。
NoData产生的原因包括：传感器故障导致的坏像元(Bad Pixel)、云层遮挡无法获取地面信息、影像裁剪后的不规则边界区域、投影变换导致的空值区域等。在空间分析中，NoData像元通常不参与计算 — 例如计算区域平均高程时，NoData像元会被自动排除，而非以0值参与统计，这对于保证分析结果的准确性至关重要。
二、四种分辨率详解
栅格数据的"分辨率"并非单一概念，而是包含空间分辨率、光谱分辨率、时间分辨率和辐射分辨率四个维度。这四种分辨率共同决定了栅格数据的信息承载能力和适用场景。
2.1 空间分辨率(Spatial Resolution)
空间分辨率指一个像元所对应的地面实际面积，通常以地面采样距离 GSD (Ground Sample Distance)来衡量。GSD越小，空间分辨率越高，能够分辨的地物越细小。例如，0.3m分辨率的卫星影像可以分辨停车场中的单辆汽车，而30m分辨率的影像只能区分不同的土地利用类型。
需要注意，空间分辨率与传感器的瞬时视场角 IFOV (Instantaneous Field of View)和轨道高度直接相关。IFOV越小、轨道越低，空间分辨率越高，但覆盖范围(Swath Width)越窄。这是遥感系统设计中的一个基本矛盾。
主流遥感卫星空间分辨率对比
卫星/传感器国家全色(Pan)多光谱(MS)幅宽(km)发射年份
WorldView-3美国0.31 m1.24 m13.12014
Pléiades Neo法国0.30 m1.2 m142021
GeoEye-1美国0.41 m1.65 m15.22008
高分二号(GF-2)中国0.8 m3.2 m452014
SPOT-7法国1.5 m6 m602014
高分一号(GF-1)中国2 m8 m602013
Sentinel-2ESA—10 m2902015
Landsat 8/9 OLI美国15 m30 m1852013/2021
MODIS (Terra/Aqua)美国—250/500/1000 m23301999/2002
参考: 空间分辨率的选择取决于应用目标。城市精细化管理需要亚米级影像(< 1m)，农业监测适用10-30m级别(Sentinel-2, Landsat)，全球尺度植被监测则使用250m-1km级别(MODIS)。分辨率越高，数据量和处理成本呈几何级数增长。
2.2 光谱分辨率(Spectral Resolution)
光谱分辨率描述传感器区分不同电磁波长的能力，由波段数量和波段宽度共同决定。根据光谱分辨率的不同，遥感影像可分为以下几类：
全色影像 (Panchromatic)
单波段，覆盖可见光范围(约0.45-0.90μm)。空间分辨率最高，但无法区分颜色。常用于与多光谱影像融合(Pan-sharpening)以提高空间细节。
多光谱影像 (Multispectral)
3-10个波段，通常包括蓝(B)、绿(G)、红(R)、近红外(NIR)。Landsat有11个波段，Sentinel-2有13个波段。波段宽度较宽(数十至上百nm)，适合地物分类和植被监测。
高光谱影像 (Hyperspectral)
数十至数百个连续窄波段(波段宽度5-10nm)。如EO-1 Hyperion有242个波段、PRISMA有240个波段。可获取近乎连续的光谱曲线，适合矿物识别、精细植被分类、水质监测等。
光谱分辨率的重要性在于：不同地物在不同波长处的反射率存在差异。例如，健康植被在近红外波段(NIR, ~0.7-1.3μm)反射率极高(约40-50%)，而在红光波段(Red, ~0.6-0.7μm)因叶绿素吸收反射率很低(约5-10%)。这种"红边效应(Red Edge)"正是NDVI等植被指数的物理基础。
2.3 时间分辨率(Temporal Resolution)
时间分辨率也称重访周期(Revisit Period)，指卫星重复观测同一地点的最短时间间隔。对于需要监测动态变化的应用（如农作物长势监测、洪水应急、森林火灾追踪），高时间分辨率至关重要。
卫星重访周期应用场景
MODIS (Terra+Aqua)1-2天全球植被变化、火点监测、海洋温度
Sentinel-2 (双星)5天农业监测、土地变化检测
Landsat 8+9 (双星)8天土地利用/覆盖制图、长时间序列分析
WorldView-3< 1天 (侧摆)应急响应、军事侦察
Planet (SkySat星座)每天日变化监测、精准农业
风云四号(FY-4A)15分钟 (静止轨道)天气预报、台风追踪
时间分辨率与空间分辨率之间通常存在权衡关系(Trade-off)：MODIS每天覆盖全球但空间分辨率仅250m-1km，而WorldView-3空间分辨率0.31m但重访周期较长。近年来的小卫星星座(如Planet Labs拥有200+颗卫星)正在打破这一限制，实现"每天全球覆盖+米级分辨率"。
2.4 辐射分辨率(Radiometric Resolution)
辐射分辨率是指传感器区分不同辐射强度差异的能力，以位深度(Bit Depth)衡量。位深度决定了像元值可能的取值范围(灰度级数)：
位深度灰度级数取值范围典型传感器
8-bit256级0 - 255早期Landsat MSS、SPOT-1
11-bit2048级0 - 2047IKONOS、QuickBird
12-bit4096级0 - 4095Landsat 8 OLI、Sentinel-2
14-bit16384级0 - 16383WorldView-3
16-bit65536级0 - 65535某些高光谱传感器
更高的辐射分辨率意味着传感器能够捕捉更微小的辐射差异。例如，8-bit传感器只能将入射辐射量化为256个等级，而12-bit传感器可以量化为4096个等级。这对于区分低对比度地物(如不同含水量的土壤、不同密度的植被)至关重要。在暗区域和亮区域同时存在的场景中，高辐射分辨率能够避免信息丢失。
三、DEM / DSM / DTM 深入对比
3.1 三种高程模型的定义
在GIS和遥感领域，高程数据是最基础的栅格数据类型之一。三种高程模型经常被提及，但它们的含义存在重要差异：
DEM (Digital Elevation Model)
数字高程模型 — 广义上是所有数字化高程数据的统称。狭义上，特指表示裸露地表(Bare Earth)高程的栅格数据，不包含建筑物和植被。在美国USGS的定义中，DEM通常等同于DTM。
DSM (Digital Surface Model)
数字表面模型 — 记录地球表面最高点的高程，包括建筑物屋顶、树冠顶部、桥梁等人工和自然地物。可以理解为"鸟瞰视角"看到的表面高度。
DTM (Digital Terrain Model)
数字地形模型 — 仅记录自然地形表面的高程，去除了所有人工建筑物和植被。在欧洲传统定义中，DTM还可能包含坡度、坡向等附加地形属性。
nDSM：归一化数字表面模型
nDSM = DSM - DTM，得到的是地物相对于地面的高度，即建筑物高度或树木高度。这在城市三维建模和森林蓄积量估算中非常有用。例如，一栋建筑物的DSM高程为85m，对应位置的DTM高程为50m，则该建筑物高度nDSM = 85 - 50 = 35m。
3.2 高程数据的获取方法
摄影测量(Photogrammetry)
通过立体像对(Stereo Pair)利用视差(Parallax)原理计算地面点的三维坐标。传统航空摄影测量精度可达分米级。卫星立体摄影测量(如ALOS PRISM、WorldView立体对)精度一般在1-5m。近年来，无人机(UAV)摄影测量结合运动恢复结构(SfM, Structure from Motion)技术，可获取厘米级精度的DSM。
LiDAR (Light Detection and Ranging)
激光雷达通过发射激光脉冲并记录回波时间来计算距离。机载LiDAR(Airborne LiDAR)可同时获取首次回波(First Return)生成DSM和末次回波(Last Return)生成DTM。点云密度通常为2-20点/m²，高程精度优于15cm(垂直RMSE)。星载LiDAR如ICESat-2采用光子计数技术，虽然覆盖范围大但为剖面式而非面状覆盖。
InSAR (Interferometric Synthetic Aperture Radar)
干涉合成孔径雷达利用两次SAR成像的相位差计算地面高程。全球SRTM DEM就是2000年航天飞机雷达地形任务通过InSAR技术获取的。InSAR的优势在于可穿透云层(微波波段不受天气影响)，适合常年多云的热带地区。但InSAR获取的是DSM而非DTM(微波难以完全穿透茂密植被)。
3.3 全球DEM数据集对比
数据集分辨率覆盖范围垂直精度(RMSE)数据源发布年份
SRTM v330m (1弧秒)60°N - 56°S~9m (全球)InSAR2000 (2014开放30m)
ASTER GDEM v330m (1弧秒)83°N - 83°S~8.5m立体摄影测量2019
ALOS World 3D30m (1弧秒)全球~5mPRISM立体对2015
COP-DEM GLO-3030m全球< 4m (相对)TanDEM-X InSAR2021
COP-DEM GLO-9090m全球< 4m (相对)TanDEM-X InSAR2021
NASADEM30m60°N - 56°S~7mSRTM再处理2020
FABDEM30m全球改进的COP-DEMCOP-DEM去除建筑物+树木2022
注意: SRTM和TanDEM-X获取的是DSM而非DTM。在森林茂密地区，SRTM高程可能比真实地面高程偏高5-20m。FABDEM基于COP-DEM通过机器学习移除了建筑物和树木影响，更接近于DTM。
3.4 应用场景
水文分析：基于DEM提取流域边界、河网、汇水面积、流向(D8/D-infinity算法)
地形分析：计算坡度(Slope)、坡向(Aspect)、曲率(Curvature)、地形湿度指数(TWI)
视域分析(Viewshed)：从DEM计算特定观察点的可视范围，应用于通信基站选址、景观规划
洪水模拟：高精度DEM是洪水淹没模拟(如HEC-RAS 2D)的关键输入
城市三维建模：DSM + nDSM用于建筑物高度提取和LOD1/LOD2三维城市模型构建
地质灾害评估：坡度、曲率、地形起伏度等DEM衍生指标是滑坡易发性评估的重要因子
四、栅格数据编码结构
原始栅格数据采用逐像元存储的方式，数据量巨大。例如，一幅10000×10000的32位单波段栅格需要约381MB存储空间。为了提高存储效率和减少I/O开销，发展出了多种压缩编码结构。
4.1 直接编码(Direct Coding)
最简单的栅格存储方式：按行或按列逐像元记录每个像元值。数据结构就是一个二维矩阵，行数×列数即为像元总数。优点是结构简单、随机访问速度快(O(1)时间复杂度)，缺点是不进行任何压缩，存储空间开销最大。
存储空间 = 行数 × 列数 × 每像元字节数
例: 10000 × 10000 × 4 bytes(Float32) = 400,000,000 bytes ≈ 381 MB
4.2 游程编码(Run-Length Encoding, RLE)
游程编码是最经典的栅格压缩方法之一。其原理是：对于同一行中连续相同的像元值，只记录"(值, 连续长度)"这样的二元组，而非重复存储。
原始数据(一行): 3 3 3 3 5 5 7 7 7 7 7 7 2 2 2
RLE编码: (3,4) (5,2) (7,6) (2,3)
原始存储: 15个值
RLE存储: 4个二元组 = 8个值
压缩率: 8/15 ≈ 53%
RLE的压缩效果取决于栅格数据的空间自相关性：连续区域越大(如大片农田、湖泊)，压缩率越高；当数据高度碎片化(如分类结果中的椒盐噪声)时，RLE可能反而增加存储量。对于典型的土地利用分类图，RLE通常能实现60%-90%的压缩率。
4.3 链码(Chain Code / Freeman Code)
链码由Freeman于1961年提出，主要用于编码区域边界。其思想是：选择一个起始点，然后用方向代码(0-7，分别代表8个方向)依次记录边界走向。
Freeman 8方向编码：
3 2 1
4 · 0
5 6 7
0=东, 1=东北, 2=北, 3=西北, 4=西, 5=西南, 6=南, 7=东南
一个矩形区域的边界链码示例：
起点(2,1), 链码: 0 0 0 6 6 4 4 4 2 2
链码适合存储以边界为主的数据(如行政区划边界)，但不适合表达内部属性信息。其变体差分链码(Differential Chain Code)只记录方向的变化量，进一步提高了压缩效率。
4.4 块码(Block Code)
块码是游程编码在二维空间上的推广。不同于RLE只在行方向压缩，块码将栅格分解为尽可能大的矩形块，每个块用"(起始行, 起始列, 行数, 列数, 值)"来表示。例如，一大片相同值的区域可以用一个矩形块来表示，而非逐行RLE编码。块码的压缩率一般优于RLE，但编码和解码算法更复杂。
4.5 四叉树编码(Quadtree Encoding)
四叉树是栅格数据最重要的层次化压缩编码结构。其原理是递归四分法：将整个栅格区域等分为四个象限(NW, NE, SW, SE)；如果某个象限内所有像元值相同，则记录为叶节点；否则继续四分，直到所有子区域均匀或达到单像元。
四叉树递归分解过程：
Level 0: 整个栅格 (如 8×8)
Level 1: 分为4个 4×4 子区域
Level 2: 不均匀的子区域继续分为4个 2×2
Level 3: 不均匀的2×2继续分为4个 1×1 (单像元，叶节点)
最大深度 = log2(max(行数, 列数))
对于 1024×1024 栅格，最大深度 = 10
四叉树的优势在于：(1) 自适应压缩 — 均匀区域用大块表示，碎片区域用小块表示；(2) 支持高效的空间查询 — 可以快速判断某点所在区域的值(O(log n)时间)；(3) 天然支持多尺度表达 — 不同层级对应不同分辨率。缺点是要求栅格尺寸为2的幂次(不足时需填充)。
各编码方式对比
编码方式压缩原理典型压缩率查询效率适用数据类型
直接编码无压缩1:1O(1)所有类型
RLE行方向连续值合并2:1 ~ 10:1O(n)分类数据、专题图
链码边界方向编码视边界复杂度而定O(n)边界数据
块码二维矩形块合并优于RLEO(n)大面积均匀区域
四叉树递归四分层次化3:1 ~ 20:1O(log n)分类数据、专题图
五、地图代数 (Map Algebra)
地图代数的概念由C. Dana Tomlin于1983年在其博士论文中首次系统提出，并在1990年出版的《Geographic Information Systems and Cartographic Modeling》一书中全面阐述。地图代数将栅格图层视为二维数学函数，定义了一套完整的基于像元的运算体系。根据运算涉及的空间范围，地图代数将运算分为四类：局部运算(Local)、焦点运算(Focal)、分区运算(Zonal)、全局运算(Global)。
参考文献: Tomlin, C. Dana (1990). Geographic Information Systems and Cartographic Modeling. Prentice-Hall. 这是GIS栅格分析的奠基性著作。
5.1 局部运算(Local Operations)
局部运算对每个像元独立地进行计算，输出像元的值仅取决于对应输入像元的值。这是最基本的地图代数运算，计算效率高(可完全并行化)。典型应用包括波段运算(Band Math)和各种遥感指数计算。
植被指数公式
NDVI (Normalized Difference Vegetation Index) — 归一化差值植被指数，是应用最广泛的植被指数，由Rouse等人于1974年提出：
NDVI = (NIR - Red) / (NIR + Red)
取值范围: [-1, +1]
典型值: 水体 < 0, 裸土 0.1-0.2, 稀疏植被 0.2-0.4, 茂密植被 0.6-0.9
其中: NIR = 近红外波段反射率, Red = 红光波段反射率
无人机获取的高分辨率NDVI影像，展示了植被健康状况的空间差异 (来源: Wikimedia Commons)
灌溉农田的NDVI影像，绿色表示健康植被，红/黄色表示植被稀疏或裸土 (来源: Wikimedia Commons)
EVI (Enhanced Vegetation Index) — 增强型植被指数，由Huete等人于2002年为MODIS传感器设计，修正了NDVI在高生物量地区饱和以及大气和土壤背景干扰的问题：
EVI = G × (NIR - Red) / (NIR + C1 × Red - C2 × Blue + L)
MODIS标准系数: G = 2.5, C1 = 6, C2 = 7.5, L = 1
其中: G = 增益系数, C1/C2 = 大气校正系数, L = 土壤背景调节系数
Blue = 蓝光波段反射率
SAVI (Soil Adjusted Vegetation Index) — 土壤调节植被指数，由Huete于1988年提出，引入了土壤亮度校正因子L，减少了土壤背景对植被指数的影响：
SAVI = (1 + L) × (NIR - Red) / (NIR + Red + L)
L = 0.5 (适用于中等植被覆盖度的经验值)
当 L = 0 时，SAVI退化为NDVI
5.2 焦点运算(Focal Operations)
焦点运算(也称邻域运算 Neighborhood Operations)以每个像元为中心，在其周围窗口(Kernel/Window)范围内进行计算。输出像元的值取决于该像元及其邻域内所有像元的值。窗口通常为3×3、5×5、7×7等奇数尺寸的方形窗口，也可以是圆形或自定义形状。
卷积滤波(Convolution Filtering)
卷积是焦点运算最核心的实现方式。一个卷积核(Kernel)在栅格上滑动，对每个位置，输出值等于窗口内像元值与对应核权重的加权和。
常用卷积核示例：
均值滤波(Mean Filter) 3×3: 拉普拉斯锐化(Laplacian) 3×3:
| 1/9 1/9 1/9 | | 0 -1 0 |
| 1/9 1/9 1/9 | | -1 4 -1 |
| 1/9 1/9 1/9 | | 0 -1 0 |
高斯滤波(Gaussian) 3×3: Sobel算子(水平梯度):
| 1/16 2/16 1/16 | | -1 0 1 |
| 2/16 4/16 2/16 | | -2 0 2 |
| 1/16 2/16 1/16 | | -1 0 1 |
坡度(Slope)计算
坡度是DEM最基本的衍生产品之一。最常用的算法是Horn (1981)三阶有限差分法，使用3×3邻域窗口：
设3×3邻域DEM值为:
| a b c |
| d e f |
| g h i |
东西方向梯度: dz/dx = [(c + 2f + i) - (a + 2d + g)] / (8 × cellsize)
南北方向梯度: dz/dy = [(g + 2h + i) - (a + 2b + c)] / (8 × cellsize)
坡度(度): Slope = arctan( sqrt( (dz/dx)² + (dz/dy)² ) ) × 180 / π
坡度(百分比): Slope% = sqrt( (dz/dx)² + (dz/dy)² ) × 100
坡向(Aspect)计算
Aspect = arctan2(dz/dy, -dz/dx) × 180 / π
结果转换为 [0°, 360°) 范围：
0°/360° = 北坡, 90° = 东坡, 180° = 南坡, 270° = 西坡
平坦区域(坡度=0)的坡向设为 -1 (无方向)
5.3 分区运算(Zonal Operations)
分区运算以一个分区图层(Zone Layer)定义区域边界，对另一个值图层(Value Layer)在每个分区内进行汇总统计。输出的每个像元值等于其所在分区内所有像元的统计值。
典型应用包括：
按行政区统计平均高程：分区图层=县级行政区(编码为整型栅格)，值图层=DEM，运算=求均值
按流域统计年降水量：分区图层=流域边界栅格，值图层=年降水栅格，运算=求和
按土地利用类型统计NDVI：分区图层=土地利用分类图，值图层=NDVI栅格，运算=求均值/标准差
常用的分区统计函数包括：ZonalMean(分区均值)、ZonalSum(分区求和)、ZonalMax/Min(分区最大/最小值)、ZonalSTD(分区标准差)、ZonalMajority(分区众数)、ZonalRange(分区极差)等。
5.4 全局运算(Global Operations)
全局运算的输出像元值可能取决于整个栅格的所有像元。这类运算通常涉及最优路径搜索、距离累积或累积效应分析，计算复杂度远高于前三类。
视域分析(Viewshed Analysis)
从DEM上的一个或多个观察点出发，沿所有方向进行视线追踪(Line of Sight)，判断每个像元是否可被观察到。视域分析的结果是一个二值栅格(0=不可见, 1=可见)。广泛应用于通信基站选址、景观视觉影响评价、军事阵地选择等。
水流方向分析(Flow Direction)
D8算法(Deterministic 8-node)是最经典的水流方向算法：对每个像元，比较其与8个邻域像元的高程差(或坡度)，水流流向最陡降坡方向。输出为方向编码(1=东, 2=东南, 4=南, 8=西南, 16=西, 32=西北, 64=北, 128=东北)。D8的缺点是水流只能流向一个方向(无法模拟分散流)，改进算法如D-infinity (Tarboton, 1997)允许流向任意角度。
成本距离分析(Cost Distance)
成本距离分析基于一个成本表面(Cost Surface)计算从源点(Source)到每个像元的最小累积成本。与简单的欧氏距离不同，成本距离考虑了穿越每个像元的"代价"(如地形坡度、土地利用类型、道路通行性)。应用包括应急救援路径规划、野生动物廊道设计、输电线路选址等。其核心算法基于Dijkstra最短路径在栅格空间上的推广。
六、重采样方法 (Resampling)
当栅格数据需要进行投影变换、几何校正、分辨率转换或配准对齐时，需要将原始像元值映射到新的网格位置上。由于新网格点通常不会精确落在原始像元中心，需要通过重采样(Resampling)方法来估算新位置的像元值。三种最常用的重采样方法各有特点：
6.1 最近邻法(Nearest Neighbor)
将输出像元的值设为距离最近的输入像元值。这是最简单、最快速的重采样方法。
最近邻法特点
优点：计算速度最快(O(1)每像元)；不改变原始像元值，保持数据的原始分类含义；适用于分类数据(土地利用类型、地质岩性等)和专题栅格。
缺点：输出影像可能出现锯齿状边缘和几何位移(最大偏移为半个像元)；不适用于连续数据，会产生明显的阶梯效应。
6.2 双线性内插(Bilinear Interpolation)
利用输出像元位置周围最近的4个(2×2)输入像元值，通过两次线性插值计算加权平均值。
双线性内插公式：
设目标点坐标为 (x, y)，周围4个像元值为 Q11, Q12, Q21, Q22
先在x方向插值：
R1 = Q11 × (x2 - x)/(x2 - x1) + Q21 × (x - x1)/(x2 - x1)
R2 = Q12 × (x2 - x)/(x2 - x1) + Q22 × (x - x1)/(x2 - x1)
再在y方向插值：
P = R1 × (y2 - y)/(y2 - y1) + R2 × (y - y1)/(y2 - y1)
双线性内插特点
优点：输出影像平滑、无锯齿；计算量适中；适用于连续数据(高程、温度、反射率)。
缺点：产生轻微模糊(低通滤波效应)；会改变原始值(不适合分类数据)；边缘细节有一定损失。
6.3 三次卷积(Cubic Convolution)
利用输出像元位置周围最近的16个(4×4)输入像元值，通过三次多项式插值计算。这是精度最高但计算量最大的方法。
三次卷积核函数 h(x):
h(x) = (a+2)|x|³ - (a+3)|x|² + 1, 当 |x| ≤ 1
h(x) = a|x|³ - 5a|x|² + 8a|x| - 4a, 当 1 < |x| ≤ 2
h(x) = 0, 当 |x| > 2
其中 a = -0.5 (Catmull-Rom样条) 或 a = -1 (更锐利)
三次卷积特点
优点：输出影像最清晰、细节保持最好；兼具平滑性和锐度；是遥感影像正射校正的首选方法。
缺点：计算量最大(需处理16个输入像元)；可能产生轻微的振铃效应(Ringing Artifact)；输出值可能略超出输入值范围。
重采样方法对比
方法采样窗口计算复杂度平滑度适用数据类型是否改变原值
最近邻1×1最低无(锯齿)分类/专题数据否
双线性内插2×2中等适中连续数据是
三次卷积4×4最高高(最清晰)连续数据(遥感影像)是
重要提示
对于分类数据(如土地利用类型图)，必须使用最近邻法。双线性内插和三次卷积会产生原始分类中不存在的中间值(例如在"森林=3"和"水体=5"之间插出"4"，但4可能代表"草地"，这种结果毫无意义)。对于连续型遥感影像，推荐使用三次卷积以获得最佳视觉质量。
七、主要栅格数据格式
栅格数据格式定义了像元值、地理参考信息、元数据等如何组织和存储在文件中。不同的应用领域和软件平台发展出了多种栅格格式，了解这些格式的特点对于数据交换和选择合适的存储方案至关重要。
Landsat 8 OLI传感器获取的迈阿密真彩色合成影像 — GeoTIFF格式是此类卫星数据的标准分发格式 (来源: NASA, Wikimedia Commons)
7.1 GeoTIFF
GeoTIFF是目前最通用的地理栅格数据格式，由TIFF(Tagged Image File Format)扩展地理参考标签而来。它在标准TIFF文件中嵌入了坐标系统(Coordinate Reference System)、仿射变换参数、投影信息等地理元数据，使得栅格数据具有空间自描述能力。
核心优势：行业标准，几乎所有GIS/遥感软件都支持；支持多波段、多种数据类型(8/16/32位整型和浮点型)；支持内部分块(Tiling)和金字塔(Overview)以提高大影像访问效率
压缩选项：LZW(无损)、Deflate(无损)、JPEG(有损)、ZSTD(高效无损)
Cloud Optimized GeoTIFF (COG)：一种特殊的GeoTIFF组织方式，通过HTTP Range请求支持直接从云端按需读取栅格数据的部分区域，无需下载整个文件。是云原生GIS(Cloud-Native GIS)时代的标准格式
文件扩展名：.tif 或 .tiff
7.2 NetCDF (Network Common Data Form)
NetCDF是气候科学和海洋学领域的主流数据格式，由Unidata开发维护。它是一种自描述的、面向数组的多维数据格式。
核心特点：原生支持多维数组(经度×纬度×时间×高度等)；内置CF(Climate and Forecast)元数据约定，描述变量单位、坐标、网格映射等
典型应用：全球气候模式输出(CMIP6)、ERA5再分析数据、海洋遥感数据(SST、海面高度)
版本：NetCDF-3(经典格式，2GB限制)、NetCDF-4(基于HDF5，支持压缩和更大文件)
文件扩展名：.nc 或 .nc4
7.3 HDF (Hierarchical Data Format)
HDF是NASA采用的标准遥感数据分发格式，由NCSA(美国国家超级计算应用中心)开发。HDF采用分层目录结构，可在一个文件中存储多种类型的数据(栅格、表格、元数据)。
HDF4：MODIS、MISR等NASA经典遥感产品的标准格式(文件扩展名 .hdf)
HDF5：更现代的版本，支持更大文件和更灵活的数据组织(文件扩展名 .h5 或 .he5)。Sentinel-5P、ICESat-2、GPM降水产品使用此格式
HDF-EOS：HDF的扩展，增加了Swath(卫星条带)、Grid(网格)、Point(点)三种地球科学数据结构
7.4 ERDAS IMG
IMG格式是ERDAS IMAGINE软件的原生栅格格式。它采用分层分块的存储结构(Hierarchical File Architecture, HFA)，支持金字塔层、多波段、属性表和统计信息。虽然由商业软件开发，但GDAL库提供了完整的读写支持。文件扩展名为 .img。
7.5 ASCII Grid (Esri Grid)
ASCII Grid是一种纯文本的栅格格式，人类可直接阅读和编辑。文件以6行头信息开始(ncols, nrows, xllcorner, yllcorner, cellsize, NODATA_value)，随后是以空格分隔的像元值矩阵。
ASCII Grid 文件示例:
ncols 5
nrows 4
xllcorner 100.0
yllcorner 30.0
cellsize 0.01
NODATA_value -9999
120 135 142 138 -9999
118 130 145 140 136
115 125 -9999 137 133
110 120 128 132 130
ASCII Grid的优点是格式简单透明、便于调试和教学；缺点是文件体积大(文本编码比二进制大3-5倍)、不支持多波段和投影信息。文件扩展名为 .asc。
栅格格式对比总览
格式类型多波段压缩多维数组主要应用领域GDAL支持
GeoTIFF二进制是LZW/Deflate/JPEG否通用GIS/遥感完整读写
COG二进制是同GeoTIFF否云原生GIS完整读写
NetCDF二进制是Deflate/Shuffle是(核心优势)气候/海洋完整读写
HDF4/5二进制是多种是NASA遥感产品完整读写
IMG二进制是内置否ERDAS生态系统完整读写
ASCII Grid文本否无否教学/数据交换完整读写
八、栅格与矢量的转换
栅格数据和矢量数据各有优势，在实际GIS项目中经常需要在两者之间转换。理解转换的原理和注意事项，对于保证数据质量和分析结果的准确性至关重要。
8.1 栅格转矢量(Raster to Vector)
栅格转矢量是将连续的像元网格转化为离散的点、线或面要素的过程。根据目标矢量类型，转换方法不同：
栅格转面(Polygonize)
将具有相同像元值的连通区域(Connected Components)转换为多边形。算法首先通过连通域标记(Connected Component Labeling)识别同值区域，然后追踪每个区域的边界生成多边形。GDAL的 gdal_polygonize 工具和ArcGIS的 Raster to Polygon 工具都实现了此功能。
转换注意事项：
阶梯效应(Staircase Effect)：转换后的多边形边界呈锯齿状(因为像元本身是方形的)。可通过后续的平滑(Smooth)操作改善
碎片问题：噪声像元会生成大量微小多边形，转换前通常需要进行众数滤波(Majority Filter)去除椒盐噪声
数据膨胀：高分辨率栅格转矢量可能产生海量多边形顶点，导致文件极大
栅格转线 — 等值线追踪(Contour Tracing)
从连续数值栅格(如DEM)中提取等值线(Contour Lines)是一种经典的栅格转矢量操作。最常用的算法是Marching Squares(二维版本的Marching Cubes)：
遍历所有2×2像元组，根据4个角点值与等值线阈值的大小关系，形成16种可能的构型(Case)
对于每种构型，根据预定义的查找表确定等值线穿过该2×2块的线段位置
通过线性插值确定等值线与像元边界的精确交点
将相邻2×2块中的线段连接成完整的等值线
GDAL的 gdal_contour 工具可以从DEM生成等高线矢量。在GIS中，等高线是DEM的经典可视化方式，也是传统纸质地形图的核心要素。
栅格转点
将每个有效像元(非NoData)转换为一个点要素，点位于像元中心，属性为像元值。这种转换数据量最大(每个像元一个点)，通常仅用于采样或结合插值分析。
8.2 矢量转栅格(Rasterization / Vector to Raster)
矢量转栅格(栅格化)是将矢量要素"烧录(Burn)"到栅格网格中的过程。核心问题是确定每个像元是否被矢量要素覆盖、以及赋予什么值。
栅格化规则
面要素栅格化：判断每个像元中心是否落在多边形内部(Point-in-Polygon测试)。若是，则该像元值设为多边形的指定属性值。也可以使用面积占比法(像元被多边形覆盖的面积比例)来处理边界像元
线要素栅格化：判断线段是否穿过像元。通常使用Bresenham直线算法来确定线段经过的所有像元，将这些像元的值设为线要素的属性值
点要素栅格化：将点所在的像元值设为点要素的属性值。当多个点落在同一像元时，可选择最大值/最小值/均值/计数等聚合方式
栅格化的关键参数
参数说明影响
输出分辨率(Cell Size)栅格像元的地面尺寸分辨率过大导致细节丢失，过小导致文件巨大且可能引入虚假精度
烧录值(Burn Value)赋予被覆盖像元的值可以是固定值(如1)或来自矢量属性字段
空间范围(Extent)输出栅格的地理范围需覆盖所有矢量要素或与其他栅格对齐
NoData值未被矢量覆盖的像元值通常设为-9999或0
ALL_TOUCHED选项是否栅格化所有被触及的像元TRUE=要素触及的所有像元都赋值; FALSE=仅中心点被覆盖的像元
转换精度与信息损失
无论是栅格转矢量还是矢量转栅格，都不可避免地存在信息损失。栅格化时，小于半个像元的矢量细节会丢失(类似于信号采样中的奈奎斯特定理 — 空间分辨率应至少为最小特征尺寸的一半)。矢量化时，原本连续的曲线边界变为阶梯状折线。因此，应尽量避免反复在两种格式之间转换，以免累积误差。在实际工程中，应根据分析目的选择最合适的数据模型，减少不必要的格式转换。
8.3 栅格与矢量的协同分析
在现代GIS分析中，栅格和矢量数据往往协同使用而非相互替代。例如：
分区统计(Zonal Statistics)：用矢量多边形(如行政区边界)作为分区，统计栅格数据(如NDVI)在每个分区内的均值、总和等。无需将矢量转栅格，GIS软件直接支持矢量-栅格叠加分析
裁剪提取(Clip/Extract by Mask)：用矢量边界裁剪栅格数据，仅保留边界内的像元值，边界外设为NoData
采样(Sample/Extract Values)：在矢量点位置提取栅格值(如在气象站点位置提取DEM高程)，常用于模型验证和地统计分析
总结
栅格数据模型是GIS空间信息表达的两大基石之一。从最基本的像元结构到复杂的地图代数运算，从四种分辨率的权衡到多种编码压缩算法，从丰富的数据格式生态到栅格-矢量转换的原理与实践，栅格数据贯穿了遥感影像处理、地形分析、环境建模、城市规划等GIS应用的方方面面。
理解栅格数据的核心要点包括：(1) 像元是栅格的基本单位，NoData是不可忽视的重要概念；(2) 空间、光谱、时间、辐射四种分辨率共同决定了数据的信息承载能力；(3) DEM/DSM/DTM虽然都是高程数据但含义不同，选择错误的高程模型会导致分析结果出现系统性偏差；(4) 四叉树等压缩编码在大数据时代仍然具有重要的实用价值；(5) Tomlin的地图代数四类运算体系(局部/焦点/分区/全局)为栅格空间分析提供了完整的理论框架；(6) 重采样方法的选择必须匹配数据类型(分类vs连续)；(7) GeoTIFF/COG是当前最通用的栅格格式，而NetCDF和HDF在专业领域不可替代；(8) 栅格与矢量的转换存在固有的信息损失，协同分析优于频繁转换。
随着云计算、AI和星载传感器技术的发展，栅格数据的规模正以指数级增长。Google Earth Engine、Microsoft Planetary Computer等云平台使得PB级栅格数据的分析成为可能。深度学习技术(如U-Net、ResNet)在遥感影像的语义分割、变化检测等方面展现出巨大潜力。理解栅格数据的底层原理，是充分利用这些新兴技术的基础。
主要参考文献：
1. Tomlin, C.D. (1990). Geographic Information Systems and Cartographic Modeling. Prentice-Hall.
2. Longley, P.A., Goodchild, M.F., Maguire, D.J., Rhind, D.W. (2015). Geographic Information Science and Systems. 4th ed. Wiley.
3. Jensen, J.R. (2015). Introductory Digital Image Processing: A Remote Sensing Perspective. 4th ed. Pearson.
4. Rouse, J.W. et al. (1974). Monitoring vegetation systems in the Great Plains with ERTS. NASA SP-351.
5. Huete, A.R. et al. (2002). Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sensing of Environment, 83, 195-213.
6. Horn, B.K.P. (1981). Hill shading and the reflectance map. Proceedings of the IEEE, 69(1), 14-47.