国产卫星定标系数与光谱响应函数：6S大气校正实操指南

简介：这套压缩包聚焦国产卫星遥感数据处理，面向从事高分（GF）系列、资源（ZY）系列等国产卫星影像处理的研究人员与开发者，解决定标系数与光谱响应函数获取难、格式不统一的问题，并支持与6S辐射传输模型结合开展大气校正和地表参数反演。包内共93个文件，约1.85MB，涵盖txt、json、xls、xlsx、doc、pdf等多种格式：txt与json提供光谱响应函数及可直接引用的参数配置，xls/xlsx保存各传感器的详细响应曲线，doc/pdf为历年绝对辐射定标系数报告，jpg为辅助元数据文件。资源按卫星型号分目录整理，覆盖GF-1/2/4/6、ZY-3、ZY-02C、HJ-1A/B等常用载荷，便于按需检索和读取。目前已有2950人学习下载，适合遥感、GIS方向的学生和工程师用于数据预处理、模型参数配置及科研项目中的精度校验。

1. 国产卫星定标系数和光谱响应函数：从“能看图”到“能量算”的第一步

做遥感定量化的人心里都有个数：国产高分辨率卫星影像下下来，DN值看起来几十到几百，看着舒服，但真要算地表反射率、植被指数、水体叶绿素浓度，第一步就是辐射定标。没有合适的定标系数，DN值就是一堆“数字玄学”；没有光谱响应函数，6S大气校正模型根本认不出这台国产传感器，只能瞎填参数。这两样东西，正是这个资源包的核心，它把某国产光学卫星的绝对辐射定标系数和分波段光谱响应函数整理成可直接读取的文件，配合6S模型做大气校正，能少走很多弯路。适合被国产卫星数据“折磨”过的遥感工程师、做定量反演的研究生，以及要给甲方交付反射率产品的数据处理团队。

2. 定标系数与光谱响应函数的底层逻辑：DN值、辐亮度与等效波段

2.1 绝对辐射定标系数：DN值到辐亮度的线性映射

先明确一个前提：国产光学卫星装载的CCD/CMOS探测器记录的是无量纲的数字量化值DN（Digital Number），它跟入瞳辐亮度之间不直接相等，中间隔着一组定标系数。常见公式是：

L = Gain × DN + Offset

其中L是入瞳等效辐亮度（单位W·m⁻²·sr⁻¹·μm⁻¹），Gain是增益，Offset是本底偏置。这组系数通常由定标场实验结合星上定标数据联合给出，随卫星在轨时间推移会发生变化——探测器老化、光学系统透过率下降，都会让同一地物在不同年份显示出不同DN值。用旧系数处理新影像，反演出的反射率会系统性偏高或偏低，时序分析时就会看到“假趋势”。

另一个方向是，有的数据产品说明里把公式写成L = (DN - Offset) / Gain，这种写法在国产卫星里不少见，尤其是线阵推扫传感器分档增益输出的情况下。拿到压缩包里的定标系数文件时，第一件事就是确认公式定义。我在某次处理一批历史影像时就吃过这个亏：系数文件表头写的是“Gain”和“Bias”，我按L = Gain × DN + Bias算完，辐亮度整体偏了大概30%，后来才知道该批次数据的说明文件里定义的是L = (DN - Bias) / Gain。别人用这批数据反演NDVI还看不出大问题，但算水温、算悬浮物浓度就直接偏到没法用。

数据包里的定标系数文件，通常按传感器、按年份组织，一个波段对应一组Gain和Bias。如果你处理的是多期影像，就要设计一个“按影像日期查表取系数”的步骤，而不是整批用同一个数值。这个我放到第三章用代码细说。

2.2 光谱响应函数的物理含义：每个波段不是一条“线”

很多初学者会把“波段”理解成一条单一波长的光线，其实每一个多光谱波段都覆盖一段波长范围。以常见的国产多光谱载荷为例，蓝波段大约覆盖450nm到520nm，近红外波段可能从760nm到890nm。在这段区间内，探测器对每个波长的敏感程度并不相同——光学镜片镀膜、分光棱镜、探测器量子效率共同决定了一条“通道响应曲线”，这就是光谱响应函数（Spectral Response Function，SRF）。

SRF最常见的表示是一个二维数据对：波长（nm）和对应响应值（0到1之间）。响应值越大，说明该波长的能量对信号贡献越多。它有两个关键用途：一是计算波段的等效太阳辐照度ESUN，二是作为6S等辐射传输模型的自定义传感器输入。为什么说“每个波段不是一条线”？因为如果忽略SRF，把波段当成某个中心波长下的单色光，计算出来的ESUN会偏差几个百分点，这个偏差直接进入反射率反演的分母，导致最终地表反射率系统性偏移。对于植被、水体这些反射率绝对值本来就不大的地物，这种误差会掩盖真实信号。

SRF文件的归一化方式也得注意。有的卫星官方发布的SRF峰值归一化到1.0，有的按积分面积归一化到某个数值。峰值归一化用起来直观，但计算等效太阳辐照度时，直接用峰值归一化的SRF乘上太阳光谱再积分，结果是相对值，还要做面积修正；面积归一化的SRF则可以“一步到位”直接参与积分。具体用哪种，取决于你是自己写ESUN计算代码，还是把SRF喂给6S去内部处理。

2.3 数据包的文件清单：解开压缩包先核对什么

这类资源包下载下来后，通常会解出四类主要文件：

我的习惯是拿到压缩包后先不着急跑代码，先花两分钟做三件事：第一，打开说明文档，确认定标系数公式用的是“L = Gain × DN + Offset”还是“L = (DN - Offset) / Gain”，这个错了后面全错。第二，看SRF文件波长列的单位是纳米还是微米，同一份数据单位写错，6S输入直接“找不到波段”。第三，看有没有按卫星型号分目录，避免把A星系数套在B星数据上。这一步虽说不难，但能救不少“数据算出来不对”的血泪现场。

提示：解压后如果看到多个年份的定标系数子目录，优先看最近一年的更新说明——国产卫星定标系数会随在轨时间调整，年份越新通常越接近该阶段真实传感器状态。

3. 辐射定标实操：从定标系数到表观反射率的一行代码与参数细节

3.1 用Python读取定标系数并换算辐亮度

到现在，资源包里的定标系数已经能落进代码了。我用Python做遥感数据处理时，习惯先把定标系数读成结构化的DataFrame，再按影像元数据中的日期来匹配对应系数。用一个实际的处理流程来展示：

import pandas as pd import numpy as np # 读取定标系数表，假设格式：sensor, band, gain, offset, valid_year cal_df = pd.read_csv('calibration_coefficient.csv') # 某景影像元数据信息 img_date = '2022-11-08' sensor_name = 'GX-1_PMS' # 按传感器和年份匹配定标系数 cal_row = cal_df[ (cal_df['sensor'] == sensor_name) & (cal_df['valid_year'] == int(img_date[:4])) ] if len(cal_row) == 0: raise ValueError('没有匹配到该年份的定标系数，请检查系数有效期') # 假设band1(蓝)的Gain=0.168, Offset=0.0012 gain = cal_row[cal_row['band'] == 'band1']['gain'].values[0] offset = cal_row[cal_row['band'] == 'band1']['offset'].values[0] print(f"band1 使用 gain={gain}, offset={offset}")

这段代码里最关键的是“按年份匹配”这个动作。定标系数通常只在某段时间内有效，跨年使用会引入系统偏差。输出一行使用参数，是为了后面排查问题时能确认“这景影像到底用了哪组系数”。常见做法是把它写进日志文件，处理大批量数据时回头能找到历史记录。

得到Gain和Offset后，辐射定标只是一行广播运算。假设原始DN影像已经读入numpy数组：

def dn_to_radiance(dn_array, gain, offset): """ DN值转辐亮度 L = Gain * DN + Offset """ radiance = gain * dn_array + offset # 对无效像元保持无效值 radiance[dn_array <= 0] = np.nan return radiance # 示例：读取band1的DN dn_b1 = np.random.randint(50, 200, size=(512, 512)) rad_b1 = dn_to_radiance(dn_b1, gain=0.168, offset=0.0012) print(f"DN范围: {dn_b1.min()}~{dn_b1.max()}") print(f"辐亮度范围: {rad_b1.min():.3f}~{rad_b1.max():.3f} W/m2/sr/um")

辐亮度值级判断有个经验：晴朗条件下，陆地多光谱各波段入瞳辐亮度通常在几十到两百W·m⁻²·sr⁻¹·μm⁻¹之间，如果算出几千甚至数万，先回头检查Gain和Offset是否用反，或者DN是不是没做转换。把无效像元直接置NaN能有效防止后续统计时被0值干扰。

3.2 从辐亮度到表观反射率：ESUN、太阳天顶角与日地距离

辐亮度是绝对物理量，但不同影像之间的太阳光照条件不一样，直接拿辐亮度做对比会产生假差异。业内统一做法是把辐亮度转成表观反射率（Top of Atmosphere Reflectance）：

import numpy as np def radiance_to_reflectance(radiance, esun, solar_zenith_deg, earth_sun_distance): """ 表观反射率 = pi * L * d^2 / (ESUN * cos(theta)) :param esun: 波段等效太阳辐照度 (W/m2/um) :param solar_zenith_deg: 太阳天顶角 (度) :param earth_sun_distance: 日地距离 (天文单位AU) """ cos_theta = np.cos(np.deg2rad(solar_zenith_deg)) reflectance = np.pi * radiance * (earth_sun_distance ** 2) / (esun * cos_theta) return reflectance # 假设band1的ESUN = 1960 W/m2/um，太阳天顶角34.2度，日地距离1.011 refl_b1 = radiance_to_reflectance(rad_b1, esun=1960, solar_zenith_deg=34.2, earth_sun_distance=1.011)

这里的ESUN怎么来？最稳妥的办法是：把资源包里的SRF与太阳光谱数据（如Wehrli或ASTM E-490标准太阳光谱）做卷积：

def calc_esun_from_srf(sun_wavelength_nm, sun_irradiance, srf_wavelength_nm, srf_value): """ 按SRF加权计算等效太阳辐照度 ESUN = ∫(E0(λ) * SRF(λ) dλ) / ∫SRF(λ) dλ """ # 把太阳光谱插值到SRF波长网格 srf_interp_e0 = np.interp(srf_wavelength_nm, sun_wavelength_nm, sun_irradiance) numerator = np.trapezoid(srf_interp_e0 * srf_value, srf_wavelength_nm) denominator = np.trapezoid(srf_value, srf_wavelength_nm) return numerator / denominator

注意：这个公式的前提是SRF已经是面积归一的权重，如果用的是峰值归一化SRF，分母的积分就是纯面积修正，结果等价。计算时波长的单位要保持一致，要么全是纳米，要么全是微米，混合使用会让ESUN大出或小出几个数量级。

trapezoid在旧版本numpy里叫trapz，新版本推荐trapezoid。如果环境里用的是老版本的NumPy，先把函数名改了再跑，这个坑我踩过，浪费了一个下午。

3.3 跨年份影像的定标系数匹配策略

数据包里通常不止一年的定标系数，实操时得维护一张“影像日期→定标系数”的映射表。我一般会用文件名里的日期字段做自动匹配，规则如下：

这个匹配动作放在批处理主循环里，比手工改参数安全得多。处理上百景影像时，一旦某景混用了错误系数，最终产品上看起来就是一条“碎掉”的反射率曲线，很难排查。我在第三章前面写了日志输出，就是为了在这种时候能快速定位是哪一景、用了哪组系数。

4. 光谱响应函数在6S模型中的接入：重采样与格式转换

4.1 6S自定义传感器：内置表里没有国产卫星

6S（Second Simulation of the Satellite Signal in the Solar Spectrum）是太阳反射波段最常用的大气辐射传输模型，做大气校正时它要求输入“传感器光谱响应”才能算出该波段的大气透过率和散射路径辐射。6S内置了一批常见传感器参数，比如Landsat TM、SPOT HRV，但没有国产高分辨率卫星。很多人在这一步卡住：要么把波段中心波长填了了事，要么找一个相近的传感器“假装”成国产卫星，结果就是大气校正结果在近红外波段明显偏绿——本质上是光谱响应差异带来的系统性误差。

正确做法是走6S的自定义传感器通道，把资源包里的SRF转成6S能读的格式。6S要求光谱响应按波长顺序逐点列出，波长范围要覆盖该波段有效响应的全部区间，步长要尽量均匀。常见的步长是10nm，也有用5nm或2.5nm的，步长越小越精细，但6S计算时间会变长。实际工程里，带宽在50nm以上的多光谱波段，用10nm步长已经足够稳定，没必要追求过细。

4.2 SRF重采样：把国产卫星数据转成6S需要的格式

原始SRF数据的波长间隔通常不均匀，官方发布文件里可能每隔0.5nm或1nm给一个响应值，也可能只在关键波长附近有数据点。6S输入需要一个固定步长的响应序列，所以先进先得做线性插值重采样：

import numpy as np from scipy.interpolate import interp1d # 原始SRF，两列：波长(nm)、响应值 srf_raw = np.loadtxt('srf_GX1_band1.txt') # 确定重采样范围：从有效响应(>0.005)的首尾取整数 valid = srf_raw[srf_raw[:, 1] > 0.005] wl_start = int(np.floor(valid[:, 0].min())) wl_end = int(np.ceil(valid[:, 0].max())) # 10nm步长 wl_new = np.arange(wl_start, wl_end + 10, 10, dtype=float) # 线性插值 interp_func = interp1d(srf_raw[:, 0], srf_raw[:, 1], kind='linear', bounds_error=False, fill_value=0.0) srf_interp = interp_func(wl_new) # 峰值归一化（6S内部会再处理，但统一到1.0更安全） srf_interp = srf_interp / np.max(srf_interp) # 输出成6S格式文本：每行一个响应值 np.savetxt('srf_GX1_band1_10nm.txt', srf_interp, fmt='%.5f', header='', comments='')

这段代码有两个参数要想清楚：一是“有效响应”的阈值，我取0.005，小于这个响应值说明该波长对信号的贡献小于0.5%，直接截断不会造成明显误差；二是步长10nm，与6S内置传感器的采样习惯保持一致。如果你后续要计算ESUN，重采样后的SRF最好再单独存一份带波长列的版本，供自己算ESUN用。

输出文件里只有响应值、没有波长列，这是6S的一种常见输入约定——程序会依据你指定的波长起止和步长来理解这一串数值。因此，文件名里一定要带上波段信息和步长，例如srf_GX1_band1_10nm.txt，否则几个月后再看文件根本想不起来里面是什么。

4.3 6S大气校正完整输入示例与输出解读

有了10nm步长的SRF文件，接下来把参数写进6S输入文件。6S的输入行比较多，我把常用参数整理成一份清单：

6S跑完会输出一个独立的结果文件，核心是“实测表观反射率”和“校正后的地表反射率”两个量，以及大气透过率、路径辐射等中间量。要注意的是，6S输出的反射率是0到1之间的小数，不是百分数。有些处理流程把6S输出直接乘以100写成百分比，后续差分、分类都能用，但存档时要保持统一，不然会闹出“反射率1.05”这类笑话。

我在实际操作中，会把6S跑完的结果与同时刻的MODIS反射率产品做交叉对比，如果偏差在0.02以内说明SRF文件和参数基本靠谱；偏差超过0.05，优先回头查SRF文件，因为几何参数和大气参数填错一般不会造成这么大且系统的反射率漂移。

5. 避坑/常见问题：定标系数与光谱响应函数应用的六个高频坑

5.1 定标公式形式看错，辐亮度整体偏移

现象：同一个DN值算出来的辐亮度，和参考数据差30%以上，且所有波段偏差方向一致。原因：说明文档里写的是L = (DN - Offset) / Gain，代码里却按L = Gain × DN + Offset执行。解决：拿到定标系数先看公式定义，再取一个已知DN值手算一遍验证。比如DN=100，Gain=0.168，Offset=0.0012，两种公式算出来的辐亮度分别是16.8012和594.04，差一个量级，一眼就能识别。

5.2 SRF重采样后未归一化，等效太阳辐照度偏小

现象：用重采样后的SRF算出ESUN只有几百W/m2/μm，反射率反演结果整体偏大，暗目标也像亮目标。原因：原始SRF在非均匀波长网格上有较好的积分特性，插值到均匀网格后总面积发生变化，没有重新归一化就参与ESUN积分。解决：重采样后做一次峰值归一化，或做面积归一化（积分值为1），再用归一化后的SRF计算ESUN。我每重采样一次SRF都会顺手把积分面积打印出来，前后对比差超过0.5%就重新检查网格范围。

5.3 6S输入波段范围与SRF有效区间对不上

现象：6S运行报错或输出的反射率在某一波段明显断层。原因：输入的波长起点和终点超出了SRF文件的有效范围，比如SRF文件里波长是450~520nm，6S输入里却填了400~550nm，插值函数把无数据区域补成0，等效于波段范围被强行“拉宽”。解决：用4.2节代码里算出来的wl_start和wl_end直接填进6S输入文件，不要手工从文档里抄波长范围，手工抄错是高频事故。

5.4 定标系数跨年份混用，反射率出现系统性漂移

现象：同一地物在不同年份的反演反射率，始终差一个固定比例，且与真实变化趋势不符。原因：批处理时把最新年份的定标系数直接套在所有影像上，没有按拍摄日期匹配。解决：建一个“影像日期→定标系数”映射表，在批处理循环里按日期查表。查不到对应年份系数时，日志里打Warning并暂停运行，不要自动沿用相邻年份。

5.5 角度单位混乱：cos函数里装了“度”而不是弧度

现象：反射率在高太阳天顶角下偏高，低纬度地区不明显，中高纬地区严重。原因：太阳天顶角是34.2度，Python的math.cos和numpy.cos默认接受弧度，直接cos(34.2)得到的是“弧度制34.2对应的余弦”，不是“34.2度的余弦”。解决：所有角度参与三角函数计算前统一deg2rad()。我给处理脚本定了规矩：所有角度值在读取后立刻转弧度存储，变量名带_rad后缀，这样后面写公式时看到cos_theta_rad就不会再犯单位错误。

5.6 大气校正后负值反射率被直接丢弃

现象：水体、阴影区域的反射率产品出现0值空洞，统计均值被拉低。原因：6S校正后的反射率是辐射传输残差的结果，暗像元校正值可能为负，代码里用refl[refl < 0] = 0直接把负值抹掉，丢失了分布信息。解决：官方产品通常不对负值做任何更改，实际应用时建议把负值保留在结果里，只在统计和可视化时处理。如果业务需要正值产品，用np.clip(refl, 0, None)前先记录负数比例，并在这个比例超过5%时检查输入辐射定标是否准确。

注意：第5.6条的“留负值”和“clip”并不矛盾。前者是存档原始输出，后者是产品后处理；很多团队把两步分开，就是为了能在排查问题时重新审视暗像元分布。

6. 进阶：多波段批量定标与SRF一致性校验的脚本化小技巧

到了这一步，你已经能把单景影像从DN值一路处理到6S大气校正。日常做项目时，真正折磨人的不是单景，是上百景的批次处理。这里分享两个我常用的脚本化技巧。

第一个是多波段批量定标。把3.1和3.2节的单波段逻辑封装成循环：

def batch_processing(dn_cube, cal_params, esun_list, solar_zenith_deg, earth_sun_distance): """ dn_cube: 多波段DN数组 (band, rows, cols) cal_params: 每波段的(gain, offset) esun_list: 每波段的ESUN """ radiance_cube = np.zeros_like(dn_cube, dtype=np.float32) reflectance_cube = np.zeros_like(dn_cube, dtype=np.float32) for band_idx in range(dn_cube.shape[0]): gain, offset = cal_params[band_idx] radiance_cube[band_idx] = gain * dn_cube[band_idx] + offset reflectance_cube[band_idx] = radiance_to_reflectance( radiance_cube[band_idx], esun_list[band_idx], solar_zenith_deg, earth_sun_distance ) return reflectance_cube

批量处理前，我会先用一小块裁剪区域跑通一遍，确认所有波段系数顺序和ESUN列表一一对应。曾经有一版代码在读取定标系数时波段排序是“蓝绿红近红外”，到ESUN列表里却按“红绿蓝近红外”排，结果就是蓝波段用红波段的ESUN，整批影像颜色都变了。从那以后，我每次批处理前都强制走一遍“波段顺序核对”，把各波段的中心波长打印出来跟定标系数表逐项对比，确认无误后再进入全流程。

第二个技巧是SRF一致性校验。数据包里的SRF文件可能来自不同年份，同一颗卫星的光谱响应在轨几年后也可能发生微小漂移。处理多期反演产品前，我会计算相邻两版SRF的交叉重叠面积：

def srf_overlap(srf1_wl, srf1_val, srf2_wl, srf2_val): # 共同波长区间取交集 wl_min = max(srf1_wl.min(), srf2_wl.min()) wl_max = min(srf1_wl.max(), srf2_wl.max()) common_wl = np.arange(wl_min, wl_max, 1.0) s1 = np.interp(common_wl, srf1_wl, srf1_val) s2 = np.interp(common_wl, srf2_wl, srf2_val) return np.trapezoid(np.minimum(s1, s2), common_wl)

重叠面积超过95%，说明两版SRF差异可以忽略，跨年份产品可以直接对比；低于90%，我建议把产品分成两个系列分别处理，否则混在一起出的反射率时序产品会在换年份的时刻出现一次“跳台阶”。

这套资源包给我的价值，不只是省了四处找系数的时间，更重要的是把“从DN到反射率”这条路上最容易出错的环节全部变成了可查的文件和可跑的脚本。希望这份拆解能帮你少踩几个我踩过的坑，让国产卫星数据的定量化处理更顺一点。

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