高分二号影像预处理全流程：从L1A到正射融合

简介：面向使用高分二号卫星影像的遥感从业者、GIS学习者及科研人员，这份PDF以问答形式系统梳理了数据版本与分辨率、WGS84坐标系、原始数据挑选标准、处理成果格式及适用软件等六类高频问题。资源共1个文件，为PDF格式，压缩包大小约183KB，内容精炼，适合随用随取。目前已有193人学习浏览。文中明确指出全色影像星下点分辨率为0.8米、多光谱为3.2米，并给出云量小于5%、分散云量不超过15%等筛选阈值；同时解释了xml、rpb、TIFF、JPG四种数据文件的用途，以及ENVI中默认显示为.dat、ArcGIS中为tiff等格式差异。此外，资源还说明数据经融合处理后分辨率提升至0.8米，并列举ArcGIS、ENVI、IDL、易康、Erdas IMAGE等适用软件，能帮助读者在预处理、融合和正射校正时少走弯路，快速定位问题。

1. 高分二号影像的问题，多数出在预处理而不是卫星本身

高分二号（GF-2）作为国产亚米级光学卫星，全色分辨率0.8米、多光谱3.2米，单景幅宽45公里，在国土测绘、农业普查、城市精细化管理里使用率极高。但很多从业者第一次拿到GF-2数据时，会被它的“原始感”吓到：多光谱和全色是两个独立文件，坐标系是WGS84但投影五花八门，影像边缘有明显黑边，调色后色彩发灰，甚至与矢量边界对不齐。这些现象并不是卫星出了问题，而是预处理链路没走完整。这篇文章按实际生产流程，从数据特性、GDAL命令、参数取舍到高频故障排查，讲清高分二号卫星影像从原始L1A级到可应用正射成果的完整路径，适合遥感数据处理工程师、GIS开发者和刚接触国产卫星影像的学生参考。

2. 高分二号数据特性与常见问题根源：从波谱、几何到存储格式

2.1 全色与多光谱分开存储，融合前必须搞清的三件事

GF-2的L1A级产品通常包含一个全色（PAN）文件和四个多光谱（MSS）波段（蓝、绿、红、近红外）。如果直接打开影像，全色是单波段灰度，多光谱则是四波段彩色合成，两者分辨率不同，空间范围也略有差异。常见做法是先检查RPC文件是否齐全，再确认各波段的中心波长和定标参数。RPC（有理多项式系数）是正射校正和区域网平差的几何基础，缺少RPC文件时，影像无法做精确的地理定位。另外，全色与多光谱的存储顺序、位深（通常为16位）和无效值（通常为0）都影响后续处理，建议先用GDAL读取元数据，确认波段数与数据类型。

以我处理过的GF-2数据为例，供应商交付的压缩包中有时会混入不同时间、不同条带的RPC文件，直接解压后gdalwarp可能提示“No RPC”或定位偏差数十米。正确的第一步是把每个影像与同名/同ID的RPC文件放在同一目录，并确保文件名字段（如“_RPC.txt”）没有缺失。整景数据中如果存在多个传感器（PMS1和PMS2），它们的RPC参数并不通用，必须对应各自的影像文件。这一步看似基础，却直接决定后续正射精度。

2.2 坐标系统与投影：为什么影像在GIS里老是对不上

高分二号的L1A级产品默认采用WGS84椭球，投影常用UTM或高斯-克吕格。很多用户在GIS里直接叠加矢量，发现影像与路口、建筑边界有几十米甚至上百米的偏移，这是因为L1A级影像本身只做了传感器校正，没有做地形校正。正射校正需要结合DEM（如SRTM、AW3D30）和RPC模型，把影像投影到目标坐标系。如果目标区域在高纬度或跨带，UTM带号选择错误会导致影像严重变形；如果数据提供方给的投影信息不完整，需要用gdalinfo检查SRS字符串。

另一个容易忽略的点是椭球与基准面。虽然GDAL默认WGS84，但地方坐标系（如CGCS2000）与WGS84在大部分区域差异极小（厘米到米级），如果项目要求CGCS2000，可以直接使用EPSG:4490或EPSG:4526系列。但在没有七参数的情况下，直接把WGS84数据转成CGCS2000，可能引入1-2米偏差，这在0.8米分辨率影像上不可接受。常见做法是先用gdalwarp将影像投影到目标UTM带，再通过地面控制点做局部配准，而不是依赖参数转换。

2.3 GF-2核心参数与典型问题对照表

下表列出GF-2关键参数与直接关联的典型问题，方便排查时对齐。

这张表在项目沟通中很有用。比如用户反馈“影像太暗”，根因不是数据问题，而是16位DN值范围在0-2047之间，默认拉伸会把低端压缩到黑。再比如“图像发绿”，多光谱波段顺序被误排为NIR、R、G，绿色物体被映射成洋红色。对照表格能快速缩小排查范围。

2.4 用gdalinfo快速定位元数据问题

拿到一景GF-2影像后，我习惯先用两行命令确认基础信息：

gdalinfo GF2_PMS1_E114.5_N30.7_20230101_L1A_MSS.tif gdalinfo GF2_PMS1_E114.5_N30.7_20230101_L1A_PAN.tif

重点关注输出中的Size、Band、NoData Value、Coordinate System、以及末尾的RPC metadata。如果MSS和PAN的RPC数值不一致，说明两份RPC来自不同处理批次，融合前必须做同名点配准。如果Coordinate System显示为未知，说明投影信息丢失，需要向数据提供方索取正确的坐标定义，或根据影像角点坐标自行推断UTM带号。查看Band信息时，注意多光谱影像是否为四个波段（BGRN），有些产品会附带一个全色合成文件，不要把它当成多光谱直接使用。

以上元数据检查通常在几分钟内完成，却可以避免后续所有步骤在错误基础上重复返工。这也是为什么我把这一步放在预处理之前：先确认输入物的状态，再决定采用哪种校正策略。

3. 用GDAL与Python跑通高分二号影像预处理的最小闭环

3.1 安装环境与读取影像元数据

处理高分二号影像的软件栈很多，常见的是ENVI+ERDAS，但工程化落地我更推荐GDAL与rasterio，因为它们可以嵌入Python脚本，批量处理多景影像。环境安装用conda比较直接：

conda create -n gis python=3.11 conda activate gis conda install -c conda-forge gdal rasterio numpy scikit-image

我一般使用rasterio读取影像，而不是直接调GDAL，因为rasterio的API更友好，且和NumPy无缝对接。以下代码读取MSS影像并打印关键信息：

import rasterio with rasterio.open("GF2_PMS1_E114.5_N30.7_20230101_L1A_MSS.tif") as src: print(src.meta) # 包含形状、波段数、数据类型、变换矩阵 print(src.crs) # 坐标系统 print(src.transform) # 仿射变换参数 print(src.band_names) # 波段名称（如有） mss_array = src.read() # 返回 (4, height, width) 的数组

src.meta输出的dtype是uint16，说明影像位深为16bit；transform给出左上角坐标和像素尺寸；crs一般是EPSG:32650或4326。如果crs为None，说明投影未定义，需要手动指定。读取到的mss_array是原始DN值，后续定标和拉伸都在这个数组上进行。注意直接调用src.read()会把所有波段加载进内存，对于1.2万×1.2万像素的影像约占用4×12000×12000×2字节，约1.1GB，运行时请确保内存充足。

3.2 辐射定标与大气校正的实用参数

辐射定标的目标是把DN值转换为大气顶部反射率。GF-2头文件或说明文档中会给出每个波段的绝对辐射定标系数（gain和offset），常见形式为L = DN * gain + offset。计算TOA反射率时还需要太阳天顶角和影像采集时的日地距离。以下代码实现从定标后辐亮度到反射率的转换：

import numpy as np def toa_reflectance(dn, gain, offset, solar_zenith, esun): radiance = dn * gain + offset cos_theta = np.cos(np.radians(solar_zenith)) reflectance = np.pi * radiance * 1.0**2 / (esun * cos_theta) return reflectance

参数说明：esun是各波段大气顶部的太阳辐照度，GF-2数据说明书中会给出，也可以在NASA的太阳光谱常数表中近似获取。solar_zenith在影像XML头文件中有记录，通常以角度为单位。日地距离1.0对应平均日地距离，实际应用中应根据成像日期修正（约±3%误差）。如果只做相对比较，可以跳过大气校正，直接使用TOA反射率；但涉及植被指数阈值或跨时相绝对对比，建议使用6S或FLAASH。GDAL本身不提供大气校正，常见做法是调用py6s库或导到ENVI里完成，这里不再展开。

3.3 正射校正与影像融合命令

正射校正是把原始L1A影像纠正到地形起伏对应的地理位置。常见做法是用gdalwarp配合RPC和DEM，命令行如下：

gdalwarp -rpc -t_srs EPSG:32650 -tr 0.8 0.8 -r cubic \ -srcnodata 0 -dstnodata 0 \ GF2_PMS1_E114.5_N30.7_20230101_L1A_MSS.tif \ GF2_mss_ortho.tif

-rpc启用RPC模型，-t_srs指定目标坐标系，这里为UTM 50N；-tr设置输出像元大小为0.8米，与全色分辨率一致，这样后续融合不需要再重采样多光谱。-r cubic使用三次卷积插值，适合地物纹理丰富的区域；如果对速度敏感，-r bilinear也可接受。-srcnodata 0和-dstnodata 0把像元值为0的地区视为无效，避免黑边被插值成灰色。注意DEM文件需要放在当前目录或通过-geoloc方式指定，GDAL会自动查找。如果RPC文件缺失，可以改用-gcp方式，但精度受控制点分布影响较大，无法在全景范围保持稳定。

多光谱与全色融合常用的方法是Gram-Schmidt谱锐化。ENVI的Gram-Schmidt工具提供了较稳定的实现，但工程化批量处理可以用Orfeo ToolBox的otbcli_Fusion命令：

otbcli_Fusion -method.simple.nbBands 4 \ -method.simple.algorithm pca \ -il GF2_mss_ortho.tif GF2_pan_ortho.tif \ -out GF2_fused.tif uint16

pca算法将多光谱做PCA变换，再把全色替换第一主成分后逆变换，速度比GS快，色彩失真略大。如果项目对光谱保真要求高，建议用gs算法。融合前务必确保两个影像已配准到同一个网格，否则融合结果会出现重影或光谱扭曲。一个快速验证方法是在QGIS中把全色影像叠加到多光谱上，调成半透明，检查道路和屋顶边缘是否有叠影。

3.4 预处理步骤与关键参数表

下面列出完整预处理链路中每一步的输入输出和关键参数，便于脚本编排：

按这个表操作，能覆盖大部分生产场景。需要注意的是，正射校正应在辐射定标之后进行，因为几何校正不会改变DN值，但重采样过程会损失辐射精度，先定标再重采样更合理。有些团队先融合后正射，虽然流程上可行，但RPC模型在融合后的影像上不再适用，就需要重新生成RPC或依赖GCP，反而增加复杂度。

4. 高分二号高频问题排查手册：黑边、云影、色彩断层与定位漂移

4.1 黑边与无效值：快速识别并掩膜

高分二号原始影像四周的黑色无效区是几乎每个使用者都会遇到的问题。如果不处理，-scale计算统计值时会把大量0值计入，导致有效地物被压缩到非常窄的亮度区间，影像看起来又黑又平。有效的处理思路是先构建二值掩膜，再使用gdal_fillnodata填充或裁剪。以下Python代码生成有效区域掩膜：

from osgeo import gdal import numpy as np ds = gdal.Open("GF2_mss_ortho.tif") band = ds.GetRasterBand(1) arr = band.ReadAsArray() mask = (arr > 10) & (arr < 65535)

这里阈值下限设为10，是为了排除极低值噪声；上限65535排除可能的饱和像元。但要注意，云影和深水体也可能被误判为无效区。如果影像中有大面积深水，建议结合其他波段判断：例如近红外波段中水体反射率低，但蓝波段通常有较高值，可以增加蓝波段的约束条件。掩膜生成后，使用gdalwarp的-cutline或-dstnodata输出无黑边影像：

gdalwarp -dstnodata 0 -cutline mask.shp input.tif output.tif

其中mask.shp可从掩膜栅格矢量化得到。另一种做法是使用gdal_fillnodata -md 50填充黑边，但这种方法会制造人造边界线，只建议用于后续拼接的过渡区域。

4.2 云检测与薄云去除的两种落地做法

云污染直接影响影像可用性，尤其在多雨地区。落地做法有两种：基于规则检测和基于深度学习分割。基于规则的检测在GF-2上很实用，因为云在蓝波段和绿波段亮度高，而近红外波段亮度相对低；云影则表现为暗且近红外低。一个简单示例：

# 假设b1蓝，b2绿，b3红，b4近红外 cloud_mask = (b1 > 1200) & (b4 < 2500) shadow_mask = (b3 < 300) & (b4 < 300)

这里阈值是针对16位DN值的经验值，不同季节和太阳高度角下会变化，需要调整。更稳定的做法是先做辐射定标，把DN值转换为反射率，再用固定阈值（如云反射率>0.35，近红外<0.2）。基于深度学习的方法更抗噪，常见的是使用S2cloudless模型，但需要把GF-2重采样到10米并与Sentinel-2特征对齐，迁移成本较高。对于单景场景，我会首选规则检测，并在结果上用形态学开运算去除孤立噪点。

薄云去除比检测复杂。一个常用方法是利用同一区域的另一期清晰影像，对云区做直方图匹配，但效果因云厚度而异。薄云区域用暗目标减法可能有效：选取云区邻近的无云暗目标，估算大气路径辐射并从云区减去。我在实际项目中不会依赖薄云去除，而是直接生成云掩膜，把云区标记为“数据缺失”，在后续时序分析中跳过。

4.3 色彩断层与直方图匹配

色彩断层通常表现为红蓝通道互换、绿色植被变成洋红，或影像呈灰蓝色。根因多半是波段顺序错误。GF-2的波段顺序是蓝（450-520nm）、绿（520-590nm）、红（630-690nm）、近红外（770-890nm）。在显示RGB时，标准组合是3,2,1波段；植被分析常用4,3,2组合。如果直接使用1,2,3组合，蓝通道被赋予红波段，红通道被赋予蓝波段，色彩完全错乱。解决方法是用gdal_translate重排波段：

gdal_translate -ot Byte -scale 0 5000 0 255 -b 3 -b 2 -b 1 input_fused.tif rgb.tif

-b 3 -b 2 -b 1从输入影像中选择红、绿、蓝三个波段，-scale设置0至5000的DN映射到0至255，这样既纠正了顺序，又完成了8位拉伸。如果影像仍然偏暗，需要把-scale的下限设为直方图第2百分位，上限设为第98百分位。可以用gdalinfo -hist或NumPy计算百分位数，再回填到命令中。

多景拼接时，各景影像因大气条件不同会形成明显晒斑。常见做法是先选择一景作为参考，在重叠区建立直方图映射。Python的scikit-image库提供了现成的直方图匹配函数：

from skimage.exposure import match_histograms import numpy as np # target和reference均为(h,w,3)的uint8数组 matched = match_histograms(target, reference, channel_axis=-1)

注意这里要求输入影像已经做了地理配准并裁剪到相同范围。匹配后色彩一致性明显改善，但会改变原始辐射值，如果后续要做定量反演，必须在匹配前保留未修改版本。

4.4 定位漂移：检查RPC文件和地面控制点

如果正射影像与高精度矢量数据叠加时仍有偏移，第一步检查RPC文件是否对应。GF-2的RPC文件通常以文本形式提供，文件名会包含影像ID，如GF2_PMS1_E114.5_N30.7_20230101_L1A_RPC.txt。如果该文件被重命名或与影像不同批次，gdalwarp读取到的RPC会失效。一个快速检查方式是用gdalinfo输出RPC块，确认八个系数是否与原始文件一致。

其次，DEM的精度直接影响定位。在平坦地区，SRTM-30m和AW3D30的差异不大；但在陡峭山地，30m和90mDEM会造成数米到十几米的水平位移。建议在山区使用AW3D30或TanDEM-X，并开启-rpc-demy参数指定DEM文件。如果仍然偏移，就需要用地面控制点校正。在QGIS的Georeferencer中导入影像，手动选取道路交叉口或明显地物，至少10个控制点并均匀分布，然后输出GCP列表。GDAL允许用-gcp在warp前追加控制点：

gdalwarp -gcp 1200 3400 500000 3400000 -gcp 2300 4500 501000 3401000 \ -t_srs EPSG:32650 input.tif output_gcp.tif

这种方法适用于影像内部存在局部变形的情况，但控制点数量不足20个时，输出边缘依然可能漂移。更可靠的做法是用区域网平差工具（如PCI、Pix4D或GDAL的-rpc联合平差），但这已经超出单景处理范围。

5. 高分二号历史影像对比的进阶技巧：从单景到时序分析

5.1 多期影像的配准与相对辐射归一化

历史卫星影像的价值在于对比分析。高分二号自2014年发射以来积累了多个时期的存档数据，但不同时相的影像，即使是同一传感器，也会因太阳高度角、大气条件、传感器状态差异而存在辐射不一致。直接相减会得到大量伪变化。我一般先做配准，再做相对辐射归一化。配准的要点是让多期影像落在同一网格上，分辨率建议统一为0.8米或1米：

gdalwarp -r cubic -tr 1 1 -t_srs EPSG:32650 source_time1.tif common_grid.tif gdalwarp -r cubic -tr 1 1 -t_srs EPSG:32650 source_time2.tif common_grid2.tif

配准后检查交叉点处的道路是否重合，如果偏移超过一个像元，就需要手动添加GCP。这一步不可跳过，因为0.8米分辨率下，一个像元的位移就会被误判为变化。

相对辐射归一化常用线性回归法，在重叠区域选择不变地物，如裸地、混凝土地面、不变的道路表面。以下代码实现基于最小二乘的归一化：

from numpy.linalg import lstsq x = ref_flat[mask] # 参考影像像素值 y = tar_flat[mask] # 待校正影像像素值 A = np.vstack([x, np.ones_like(x)]).T k, b = lstsq(A, y, rcond=None)[0] tar_norm = k * tar + b

选择样本时，可以用坡度小于5度的区域、NDVI在0.1以下且不变的区域来排除植被和地形阴影。样本数量没有硬性要求，但至少覆盖5000个像元，回归才稳定。归一化后可以计算NDVI等指数，再进行差值分析。

5.2 用波段运算做变化检测的注意点

变化检测时，我会先用NDVI差值区分植被减少与增加，再用R-NIR散点图排除噪声。光谱指数的阈值需要根据影像的时间分布调整，例如秋冬季节植被叶绿素下降，NDVI普遍比夏季低0.1左右，因此不能用一个固定阈值跨季度比较。合理做法是先把影像按季节分组，只比较同月或同季的影像。

另外，高分二号的重访周期虽短，但云覆盖和侧摆角度会造成有效影像稀疏。做长时序分析时，我宁可每季度选一景整体干净的影像，而不是用多景含云的影像拼凑。云掩膜后的空洞区域，在时序曲线上应标记为数据缺失，不要用插值填补，否则会产生虚假趋势。对于变化区域的验证，可以用混淆矩阵。随机选取200个像元，人工目视解释为“变化/未变化”，对比检测结果，计算总体精度和Kappa系数。如果精度低于85%，建议调整阈值或加入纹理特征，例如用灰度共生矩阵的对比度作为第二维输入。这一步把历史影像对比从“看图说话”推向可量化的评估，也更容易在项目中落地可信度指标。

本文还有配套的精品资源，点击获取