简介这份资源围绕GF3卫星影像的完整处理流程展开面向遥感影像处理初学者及需要上手SARscape工具的从业者帮助解决从原始数据到成果输出各环节的参数设置与操作规范问题。包内共1个docx文档压缩包约4.58MB以图文步骤形式组织内容便于对照软件界面逐步操作。文档依次覆盖系统设置、数据导入、多视处理、单通道强度数据滤波等核心环节并延伸至地理编码与辐射定标对极化方式选择、制图分辨率设定、Frost滤波窗口参数等关键配置均有说明还给出入射角、采样间隔等元数据读取与分辨率换算思路。目前已有1577人学习下载适合希望系统掌握GF3影像预处理流程、减少参数试错成本的读者参考。1. GF3卫星影像处理从原始数据到可用成果的完整链路拿到一景GF3卫星影像很多人第一反应是直接拖进软件里开始调色。但如果真这么干大概率会在某个环节卡住——要么是数据打不开要么是坐标对不上要么是做完正射校正后发现边缘区域拉伸得没法看。GF3作为C波段合成孔径雷达卫星它的影像处理和光学影像完全是两套逻辑。光学影像处理的核心是辐射校正和大气校正而GF3影像处理的核心是聚焦、辐射定标、斑点滤波、几何校正和地形校正这一整条链路。这篇文章面向的是已经拿到GF3数据、需要把它变成可分析成果的从业者不管你是做地表形变监测、地物分类还是变化检测这套处理流程都是绕不开的基本功。我会按实际操作顺序把每一步的参数设置、常见报错和踩坑经验讲清楚。2. GF3影像处理的前置准备与数据理解2.1 GF3数据的产品分级与选型逻辑GF3卫星影像按处理程度分为多个级别常见的有原始回波数据、单视复数数据、地距探测数据、地理编码数据等。不同级别的数据对应不同的处理起点选错了级别要么做重复工作要么根本没法往下走。产品级别内容适用场景处理起点原始回波未聚焦的原始信号需要自定义聚焦算法聚焦处理单视复数聚焦后的复数数据干涉测量、极化分析辐射定标滤波地距探测检测后的幅度数据强度分析、分类滤波几何校正地理编码投影到地理坐标直接叠加分析精校正后处理我一般建议如果你的目标是做地物分类或变化检测直接从地距探测数据开始最省事如果要做干涉测量必须从单视复数数据起步。选错级别是新手最容易翻车的地方——拿地理编码数据去做干涉等于把已经投影过的数据再投影一次精度损失不可逆。2.2 处理环境的搭建与依赖检查GF3影像处理常用的工具链包括欧空局的开源工具箱、商业遥感软件以及Python生态中的遥感处理库。不管用哪套工具环境搭建时最容易忽略的是依赖库版本冲突。# 以Python生态为例创建独立环境避免版本冲突 conda create -n gf3_process python3.9 conda activate gf3_process # 安装核心依赖注意版本锁定 pip install numpy1.23.5 pip install scipy1.10.1 pip install gdal3.6.2 pip install rasterio1.3.6 pip install snapista1.0.0 # 用于调用欧空局工具箱的Python接口这段脚本的关键在于版本锁定。GDAL和rasterio的版本必须匹配否则读写GF3的GeoTIFF时会报“not recognized as a supported file format”。snapista是调用欧空局工具箱的Python封装如果你习惯用图形界面操作可以跳过这一步但批量处理时命令行效率高得多。参数说明python3.9是经过验证的稳定版本3.10以上在部分遥感库中会有兼容性问题。numpy锁定1.23.5是因为更高版本与scipy 1.10.1存在ABI冲突。gdal 3.6.2对GF3的元数据解析支持最完整。注意如果你在Windows下操作GDAL的安装建议用conda而不是pippip安装的GDAL经常找不到proj数据库路径。2.3 元数据解析与轨道信息提取GF3影像的元数据文件里藏着处理所需的关键参数成像时间、轨道方向、入射角、极化方式、采样间隔等。这些参数直接决定后续辐射定标和几何校正的系数。import xml.etree.ElementTree as ET def parse_gf3_metadata(xml_path): 解析GF3元数据文件提取关键处理参数 tree ET.parse(xml_path) root tree.getroot() # 命名空间处理GF3元数据通常带命名空间 ns {ns: root.tag.split(})[0].strip({)} if } in root.tag else {} meta {} # 提取成像时间 meta[acquisition_time] root.findtext(.//ns:imagingTime, namespacesns) # 提取轨道方向 meta[orbit_direction] root.findtext(.//ns:orbitDirection, namespacesns) # 提取中心入射角 meta[incidence_angle] float(root.findtext(.//ns:incidenceAngle, namespacesns)) # 提取极化方式 meta[polarization] root.findtext(.//ns:polarization, namespacesns) # 提取距离向和方位向采样间隔 meta[range_spacing] float(root.findtext(.//ns:rangeSpacing, namespacesns)) meta[azimuth_spacing] float(root.findtext(.//ns:azimuthSpacing, namespacesns)) return meta # 调用示例 meta parse_gf3_metadata(GF3_XXX_meta.xml) print(f入射角: {meta[incidence_angle]}度) print(f极化方式: {meta[polarization]})这段代码的逻辑是遍历XML树按标签名提取参数。命名空间处理是关键——GF3元数据文件通常带命名空间前缀直接按标签名查找会返回None。入射角参数用于后续辐射定标中的入射角校正采样间隔用于计算像元大小和几何校正的网格间距。如果解析报错“标签未找到”先用文本编辑器打开XML文件确认实际的标签名和命名空间。不同批次的数据元数据格式可能有细微差异我遇到过同一颗卫星不同年份的数据标签名从“incidenceAngle”变成“incidence_angle”的情况。3. GF3影像的核心处理步骤与参数设置3.1 辐射定标把DN值转成后向散射系数GF3影像的原始像元值是DN值没有物理意义。辐射定标就是把DN值转换成后向散射系数让不同时间、不同轨道的影像可以对比。import numpy as np import rasterio def radiometric_calibration(dn_path, output_path, calib_factor, incidence_angle): GF3辐射定标DN值转后向散射系数 with rasterio.open(dn_path) as src: dn src.read(1).astype(np.float32) profile src.profile.copy() # 避免除零 dn[dn 0] np.nan # 定标公式sigma0 (DN^2 / calib_factor) * sin(incidence_angle) # calib_factor来自元数据中的定标常数 sigma0 (dn ** 2) / calib_factor # 入射角校正 sigma0 sigma0 * np.sin(np.radians(incidence_angle)) # 转dB sigma0_db 10 * np.log10(sigma0) # 写入输出 profile.update(dtyperasterio.float32, count1, nodatanp.nan) with rasterio.open(output_path, w, **profile) as dst: dst.write(sigma0_db, 1) return sigma0_db # 调用示例 calib_factor 1.234e5 # 从元数据中读取 sigma0 radiometric_calibration(gf3_dn.tif, gf3_sigma0.tif, calib_factor, 35.2)定标公式中的calib_factor是定标常数必须从元数据中读取不能自己猜。入射角校正这一步很多人会漏掉——如果不做入射角校正同一景影像中近距端和远距端的后向散射系数会有系统性偏差做分类时边缘区域会大量误分。参数说明calib_factor的量级通常在10^4到10^6之间具体值取决于卫星的发射功率和增益设置。入射角从元数据中提取单位是度。转dB是为了压缩动态范围方便后续滤波和显示。注意如果定标后的影像出现大量NaN检查原始DN值中是否有零值区域。GF3影像的边缘区域经常有零值填充这些区域在定标前就要标记为无效。3.2 斑点滤波在去噪和保边之间找平衡合成孔径雷达影像天生有斑点噪声这是相干成像的物理机制决定的。斑点滤波的目标是抑制噪声的同时保留边缘和纹理信息。常用的滤波器有Lee滤波、Frost滤波、Gamma MAP滤波等。from scipy.ndimage import uniform_filter, median_filter def lee_filter(image, window_size7, noise_varNone): Lee滤波基于局部统计的自适应斑点滤波 if noise_var is None: # 估计噪声方差通常用均匀区域的方差 noise_var np.nanvar(image[100:200, 100:200]) # 局部均值和方差 local_mean uniform_filter(image, window_size) local_sq_mean uniform_filter(image ** 2, window_size) local_var local_sq_mean - local_mean ** 2 # Lee滤波权重 weight local_var / (local_var noise_var) weight np.clip(weight, 0, 1) # 滤波结果 filtered local_mean weight * (image - local_mean) return filtered # 调用示例 filtered lee_filter(sigma0, window_size7)Lee滤波的核心思想是用局部方差和噪声方差的比值来决定平滑程度。在均匀区域局部方差接近噪声方差权重趋近于0.5滤波强度大在边缘区域局部方差远大于噪声方差权重趋近于1滤波强度小边缘得以保留。窗口大小的选择是个经验活。7×7是常用值窗口太小去噪不干净窗口太大边缘模糊。我一般会先用7×7跑一遍看效果如果斑点残留明显再调到9×9。noise_var的估计很关键——估大了边缘模糊估小了去噪不足。建议选影像中一块均匀区域比如平静水面或裸地来计算方差。3.3 几何校正与地形校正让影像和地图对得上GF3影像的几何校正分两步先做系统几何校正用轨道参数和地球模型把影像投影到地理坐标再做精校正用地面控制点或参考影像消除残余误差。地形校正则是消除地形起伏引起的辐射畸变。def geometric_correction(input_path, output_path, dem_path, gcp_listNone): GF3几何校正与地形校正 import subprocess # 使用欧空局工具箱的命令行接口 # 第一步应用轨道文件进行系统几何校正 cmd_orbit [ gpt, Apply-Orbit-File, -Ssource input_path, -PorbitTypeSentinel Precise, -t, temp_orbit.tif ] subprocess.run(cmd_orbit, checkTrue) # 第二步地形校正需要外部DEM cmd_terrain [ gpt, Terrain-Correction, -Ssourcetemp_orbit.tif, -PdemNameExternal DEM, -PexternalDEMFile dem_path, -PpixelSpacingInMeter10.0, -PmapProjectionEPSG:4326, -t, output_path ] subprocess.run(cmd_terrain, checkTrue) return output_path # 调用示例 geometric_correction(gf3_sigma0_filtered.tif, gf3_terrain_corrected.tif, srtm_dem.tif)地形校正的关键参数是pixelSpacingInMeterGF3的常用输出分辨率是10米。DEM的选择直接影响校正精度——SRTM 30米DEM在平坦地区够用山区建议用ASTER GDEM或更高精度的DEM。mapProjection设为EPSG:4326是经纬度投影如果需要等面积或等距离分析改成对应的投影坐标系。如果校正后的影像出现明显的几何畸变检查轨道文件是否匹配。GF3的轨道文件需要从数据分发方获取用错轨道文件会导致系统性偏移。我遇到过用相邻轨道的轨道文件做校正结果影像整体偏移了200多米。4. GF3影像处理的避坑与常见问题排查4.1 定标后影像全黑或全白现象辐射定标后影像显示全黑或全白直方图集中在两端。原因定标常数用错量级或者DN值本身是16位而读取时按8位处理了。解决先检查原始DN值的范围用np.percentile看1%和99%分位数。如果DN值范围在0到65535之间说明是16位数据读取时要指定dtypenp.uint16。定标常数从元数据中核对不要用其他卫星的常数替代。4.2 斑点滤波后边缘出现亮线现象滤波后影像中道路、建筑边缘出现亮线或暗线。原因Lee滤波的窗口跨越了强边缘局部均值被边缘两侧的不同地物拉偏。解决改用边缘保持更好的滤波器比如Frost滤波或非局部均值滤波。或者先做边缘检测在边缘区域减小滤波窗口。我一般会在滤波前先做一次3×3的中值滤波把椒盐噪声去掉再用Lee滤波边缘亮线会明显减轻。4.3 地形校正后山区出现黑色阴影现象地形校正后山体背坡面出现大面积黑色区域。原因地形校正时入射角计算错误或者DEM与影像配准有偏移。解决检查DEM和影像的配准精度用同名点验证。如果偏移超过一个像元先做DEM和影像的配准再校正。另外背坡面的雷达阴影是物理现象不是校正错误——这些区域本来就没有回波信号校正后仍然是无效值不要试图填充。4.4 批量处理时内存溢出现象处理大区域多景影像时程序崩溃报MemoryError。原因一次性把整景影像读入内存GF3单景数据量在1GB以上多景叠加就爆了。解决用分块处理策略每次只读一个窗口。def block_process(input_path, output_path, block_size1024): 分块处理大影像避免内存溢出 with rasterio.open(input_path) as src: profile src.profile.copy() profile.update(dtyperasterio.float32) with rasterio.open(output_path, w, **profile) as dst: for i in range(0, src.height, block_size): for j in range(0, src.width, block_size): # 计算当前块的范围 window rasterio.windows.Window( j, i, min(block_size, src.width - j), min(block_size, src.height - i) ) # 只读取当前块 block src.read(1, windowwindow).astype(np.float32) # 处理当前块 processed lee_filter(block, window_size7) # 写入输出 dst.write(processed, 1, windowwindow)分块大小的选择1024×1024是经验值太小则I/O次数多太大则内存占用高。如果处理环境内存充足32GB以上可以调到2048×2048。4.5 坐标系统不匹配导致叠加错位现象GF3影像和光学影像或矢量数据叠加时错位。原因GF3地理编码产品默认用WGS84地理坐标系而光学影像常用UTM投影坐标系。解决统一坐标系后再叠加。用rasterio的warp功能做投影转换from rasterio.warp import calculate_default_transform, reproject, Resampling def reproject_to_utm(src_path, dst_path, dst_crsEPSG:32650): 将GF3影像重投影到UTM坐标系 with rasterio.open(src_path) as src: transform, width, height calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds ) kwargs src.meta.copy() kwargs.update({ crs: dst_crs, transform: transform, width: width, height: height }) with rasterio.open(dst_path, w, **kwargs) as dst: for i in range(1, src.count 1): reproject( sourcerasterio.band(src, i), destinationrasterio.band(dst, i), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, resamplingResampling.bilinear )重采样方法选bilinear适合连续数据如果是分类结果用nearest保留原始类别值。5. 从处理完的影像到可交付成果进阶技巧与验证方法处理完的GF3影像怎么判断质量合不合格我一般用三个指标来验证。第一是辐射一致性选影像中一块均匀区域计算其均值和标准差标准差越小说明斑点滤波效果越好但均值不能有明显偏移。第二是几何精度用至少5个检查点验证校正后的位置误差平坦地区应该在1到2个像元以内山区可以放宽到3个像元。第三是地形校正效果比较校正前后山体阴影区的后向散射系数变化校正后同一坡向的均值应该趋于一致。进阶用法里我比较常用的是多时相堆叠分析。把同一区域不同时间的GF3影像处理成统一的后向散射系数尺度然后做时序分析。这里有个容易忽略的点不同时间的入射角不同即使做了入射角校正残余的入射角差异仍然会影响时序对比。我的做法是选入射角差异在3度以内的影像做时序超过3度的单独处理。def time_series_stack(image_list, output_path): 多时相GF3影像堆叠用于时序分析 import rasterio from rasterio.merge import merge # 读取所有影像 src_files [rasterio.open(f) for f in image_list] # 合并为多波段影像 mosaic, out_transform merge(src_files) # 写入输出 out_meta src_files[0].meta.copy() out_meta.update({ driver: GTiff, height: mosaic.shape[1], width: mosaic.shape[2], transform: out_transform, count: len(src_files) }) with rasterio.open(output_path, w, **out_meta) as dst: dst.write(mosaic) # 关闭文件句柄 for src in src_files: src.close() return output_path堆叠后的多波段影像可以直接输入到时序分析算法中比如变化检测、趋势分析等。如果要做分类建议先做特征提取把后向散射系数、纹理特征、极化分解参数组合成特征向量再输入分类器。最后一个技巧关于成果输出GF3处理完的影像如果要做成图件建议在输出前做一次直方图拉伸用2%到98%的截断拉伸视觉效果比线性拉伸好得多。如果要做定量分析保留原始dB值不要拉伸。我自己的习惯是每处理完一景影像先存一份原始dB值的GeoTIFF作为存档再另存一份拉伸后的用于目视检查。这样后面要做定量分析时不用重新跑一遍流程。这个习惯帮我省过好几次返工的时间——有一次做变化检测时发现拉伸后的影像对比度不对回头用存档的原始数据重新分析发现是拉伸参数设错了原始数据没问题。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/10 6:53:24 · 阅读 1.2k GF3卫星影像处理全流程：从原始数据到可用成果的完整链路
简介这份资源围绕GF3卫星影像的完整处理流程展开面向遥感影像处理初学者及需要上手SARscape工具的从业者帮助解决从原始数据到成果输出各环节的参数设置与操作规范问题。包内共1个docx文档压缩包约4.58MB以图文步骤形式组织内容便于对照软件界面逐步操作。文档依次覆盖系统设置、数据导入、多视处理、单通道强度数据滤波等核心环节并延伸至地理编码与辐射定标对极化方式选择、制图分辨率设定、Frost滤波窗口参数等关键配置均有说明还给出入射角、采样间隔等元数据读取与分辨率换算思路。目前已有1577人学习下载适合希望系统掌握GF3影像预处理流程、减少参数试错成本的读者参考。1. GF3卫星影像处理从原始数据到可用成果的完整链路拿到一景GF3卫星影像很多人第一反应是直接拖进软件里开始调色。但如果真这么干大概率会在某个环节卡住——要么是数据打不开要么是坐标对不上要么是做完正射校正后发现边缘区域拉伸得没法看。GF3作为C波段合成孔径雷达卫星它的影像处理和光学影像完全是两套逻辑。光学影像处理的核心是辐射校正和大气校正而GF3影像处理的核心是聚焦、辐射定标、斑点滤波、几何校正和地形校正这一整条链路。这篇文章面向的是已经拿到GF3数据、需要把它变成可分析成果的从业者不管你是做地表形变监测、地物分类还是变化检测这套处理流程都是绕不开的基本功。我会按实际操作顺序把每一步的参数设置、常见报错和踩坑经验讲清楚。2. GF3影像处理的前置准备与数据理解2.1 GF3数据的产品分级与选型逻辑GF3卫星影像按处理程度分为多个级别常见的有原始回波数据、单视复数数据、地距探测数据、地理编码数据等。不同级别的数据对应不同的处理起点选错了级别要么做重复工作要么根本没法往下走。产品级别内容适用场景处理起点原始回波未聚焦的原始信号需要自定义聚焦算法聚焦处理单视复数聚焦后的复数数据干涉测量、极化分析辐射定标滤波地距探测检测后的幅度数据强度分析、分类滤波几何校正地理编码投影到地理坐标直接叠加分析精校正后处理我一般建议如果你的目标是做地物分类或变化检测直接从地距探测数据开始最省事如果要做干涉测量必须从单视复数数据起步。选错级别是新手最容易翻车的地方——拿地理编码数据去做干涉等于把已经投影过的数据再投影一次精度损失不可逆。2.2 处理环境的搭建与依赖检查GF3影像处理常用的工具链包括欧空局的开源工具箱、商业遥感软件以及Python生态中的遥感处理库。不管用哪套工具环境搭建时最容易忽略的是依赖库版本冲突。# 以Python生态为例创建独立环境避免版本冲突 conda create -n gf3_process python3.9 conda activate gf3_process # 安装核心依赖注意版本锁定 pip install numpy1.23.5 pip install scipy1.10.1 pip install gdal3.6.2 pip install rasterio1.3.6 pip install snapista1.0.0 # 用于调用欧空局工具箱的Python接口这段脚本的关键在于版本锁定。GDAL和rasterio的版本必须匹配否则读写GF3的GeoTIFF时会报“not recognized as a supported file format”。snapista是调用欧空局工具箱的Python封装如果你习惯用图形界面操作可以跳过这一步但批量处理时命令行效率高得多。参数说明python3.9是经过验证的稳定版本3.10以上在部分遥感库中会有兼容性问题。numpy锁定1.23.5是因为更高版本与scipy 1.10.1存在ABI冲突。gdal 3.6.2对GF3的元数据解析支持最完整。注意如果你在Windows下操作GDAL的安装建议用conda而不是pippip安装的GDAL经常找不到proj数据库路径。2.3 元数据解析与轨道信息提取GF3影像的元数据文件里藏着处理所需的关键参数成像时间、轨道方向、入射角、极化方式、采样间隔等。这些参数直接决定后续辐射定标和几何校正的系数。import xml.etree.ElementTree as ET def parse_gf3_metadata(xml_path): 解析GF3元数据文件提取关键处理参数 tree ET.parse(xml_path) root tree.getroot() # 命名空间处理GF3元数据通常带命名空间 ns {ns: root.tag.split(})[0].strip({)} if } in root.tag else {} meta {} # 提取成像时间 meta[acquisition_time] root.findtext(.//ns:imagingTime, namespacesns) # 提取轨道方向 meta[orbit_direction] root.findtext(.//ns:orbitDirection, namespacesns) # 提取中心入射角 meta[incidence_angle] float(root.findtext(.//ns:incidenceAngle, namespacesns)) # 提取极化方式 meta[polarization] root.findtext(.//ns:polarization, namespacesns) # 提取距离向和方位向采样间隔 meta[range_spacing] float(root.findtext(.//ns:rangeSpacing, namespacesns)) meta[azimuth_spacing] float(root.findtext(.//ns:azimuthSpacing, namespacesns)) return meta # 调用示例 meta parse_gf3_metadata(GF3_XXX_meta.xml) print(f入射角: {meta[incidence_angle]}度) print(f极化方式: {meta[polarization]})这段代码的逻辑是遍历XML树按标签名提取参数。命名空间处理是关键——GF3元数据文件通常带命名空间前缀直接按标签名查找会返回None。入射角参数用于后续辐射定标中的入射角校正采样间隔用于计算像元大小和几何校正的网格间距。如果解析报错“标签未找到”先用文本编辑器打开XML文件确认实际的标签名和命名空间。不同批次的数据元数据格式可能有细微差异我遇到过同一颗卫星不同年份的数据标签名从“incidenceAngle”变成“incidence_angle”的情况。3. GF3影像的核心处理步骤与参数设置3.1 辐射定标把DN值转成后向散射系数GF3影像的原始像元值是DN值没有物理意义。辐射定标就是把DN值转换成后向散射系数让不同时间、不同轨道的影像可以对比。import numpy as np import rasterio def radiometric_calibration(dn_path, output_path, calib_factor, incidence_angle): GF3辐射定标DN值转后向散射系数 with rasterio.open(dn_path) as src: dn src.read(1).astype(np.float32) profile src.profile.copy() # 避免除零 dn[dn 0] np.nan # 定标公式sigma0 (DN^2 / calib_factor) * sin(incidence_angle) # calib_factor来自元数据中的定标常数 sigma0 (dn ** 2) / calib_factor # 入射角校正 sigma0 sigma0 * np.sin(np.radians(incidence_angle)) # 转dB sigma0_db 10 * np.log10(sigma0) # 写入输出 profile.update(dtyperasterio.float32, count1, nodatanp.nan) with rasterio.open(output_path, w, **profile) as dst: dst.write(sigma0_db, 1) return sigma0_db # 调用示例 calib_factor 1.234e5 # 从元数据中读取 sigma0 radiometric_calibration(gf3_dn.tif, gf3_sigma0.tif, calib_factor, 35.2)定标公式中的calib_factor是定标常数必须从元数据中读取不能自己猜。入射角校正这一步很多人会漏掉——如果不做入射角校正同一景影像中近距端和远距端的后向散射系数会有系统性偏差做分类时边缘区域会大量误分。参数说明calib_factor的量级通常在10^4到10^6之间具体值取决于卫星的发射功率和增益设置。入射角从元数据中提取单位是度。转dB是为了压缩动态范围方便后续滤波和显示。注意如果定标后的影像出现大量NaN检查原始DN值中是否有零值区域。GF3影像的边缘区域经常有零值填充这些区域在定标前就要标记为无效。3.2 斑点滤波在去噪和保边之间找平衡合成孔径雷达影像天生有斑点噪声这是相干成像的物理机制决定的。斑点滤波的目标是抑制噪声的同时保留边缘和纹理信息。常用的滤波器有Lee滤波、Frost滤波、Gamma MAP滤波等。from scipy.ndimage import uniform_filter, median_filter def lee_filter(image, window_size7, noise_varNone): Lee滤波基于局部统计的自适应斑点滤波 if noise_var is None: # 估计噪声方差通常用均匀区域的方差 noise_var np.nanvar(image[100:200, 100:200]) # 局部均值和方差 local_mean uniform_filter(image, window_size) local_sq_mean uniform_filter(image ** 2, window_size) local_var local_sq_mean - local_mean ** 2 # Lee滤波权重 weight local_var / (local_var noise_var) weight np.clip(weight, 0, 1) # 滤波结果 filtered local_mean weight * (image - local_mean) return filtered # 调用示例 filtered lee_filter(sigma0, window_size7)Lee滤波的核心思想是用局部方差和噪声方差的比值来决定平滑程度。在均匀区域局部方差接近噪声方差权重趋近于0.5滤波强度大在边缘区域局部方差远大于噪声方差权重趋近于1滤波强度小边缘得以保留。窗口大小的选择是个经验活。7×7是常用值窗口太小去噪不干净窗口太大边缘模糊。我一般会先用7×7跑一遍看效果如果斑点残留明显再调到9×9。noise_var的估计很关键——估大了边缘模糊估小了去噪不足。建议选影像中一块均匀区域比如平静水面或裸地来计算方差。3.3 几何校正与地形校正让影像和地图对得上GF3影像的几何校正分两步先做系统几何校正用轨道参数和地球模型把影像投影到地理坐标再做精校正用地面控制点或参考影像消除残余误差。地形校正则是消除地形起伏引起的辐射畸变。def geometric_correction(input_path, output_path, dem_path, gcp_listNone): GF3几何校正与地形校正 import subprocess # 使用欧空局工具箱的命令行接口 # 第一步应用轨道文件进行系统几何校正 cmd_orbit [ gpt, Apply-Orbit-File, -Ssource input_path, -PorbitTypeSentinel Precise, -t, temp_orbit.tif ] subprocess.run(cmd_orbit, checkTrue) # 第二步地形校正需要外部DEM cmd_terrain [ gpt, Terrain-Correction, -Ssourcetemp_orbit.tif, -PdemNameExternal DEM, -PexternalDEMFile dem_path, -PpixelSpacingInMeter10.0, -PmapProjectionEPSG:4326, -t, output_path ] subprocess.run(cmd_terrain, checkTrue) return output_path # 调用示例 geometric_correction(gf3_sigma0_filtered.tif, gf3_terrain_corrected.tif, srtm_dem.tif)地形校正的关键参数是pixelSpacingInMeterGF3的常用输出分辨率是10米。DEM的选择直接影响校正精度——SRTM 30米DEM在平坦地区够用山区建议用ASTER GDEM或更高精度的DEM。mapProjection设为EPSG:4326是经纬度投影如果需要等面积或等距离分析改成对应的投影坐标系。如果校正后的影像出现明显的几何畸变检查轨道文件是否匹配。GF3的轨道文件需要从数据分发方获取用错轨道文件会导致系统性偏移。我遇到过用相邻轨道的轨道文件做校正结果影像整体偏移了200多米。4. GF3影像处理的避坑与常见问题排查4.1 定标后影像全黑或全白现象辐射定标后影像显示全黑或全白直方图集中在两端。原因定标常数用错量级或者DN值本身是16位而读取时按8位处理了。解决先检查原始DN值的范围用np.percentile看1%和99%分位数。如果DN值范围在0到65535之间说明是16位数据读取时要指定dtypenp.uint16。定标常数从元数据中核对不要用其他卫星的常数替代。4.2 斑点滤波后边缘出现亮线现象滤波后影像中道路、建筑边缘出现亮线或暗线。原因Lee滤波的窗口跨越了强边缘局部均值被边缘两侧的不同地物拉偏。解决改用边缘保持更好的滤波器比如Frost滤波或非局部均值滤波。或者先做边缘检测在边缘区域减小滤波窗口。我一般会在滤波前先做一次3×3的中值滤波把椒盐噪声去掉再用Lee滤波边缘亮线会明显减轻。4.3 地形校正后山区出现黑色阴影现象地形校正后山体背坡面出现大面积黑色区域。原因地形校正时入射角计算错误或者DEM与影像配准有偏移。解决检查DEM和影像的配准精度用同名点验证。如果偏移超过一个像元先做DEM和影像的配准再校正。另外背坡面的雷达阴影是物理现象不是校正错误——这些区域本来就没有回波信号校正后仍然是无效值不要试图填充。4.4 批量处理时内存溢出现象处理大区域多景影像时程序崩溃报MemoryError。原因一次性把整景影像读入内存GF3单景数据量在1GB以上多景叠加就爆了。解决用分块处理策略每次只读一个窗口。def block_process(input_path, output_path, block_size1024): 分块处理大影像避免内存溢出 with rasterio.open(input_path) as src: profile src.profile.copy() profile.update(dtyperasterio.float32) with rasterio.open(output_path, w, **profile) as dst: for i in range(0, src.height, block_size): for j in range(0, src.width, block_size): # 计算当前块的范围 window rasterio.windows.Window( j, i, min(block_size, src.width - j), min(block_size, src.height - i) ) # 只读取当前块 block src.read(1, windowwindow).astype(np.float32) # 处理当前块 processed lee_filter(block, window_size7) # 写入输出 dst.write(processed, 1, windowwindow)分块大小的选择1024×1024是经验值太小则I/O次数多太大则内存占用高。如果处理环境内存充足32GB以上可以调到2048×2048。4.5 坐标系统不匹配导致叠加错位现象GF3影像和光学影像或矢量数据叠加时错位。原因GF3地理编码产品默认用WGS84地理坐标系而光学影像常用UTM投影坐标系。解决统一坐标系后再叠加。用rasterio的warp功能做投影转换from rasterio.warp import calculate_default_transform, reproject, Resampling def reproject_to_utm(src_path, dst_path, dst_crsEPSG:32650): 将GF3影像重投影到UTM坐标系 with rasterio.open(src_path) as src: transform, width, height calculate_default_transform( src.crs, dst_crs, src.width, src.height, *src.bounds ) kwargs src.meta.copy() kwargs.update({ crs: dst_crs, transform: transform, width: width, height: height }) with rasterio.open(dst_path, w, **kwargs) as dst: for i in range(1, src.count 1): reproject( sourcerasterio.band(src, i), destinationrasterio.band(dst, i), src_transformsrc.transform, src_crssrc.crs, dst_transformtransform, dst_crsdst_crs, resamplingResampling.bilinear )重采样方法选bilinear适合连续数据如果是分类结果用nearest保留原始类别值。5. 从处理完的影像到可交付成果进阶技巧与验证方法处理完的GF3影像怎么判断质量合不合格我一般用三个指标来验证。第一是辐射一致性选影像中一块均匀区域计算其均值和标准差标准差越小说明斑点滤波效果越好但均值不能有明显偏移。第二是几何精度用至少5个检查点验证校正后的位置误差平坦地区应该在1到2个像元以内山区可以放宽到3个像元。第三是地形校正效果比较校正前后山体阴影区的后向散射系数变化校正后同一坡向的均值应该趋于一致。进阶用法里我比较常用的是多时相堆叠分析。把同一区域不同时间的GF3影像处理成统一的后向散射系数尺度然后做时序分析。这里有个容易忽略的点不同时间的入射角不同即使做了入射角校正残余的入射角差异仍然会影响时序对比。我的做法是选入射角差异在3度以内的影像做时序超过3度的单独处理。def time_series_stack(image_list, output_path): 多时相GF3影像堆叠用于时序分析 import rasterio from rasterio.merge import merge # 读取所有影像 src_files [rasterio.open(f) for f in image_list] # 合并为多波段影像 mosaic, out_transform merge(src_files) # 写入输出 out_meta src_files[0].meta.copy() out_meta.update({ driver: GTiff, height: mosaic.shape[1], width: mosaic.shape[2], transform: out_transform, count: len(src_files) }) with rasterio.open(output_path, w, **out_meta) as dst: dst.write(mosaic) # 关闭文件句柄 for src in src_files: src.close() return output_path堆叠后的多波段影像可以直接输入到时序分析算法中比如变化检测、趋势分析等。如果要做分类建议先做特征提取把后向散射系数、纹理特征、极化分解参数组合成特征向量再输入分类器。最后一个技巧关于成果输出GF3处理完的影像如果要做成图件建议在输出前做一次直方图拉伸用2%到98%的截断拉伸视觉效果比线性拉伸好得多。如果要做定量分析保留原始dB值不要拉伸。我自己的习惯是每处理完一景影像先存一份原始dB值的GeoTIFF作为存档再另存一份拉伸后的用于目视检查。这样后面要做定量分析时不用重新跑一遍流程。这个习惯帮我省过好几次返工的时间——有一次做变化检测时发现拉伸后的影像对比度不对回头用存档的原始数据重新分析发现是拉伸参数设错了原始数据没问题。希望帮到你。本文还有配套的精品资源点击获取