2019, Vol. 39
2. 北京航天宏图信息技术股份有限公司，北京 100195
2. Beijing PIESAT Information Technology Co. Ltd., Beijing 100195, China
多光谱遥感影像彩色合成方法主要有2种：真彩色合成和假彩色合成。真彩色合成一般是指将多光谱影像的红、绿、蓝波段数据与R、G、B通道一一对应合成，合成后的遥感影像上地物的颜色与实际地物颜色接近或一致；R、G、B通道没有对应影像红、绿、蓝波段数据的彩色合成则为假彩色合成，合成后的影像颜色与实际地物颜色差异较大[1]。真彩色合成过程不涉及大计算量的反演算法，无需其他的遥感产品，在卫星数据接收后短时间内即可合成出来。由于与人眼直观观测具有一致性，业务人员稍加培训即可完成对该类影像的解译。因此，真彩色合成提高了遥感数据的应用效率，特别是在减灾应急领域中，能够及时有效地反映灾区的相关情况，便于开展抢险救灾工作。除此之外，真彩色合成影像对影像信息的还原程度好，也可用于检验相关反演产品的准确性，应用范围较广。
真彩色图像可以真实地反映大气和地球表面特征。在真彩色图像中，云和雪是明亮的白色，霾和烟是半透明的灰色。在晴朗少云的地区，可以看到森林和作物呈现绿色，沙漠为棕色，海洋和湖泊为蓝色。
第一景卫星遥感真彩色影像由美国静止轨道遥感卫星ATS-Ⅲ的多色旋转扫描云遮相机(multicolor spin-scan cloud cover camera，MSSCC)传感器于1967年拍摄[2]。1997年发射的极轨卫星SeaStar/OrbView-2搭载的海洋宽视场扫描仪(sea-viewing wide field-of-view sensor，SeaWiFS)会针对特定任务有选择地生成真彩图产品[3]。Terra/Aqua卫星搭载的中分辨率成像光谱仪(moderate-resolution imaging spectroradiometer，MODIS)传感器将分辨率为500 m的绿波段和蓝波段重采样到250 m后，分别提供上午和下午的250 m业务化的真彩色影像，由WorldView应用向外发布全球拼接后的真彩色影像图，同时由LANCE(land, atmosphere near real-time capability for EOS)提供准实时的单景真彩色影像[4]。2011年发射的Suomi NPP(national polar-orbiting partnership)卫星搭载的可见光红外成像辐射计(visible infrared imaging radiometer suite，VIIRS)的第5、4和3通道分别为红、绿和蓝可见光波段，可以用于真彩色影像的合成，NEO(NASA Earth observation)对外发布其每日全球拼接真彩图，分辨率为750 m[5]。日本气象厅于2014年发射葵花-8(Himawari-8)卫星，其搭载的先进葵花成像仪(advanced himawari imager，AHI)具有红、绿和蓝三个可见光波段，由大气研究合作研究所(the cooperative institute for research in the atmosphere，CIRA)为其开发GeoColor算法合成了真彩色影像并对外发布[6]。美国新一代静止轨道气象卫星GOES-R搭载的先进基准图像仪(advanced baseline imager，ABI)，在可见光波段只有红和蓝两个波段，通过合成绿通道方法实现了真彩色影像合成，并作为实时产品对外发布。ABI有两套真彩色合成方案，1)GeoColor方案：采用16天MODIS地表反射率产品建立红波段(0.64 μm)、近红外波段(0.86 μm)和蓝波段(0.47 μm)查找表，模拟ABI绿波段反射率[7]；2)NatureColor方案：使用近红外、红光和蓝光波段与绿波段反射率的统计关系，构建回归公式，估计ABI绿波段反射率，进而合成真彩色影像[8]。
综上所述，气象遥感卫星和对地观测卫星均把真彩色影像长期作为一项重要的业务化产品。尽管新一代遥感设备通道数量越来越多、分辨率越来越高，但是真彩色影像的直观、快速、易理解的特点使得该项产品仍然是这两类卫星不可或缺的重要产品之一。
风云三号D星(FY-3D)是中国的第二代极轨气象卫星，于2017年底发射升空，每天可以完成对全球区域的覆盖观测，其搭载的MERSI-Ⅱ传感器具有红、绿、蓝三个250 m的可见光通道，满足制作真彩色影像的基本要求。本文以MERSI-Ⅱ数据为例介绍真彩色影像合成的具体算法并展示说明真彩色影像在相关业务领域所具有的重要作用。
1 数据FY-3D是中国的第二代极轨气象卫星，于2017年11月发射升空，为下午轨道卫星，在高度830 km、倾角98.75°的轨道上每天绕地球南北极飞行14圈，每圈大约用时102 min，卫星每天可以完成一次对全球的完整覆盖观测。星上搭载了10台先进的遥感仪器[9]。作为FY-3D星搭载的核心光学仪器，中分辨率光谱成像仪(medium resolution spectral imager Ⅱ，MERSI-Ⅱ)具有6个可见光波段通道，10个可见光/近红外波段通道，3个短波红外波段通道和6个中长波红外波段通道(表 1)，地面分辨率有1 000 m和250 m不等，光谱分辨率有20 nm、50 nm、1.0 μm不等[9]，可以同时获取丰富的地气辐射影像，是目前国际上最先进的宽幅成像遥感仪器之一[10]。该载荷是在风云三号卫星前三颗星配置的两台成像仪器——扫描辐射计和中分辨率光谱成像仪的基础上升级而来的。与上一代的光谱成像仪相比，MERSI-Ⅱ新增了6个红外通道和地面分辨率达250 m的红外分裂窗通道，仪器定标精度和探测灵敏度指标得到了全面提升。MERSI-Ⅱ通过三个250 m可见光通道可以每日无缝隙地获取全球真彩色遥感影像，可以为全球生态环境、灾害监测和气候评估提供观测方案[11]。
|
表 1 MERSI-Ⅱ传感器波段参数 Table 1 Band and parameter of MERSI-Ⅱ |
风云三系列极轨卫星地面数据处理系统将接收到的卫星数据处理成L0、L1和L2等不同级别的数据产品。其中，L0为星上遥感仪器原始测量资料，L1为地理定位和辐射定标后的星上遥感仪器辐射测量资料，L2为经反演算法处理生成的反映大气、云、地表、海表和空间环境变化的各种地球和大气物理参数产品。本算法使用MERSI-Ⅱ L1级数据。MERSI-Ⅱ L1级数据分为250 m和1 000 m两个分辨率的文件，以观测时间5 min为限分块保存为单独的HDF5文件。每个观测数据文件对应两个辅助数据文件，分别为GEOQK和GEO1K。GEOQK文件包含250 m像元经纬度定位数据；GEO1K文件包含1 km像元经纬度定位数据和观测几何数据。真彩色合成所用到的数据文件与数据集如表 2所示。
|
|
表 2 MERSI-Ⅱ L1数据文件与数据集 Table 2 Description of MERSI-Ⅱ L1 datase |
本文所介绍的真彩色影像全球拼图产品制作流程包含三个主要部分：第一部分是快速大气校正，去除由于大气分子散射和吸收对遥感数据造成的影响；第二部分是分段线性拉伸，用于突出弱反射率地物的纹理细节；第三部分是相邻轨道间影像无缝拼接，用于从视觉上实现不同轨道间数据的平滑过渡。
2.1 大气校正大气校正是遥感影像处理中不可或缺的重要环节。可见光波段的大气校正主要用于去除大气分子散射和吸收对遥感影像造成的影响，以获得地表的准确信息，更加客观准确地反映研究区域地物特征。目前大气校正方法主要有:不变目标法(invariable-object methods)[12]、直方图匹配法(histogram matching method)[13]、黑暗目标法(dark object method)[14]、对比度衰减法(contrast reduction methods)[15]，以及基于辐射传输的模型，如：LOWTRAN[16]、MODTRAN[17]、ATREM[18]、6S(second simulation of the satellite signal in the solar spectrum)模型[19]等。其中，6S模型建立在辐射传输理论基础之上，综合了地形、气象、光谱等多种参数，而且适用于多种卫星传感器的不同波段范围，不受研究区特点及目标类型等的影响，在大气辐射和遥感等学科应用广泛[20]。
本文首先基于6S辐射传输模型，对MERSI-Ⅱ 470 nm、550 nm、650 nm三个波段数据进行大气校正。将地表近似为朗伯面进行大气校正。大气校正过程中仅考虑大气分子的贡献，大气气溶胶造成的影响不在本文讨论的范围内。根据6S辐射传输模型的公式[19]，大气校正涉及到的大气散射、大气透过率与地表反射率之间的关系计算可以表达为：
| $ \begin{array}{l} {\rho ^{{\rm{TOA}}}}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}, \varphi } \right) = {T^o}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}} \right)\left[ {{\rho _{\rm{R}}}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}, \varphi } \right) + } \right.\\ \left. {\frac{{{\rho _{\rm{t}}}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}, \varphi } \right)}}{{1 - {\rho _{\rm{t}}}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}, \varphi } \right)S}}{T^H}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}} \right)T{ \downarrow _{\rm{R}}}\left( {{\mu _{\rm{s}}}} \right)T{ \uparrow _{\rm{R}}}\left( {{\mu _{\rm{v}}}} \right)} \right] \end{array} $ | (1) |
式(1)中，μs=cosθs，μv=cosθv，φ=φs-φv，ρTOA(μs, μv, φ)为卫星传感器接收到的信号经过辐射校正与太阳天顶角校正后的大气层顶表观反射率，ρRμs, μv, φ为大气分子散射(Rayleigh scattering)所构成的路径辐射反射率；ρt(μs, μv, φ)为地表反射率；TO(μs, μv)为臭氧吸收造成的大气透过率；TH(μs, μv)为大气水汽透过率；T↓R(μs)和T↑R(μv)分别为大气分子下行辐射透过率和大气分子上行辐射透过率；μs和μv分别为太阳天顶角和观测天顶角的余弦值，φ为太阳与传感器之间的相对方位角；θs、θv、φs和φv分别为太阳天顶角、观测天顶角、太阳方位角和观测方位角；S为大气球形反照率。大气校正就是通过卫星观测参数、大气相关参数以及波段信息计算地表反射率ρt的过程。
2.1.1 大气路径辐射项反射率计算式(1)中大气路径辐射项反射率ρR(μs, μv, φ)需要对单次散射和多次散射求和得到，具体计算公式[21]如下：
| $ \begin{array}{l} {\rho _{\rm{R}}}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}, \varphi } \right) = \sum\limits_{m = 0}^2 {\left( {2 - {\delta _{0, m}}} \right)} \times \rho _l^m\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}, \tau } \right) \times \\ \cos [m(\varphi )] + \left( {1 - {e^{ - \tau /{\mu _{\rm{s}}}}}} \right) \times \left( {1 - {e^{ - \tau /{\mu _{\rm{v}}}}}} \right) \times \\ \sum\limits_{m = 0}^2 {\left( {2 - {\delta _{0, m}}} \right)} \times {\mathit{\Delta }^m}(\tau ) \times {P^m}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}} \right) \times \cos [m(\varphi )] \end{array} $ | (2) |
式(2)中，τ为目标地物大气分子散射在对应波段的光学厚度(表 3)，该值与目标区域的大气压强(p)和标准大气压(p0)下大气分子在该波段的光学厚度(τ0)有关，τ0的计算采用6S辐射传输模式，通过输入MERSI传感器对应波段的光谱响应函数计算得到；δ0, m为克罗内克函数；Pm(μs, μv)为大气分子散射相函数经过傅里叶展开后第m项数值，前三项计算公式如式(3)；ρlm(μs, μv, τ)为第m项散射相函数对应的单次散射反射率，计算公式如式(4)；Δm(τ)为大气分子光学厚度对应多次散射的订正项，前三次项的计算公式如式(5)。
|
|
表 3 风云三D星MERSI-Ⅱ传感器前三波段臭氧吸收系数 Table 3 Ozone absorption coefficient in the first three bands of FY-3D MERSI-Ⅱ |
| $ \begin{array}{l} {P^0}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}} \right) = 1 + \left( {3\mu _{\rm{s}}^{\rm{2}} - 1} \right)\left( {3\mu _{\rm{v}}^2 - 1} \right) \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{1 - \frac{\delta }{{2 - \delta }}}}{{1 + 2\frac{\delta }{{2 - \delta }}}} \times \frac{1}{8}\\ {P^1}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}} \right) = - {\mu _{\rm{s}}}{\mu _{\rm{v}}}\left( {1 - \mu _{\rm{s}}^{\rm{2}}} \right)\frac{1}{2}\left( {1 - \mu _{\rm{v}}^{\rm{2}}} \right)\frac{1}{2} \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{{1 - \frac{\delta }{{2 - \delta }}}}{{1 + 2\frac{\delta }{{2 - \delta }}}} \times \beta \times 1.5\\ {P^2}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}} \right) = \left( {1 - \mu _{\rm{s}}^2} \right)\left( {1 - \mu _{\rm{v}}^2} \right) \times \\ \;\;\;\;\;\;\;\;\frac{{1 - \frac{\delta }{{2 - \delta }}}}{{1 + 2\frac{\delta }{{2 - \delta }}}} \times \beta \times 0.375 \end{array} $ | (3) |
其中，δ=0.027 9，β=0.5。
| $ \begin{array}{l} \rho _l^{\rm{m}}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}, \tau } \right) = {P^{\rm{m}}}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}} \right) \times \\ \left( {1 - {e^{ - \tau \left( {\frac{1}{{{\mu _{\rm{s}}}}} + \frac{1}{{{\mu _{\rm{v}}}}}} \right)}}} \right) \times \frac{1}{{4\left( {{\mu _{\rm{s}}} + {\mu _{\rm{v}}}} \right)}} \end{array} $ | (4) |
| $ {\mathit{\Delta }^{\rm{m}}}(\tau ) = {a^{\rm{m}}} + {b^{\rm{m}}} \times \ln (\tau ) $ | (5) |
式(5)中，当m取0时，a0和b0的计算公式如下：
| $ \begin{array}{l} {a^0} = a_0^0 + a_1^0{\mu _{\rm{s}}}{\mu _{\rm{v}}} + a_2^0{\left( {{\mu _{\rm{s}}}{\mu _{\rm{v}}}} \right)^2} + \\ \;\;\;\;\;a_3^0\left( {{\mu _{\rm{s}}} + {\mu _{\rm{v}}}} \right) + a_4^0\left( {\mu _{\rm{s}}^2 + \mu _{\rm{v}}^2} \right) \end{array} $ | (6) |
| $ \begin{array}{l} {b^0} = b_0^0 + b_1^0{\mu _{\rm{s}}}{\mu _{\rm{v}}} + b_2^0{\left( {{\mu _{\rm{s}}}{\mu _{\rm{v}}}} \right)^2} + \\ \;\;\;\;b_3^0\left( {{\mu _{\rm{s}}} + {\mu _{\rm{v}}}} \right) + b_4^0\left( {\mu _{\rm{s}}^2 + \mu _{\rm{v}}^2} \right) \end{array} $ | (7) |
式(6)中，a0-40的取值分别为0.332 438、-0.103 244、0.162 854、0.114 933、-0.309 248；式(7)中，b0-40的取值分别为-0.067 771、0.032 417、0.001 577、-0.035 037、-0.012 409。式(5)中，当m分别取1和2时，对应a1、a2、b1和b2分别取值0.196 66、0.145 459、-0.054 391和-0.029 108。
2.1.2 大气透过率计算真彩色影像合成使用了MERSI传感器前三个波段，不考虑大气气溶胶的影响，这三个波段的大气透过率可以近似认为由臭氧透过率TO(μs, μv)、大气水汽透过率TH(μs, μv)和大气分子透过率(包括大气下行辐射透过率T↓R(μs)和上行辐射透过率T↑R(μv))组成的。臭氧透过率TO(μs, μv)计算公式如下：
| $ {T^o}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}} \right) = \exp \left[ { - \left( {\frac{1}{{{\mu _{\rm{s}}}}} + \frac{1}{{{\mu _{\rm{v}}}}}} \right) \times U \times A} \right] $ | (8) |
式(8)中U为大气层臭氧含量，本文取0.319 cm；A为臭氧吸收系数，可以通过KNEIZYS et al.[16]的文献中检索得到，查到的结果如表 3所示。
大气水汽透过率TH(μs, μv)的计算公式如下[14]：
| $ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{T^H}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}} \right) = \\ \exp \left\{ { - \exp \left[ {{A^H} + {B^H}\ln \left( {\frac{1}{{{\mu _{\rm{s}}}}}{U^H} + \frac{1}{{{\mu _{\rm{v}}}}}{U^H}} \right)} \right]} \right\} \end{array} $ | (9) |
式(9)中UH为大气水汽含量，这里取2.93 g·cm-2；AH和BH为水汽吸收系数，可以通过表 4查得。
|
|
表 4 风云三D星MERSI-Ⅱ传感器前三波段水汽吸收系数 Table 4 Water vapor absorption coefficient in the first three bands of FY-3D MERSI-Ⅱ |
大气分子散射没有吸收且不对称因子为0，大气下行辐射透过率T↓R(μs)和上行辐射透过率T↑R(μv)可以通过下式求得：
| $ T{ \downarrow _{\rm{R}}}\left( {{\mu _{\rm{s}}}} \right) = \frac{{\left( {\frac{2}{3} + {\mu _{\rm{s}}}} \right) - \left( {\frac{2}{3} - {\mu _{\rm{s}}}} \right){e^{ - \tau /{\mu _{\rm{s}}}}}}}{{\frac{4}{3} + {\mu _{\rm{s}}}}} $ | (10) |
| $ T{ \uparrow _{\rm{R}}}\left( {{\mu _{\rm{v}}}} \right) = \frac{{\left( {\frac{2}{3} + {\mu _{\rm{v}}}} \right) - \left( {\frac{2}{3} - {\mu _{\rm{v}}}} \right){e^{ - \tau /{\mu _{\rm{v}}}}}}}{{\frac{4}{3} + {\mu _{\rm{v}}}}} $ | (11) |
针对完全散射且无吸收的分子散射，其球形反照率S的计算公式如下：
| $ S = \frac{1}{{4 + 3\tau }}\left[ {3\tau - 4{E_3}(\tau ) + 6{E_4}(\tau )} \right] $ | (12) |
式中En计算公式如下：
| $ {E_n}(\tau ) = \mathop \smallint \limits_1^\infty \frac{{{e^{ - \tau t}}}}{{{t^{\rm{n}}}}}{\rm{d}}t $ | (13) |
式中τ为大气分子散射的光学厚度。
2.1.4 大气校正计算大气校正处理是逐像元进行的，其流程如图 1所示。读取MERSI-Ⅱ L1数据470 nm、550 nm、650 nm三个波段的像素值，根据波段的定标系数计算成反射率数据，将此反射率数据除以太阳天顶角的余弦值μs得到大气层顶表观反射率ρTOA(μs, μv, φ)。读取GEOQK文件中对应位置像元的太阳天顶角θs、太阳方位角φs、观测天顶角θv、观测方位角φv和数字高程模型(DEM)数据。根据DEM数据可以计算出目标像元修正后的大气分子光学厚度τ。结合大气分子光学厚度数据、太阳与观测角度数据计算出目标像元处的大气路径辐射项反射率ρR(μs, μv, φ)、大气透过率(包括臭氧透过率TO(μs, μv)、大气水汽透过率TH(μs, μv)、大气下行辐射透过率T↓R(μs)和大气上行辐射透过率T↑R(μv))和大气球形反照率S，将上述数据带入式(14)、(15)计算出像元的地表反射率ρt(μs, μv, φ)：
|
图 1 大气校正流程图 Fig.1 Flow chart of atmospheric correction |
| $ \begin{array}{l} {\rho _{\rm{t}}}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}, \varphi } \right):\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rho _{\rm{t}}}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}, \varphi } \right) = t{\rm{/}}(1 + tS) \end{array} $ | (14) |
| $ \begin{array}{l} t = \left( {{\rho ^{{\mathop{\rm TOA}\nolimits} }}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}, \varphi } \right)/{T^O}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}} \right) - {\rho _{\rm{R}}}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}, \varphi } \right)} \right)/\\ \left. {\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {{T^H}\left( {{\mu _{\rm{s}}}, {\mu _{\rm{v}}}} \right)T{ \downarrow _{\rm{R}}}} \right.\left( {{\mu _{\rm{s}}}} \right)T{ \uparrow _{\rm{R}}}\left( {{\mu _{\rm{v}}}} \right)} \right) \end{array} $ | (15) |
2.1节中计算得到的反射率数据浮点型数据，值的范围在0.0~1.0之间。为了合成RGB影像，需要将反射率值转换为0~255的整数。我们将反射率数据采用线性拉伸的方式转换为整数，转换公式如下：
| $ y = {\mathop{\rm int}} \left( {255 \times {\rho _t}} \right) $ | (16) |
式中y表示输出的8位整型数值；int表示取整函数；ρt为大气校正得到的地表反射率。采用650 nm、550 nm和470 nm三个可见光波段线性拉伸后的数据，对应R、G、B三个通道直接合成为真彩色影像(图 2a)，但是对于暗目标地物，例如浓密植被区和水体，上述方法合成的影像难以体现其具体的纹理和灰度信息。为了突出暗目标地物的特征，本方法使用GUMLEY et al.[4]开发的非线性亮度增强表(表 5)对转换后每个波段的整型数据进行进一步非线性拉伸，然后再进行真彩色合成。非线性增强后合成的真彩色图像如图 2b所示。可以看出，与图 2a相比，图 2b中较暗的地物亮度有所增加。
|
图 2 线性拉伸(a)与非线性亮度增强(b)真彩色合成影像对比 Fig.2 Comparison of true color composite image by linear stretch (a) and nonlinear brightness enhancement (b) |
日常业务应用中，由于极轨卫星的轨道特征限制，单景五分钟影像很难完全覆盖特定任务区域，为了覆盖完整的任务区域需要将多景真彩色影像拼接到一起以供使用。本小节涉及到的真彩色影像除非特殊说明，均为经过初步几何校正，坐标系为等经纬度坐标系的影像。
将多景真彩色影像的拼接分为以下两类：1)邻近两个五分钟段的影像拼接，称为“同轨拼接”；2)时间上不是相邻时次但空间上有重叠区域的两个影像拼接，称为“不同轨拼接”。
2.3.1 同轨拼接同轨数据拼接涉及的两景或者多景影像没有重叠区域，因此只需要通过逐像元的经纬度信息即可将多景影像无缝地拼接到一起。由于同轨的多景影像处于同一个连续变化的光照条件，因此拼接后的影像过渡自然没有任何拼接痕迹(图 3)。
|
图 3 2018年12月5日同轨相邻两景影像拼接效果(a. 05:20 UTC影像；b. 05:15 UTC影像；c.拼接结果) Fig.3 Mosaic effect of two adjacent images in the same orbit on 5 December 2018 (a. image at 05:20 UTC; b. image at 05:15 UTC; c. the mosaic result) |
相邻不同轨的影像存在重叠区域，且光照条件并不连续。以FY-3D为例，同一纬度相邻两轨数据的时间间隔在100 min左右。由于地表存在BRDF特性以及光照条件并不连续的原因，直接使用像素经纬度信息拼接后的数据会存在明显的不连续过渡特征，另一方面由于存在重叠区域，如何选择镶嵌线的位置也是需要考虑的问题。
为了解决不同轨图像拼接出现的明显间断现象，在视觉上尽可能实现不同轨数据间的平滑过渡，我们提出了距离加权图像融合拼接方法(distance weighted pixel blending, DWPB)。具体原理如下，在两景相邻的不同轨影像重叠区，以拼接线为中心向各自图像两侧做加权的RGB像素值合成。RGB像素值的合成公式如下：
| $ \left\{ \begin{array}{l} DN = D{N_{{\rm{left}}}}{W_{{\rm{left}}}} + D{N_{{\rm{right}}}}{W_{{\rm{right}}}}\\ \;\;\;\;\;\;\;\;\;\;\;{W_{{\rm{left}}}} + {W_{{\rm{right}}}} = 1 \end{array} \right. $ | (17) |
式中DN为两景影像重叠区合成后的像素值；DNleft为左侧景影像的像素值；DNright为右侧景影像的像素值；Wleft和Wright为左侧和右侧景影像像素值的权重，取值范围在0.0~1.0之间，权重值的计算方法如下：
| $ \left\{ \begin{array}{l} {W_{{\rm{ left }}}} = \frac{{\left( {x - {x_0}} \right)}}{{2\left( {{x_1} - {x_0}} \right)}} + 0.5, \left( {x < {x_0}} \right)\\ {W_{{\rm{ right }}}} = \frac{{\left( {x - {x_0}} \right)}}{{2\left( {{x_{\rm{r}}} - {x_0}} \right)}} + 0.5, \left( {x \ge {x_0}} \right) \end{array} \right. $ | (18) |
式中x0为拼接线的x坐标值；xl和xr分别为拼接区域最左端x坐标值和最右端x坐标值。x0的确定可以通过计算重叠区相同y坐标值对应的最左和最右端x值的中间值得到；xl和xr分别为x0向左和向右平移的距离。这两个距离如果太小会造成拼接线处过渡生硬，如果太大会造成重叠区域比较模糊。经过多次试验，我们选定xl和xr的取值是通过计算重叠区x方向宽度除以2得到，但xl和xr的最大值不超过200个像素(图 4)。
|
图 4 DWPB方法示意图 Fig.4 Illustration of DWPB method |
|
图 5 台风区域DWPB拼接(a)和普通拼接(b)结果对比 Fig.5 Comparison of DWPB method (a) and normal method (b) in typhoon area |
|
图 6 陆地区域采用DWPB拼接(a)和普通拼接(b)结果对比 Fig.6 Comparison of DWPB method (a) and normal method (b) in land area |
采用上述真彩色合成方法，利用MERSI-Ⅱ数据，分别针对2018年的几大天气和环境变化现象，获取合成其真彩色图并加以分析和解释。
3.1 东欧沙尘与积雪2018年3月下旬，东欧与俄罗斯接壤的地区地面积雪有明显的橙色。这种颜色来自撒哈拉沙漠的大量沙尘，这些沙尘被强风吹到地中海上空，堆积在保加利亚、罗马尼亚、摩尔多瓦、乌克兰和俄罗斯地区。FY-3D卫星于2018年3月22日获得了一景地中海区域的真彩图(图 7a)，这幅真彩色影像清晰地还原了来自北非的尘埃穿越地中海到达东欧的过程。22日早些时候的一景真彩色影像显示出东欧与俄罗斯接壤地区的积雪尚未受到沙尘的污染，仍然呈现雪白色(图 7b)。几天以后FY-3D卫星于2018年3月25日获得了该地区的真彩色影像，该景影像显示出由于沙尘的沉降导致该地区地表积雪呈现出橘黄色(图 7c)。
|
图 7 2018年3月22日地中海上空沙尘(a)、东欧地区白色积雪(b)和2018年3月25日东欧地区橘黄色积雪(c) Fig.7 Sand and dust transport over Mediterranean Sea (a) on 22 March 2018, white snow cover in Eastern Europe (b) on 22 March 2018, and orange snow cover in Eastern Europe (c) on 25 March 2018 |
2018年6月2日，新疆出现沙尘天气。12:20(北京时)的FY-3D气象卫星监测(图 8)显示，沙尘区主要位于新疆南疆盆地西部，部分沙尘区上空有云覆盖。据估算，卫星可视的沙尘区面积约为19万km2。受沙尘天气影响，上述地区空气质量、能见度下降。
|
图 8 新疆地区沙尘暴 Fig.8 Sand and dust storm in Xinjiang region |
2018年5月6日，FY-3D的MERSI-Ⅱ获得了北海(North Sea)浮游植物大量爆发的真彩色图像(图 9)。由图 9可以看到，浮游植物正在向西移动。水体呈现淡乳白色的区域可能含有颗石藻类(Coccolithophores)，呈现绿色的水域可能含有硅藻类(Diatoms)。除此之外，真彩色图的颜色的亮度可以反映浮游植物的密度，而不同的漩涡和形状可以反映出洋流、涡流和潮汐的复杂运动。北海全年在春末和夏初海水营养含量最高，浮游植物数量最多，充足的阳光、上升的水温和丰富的营养物质，引发浮游植物大量繁殖[22]。
|
图 9 大量繁殖的北海浮游植物 Fig.9 Phytoplankton bloom in the North Sea |
FY-3D气象卫星2018年7月28日21:45(UTC)资料监测显示，美国加利福尼亚州、爱达荷州、内华达州、犹他州等四个州出现11处较大火点，集中分布于5个区域，影响范围约298.8 km2，经估算明火区面积约127.7 km2(图 10)。
|
图 10 美国西部大火MERSI-Ⅱ真彩色影像(红色区域为着火点) Fig.10 MERSI-Ⅱ true color image of big fires in the west of the United States of America (red area for active fire area) |
FY-3D气象卫星于2018年9月16日05:40(UTC)“捕捉”到1822号台风“山竹”(强台风级)接近我国东部沿海的过程。据中国气象局国家卫星中心9月16日发布的2018年第64期《气象卫星台风监测报告》[23]，10:00，台风山竹的中心位于(20.8°N，115.0°E)，中心附近最大风速48 m·s-1，中心最低气压945 hPa。FY-3D监测图像(图 11)显示：“山竹”台风整体结构完整，台风眼依稀可见，外围螺旋云带已经开始影响广东大部。
|
图 11 台风山竹2018年9月16日MERSI-Ⅱ真彩色影像 Fig.11 MERSI-Ⅱ true color image of Typhoon MANGKHUT on 16 September 2018 |
本文介绍了针对FY-3D卫星MERSI-Ⅱ数据的真彩色影像合成算法，其中发展了基于6S模式的大气校正算法、非线性加强算法与加权融合拼接算法。在大气校正算法中，包括臭氧和水汽透过率，上行和下行辐射透过率以及球形反照率的校正。非线性增强旨在区分暗目标的纹理和详细特征。加权融合拼接算法能够以在相邻轨道和不同扫描线之间创建平滑过渡的真彩色图像。应用合成后的真彩色影像，分析了在不同地表和天气现象如：沙尘、积雪、蓝藻、火点烟羽和台风等，真彩色影像起到的重要作用。样例分析结果表明：真彩图提高遥感影像的对比度，还原了地表的真实色彩；在沙尘、蓝藻爆发、火点识别和台风等方面可以提供便于识别和判读的数据资料。在以后的研究工作中，要利用真彩色图技术开发出高质量的业务产品，进而更好地服务于社会大众。
| [1] |
日本遥感研究会.遥感精解(修订版)[M].刘勇卫, 译.2版.北京: 测绘出版社, 2010: 192.
|
| [2] |
GREAVES J R, SHENK W E.The development of the geosynchronous weather satellite system[M]//SCHNAPF A.Progress in Astronautics and Aeronautics Volume 97: Monitoring earth's ocean, land, and atmosphere from space-sensors, systems, and applications.Reston, Virginia: American Institute of Aeronautics and Astronautics, 1985: 150-181. DOI: 10.2514/5.9781600865725.0150.0181.
|
| [3] |
HOOKER S B, ESAIAS W E, FELDMAN G C, et al. NASA Technical Memorandum 104566, Volume1: An Overview of SeaWiFS and ocean color[R]//HOOKER S B, FIRESTONE E R. SeaWiFS Technical Report Series. Greenbelt, Maryland: NASA Goddard Space Flight Center, 1992: 1-24.
|
| [4] |
GUMLEY L, DESCLOITRES J, SCHMALTZ J. Creating reprojected true color MODIS images: A tutorial[R/OL].(2010-01-14)[2019-03-16].Madison: Space Science and Engineering Center, University of Wisconsin-Madison. http://gis-lab.info/docs/modis_true_color.pdf.
|
| [5] |
HILLGER D, SEAMAN C, LIANG C, et al. Suomi NPP VIIRS imagery evaluation[J]. J Geophys Res: Atmos, 2014, 119(11): 6440-6455. DOI:10.1002/2013JD021170 |
| [6] |
BIKOS D, MILLER S, LINDSEY D. GeoColor product: Quick guide[R/OL]. (2017-10-19)[2019-03-16].Washington, DC: NOAA, 2017. http://rammb.cira.colostate.edu/research/goes-r/proving_ground/cira_product_list/docs/QuickGuide_CIRA_Geocolor_20171019.pdf.
|
| [7] |
GROSSBERG M D, SHAHRIAR F, GLADKOVA I, et al. Estimating true color imagery for GOES-R[C]//SHEN S S, LEWIS P E. Proceedings of SPIE, Volume 8048: Algorithms and technologies for multispectral, hyperspectral, and ultraspectral imagery XVII. Bellingham, Washington: The International Society for Optical Engineering, 2011: 80481A. DOI: 10.1117/12.884020.
|
| [8] |
MILLER S D, SCHMIT T L, SEAMAN C J, et al. A sight for sore eyes: The return of true color to geostationary satellites[J]. Bull Amer Meteor Soc, 2016, 97(10): 1803-1816. DOI:10.1175/BAMS-D-15-00154.1 |
| [9] |
朱爱军, 胡秀清, 贾树泽, 等. 风云三号D气象卫星全球数据获取方法及数据分发[J]. 海洋气象学报, 2018, 38(3): 1-10. |
| [10] |
中国科学院.风云三号D星升空技物所贡献两大光学主载荷[EB/OL]. (2017-11-16)[2019-03-16]. http://www.sitp.cas.cn/xwzx/kydt/201711/t20171116_4895184.html.
|
| [11] |
中国科学院.中科院七大载荷助力风云三号D星观测"风云"[EB/OL]. (2017-11-16)[2019-03-16].http://www.cas.cn/syky/201711/t20171117_4622313.shtml.
|
| [12] |
ZHENG W, LIU C, ZENG Z Y, et al. A feasible atmospheric correction method to TM image[J]. J China U Min Techno, 2007, 17(1): 112-115. DOI:10.1016/S1006-1266(07)60024-8 |
| [13] |
RICHTER R. A spatially adaptive fast atmospheric correction algorithm[J]. Int J Remote Sens, 1996, 17(6): 1201-1214. DOI:10.1080/01431169608949077 |
| [14] |
TEILLET P M, FEDOSEJEVS G. On the dark target approach to atmospheric correction of remotely sensed data[J]. Can J Remote Sens, 1995, 21(4): 374-387. DOI:10.1080/07038992.1995.10855161 |
| [15] |
TANRE D, LEGRAND M. On the satellite retrieval of Saharan dust optical thickness over land: Two different approaches[M]. J Geophys Res: Atmos, 1991: 5221-5227.
|
| [16] |
KNEIZYS F X, SHETTLE E, GALLERY W, et al. Atmospheric transmittance/radiance: Computer code LOWTRAN 5[M]. Bedford, Massachusetts: Air Force Geophysics Laboratory, 1980.
|
| [17] |
XU Y L, WANG R S, LIU S W, et al. Atmospheric correction of hyperspectral data using MODTRAN model[C]//TONG Q X. Proceedings of SPIE, Volume 7123: Remote Sensing of the Environment: 16th National Symposium on Remote Sensing of China. Bellingham, Washington: Society of Photo-optical Instrumentation Engineers, 2008: 712306. DOI: 10.1117/12.815552.
|
| [18] |
Gao B C, Montes M J, Ahmad Z, et al. Atmospheric correction algorithm for hyperspectral remote sensing of ocean color from space[J]. Appl Opt, 2000, 39(6): 887-896. DOI:10.1364/AO.39.000887 |
| [19] |
VERMOTE E F, TANRé D, DEUZé J L, et al. Second simulation of the satellite signal in the solar spectrum, 6S: An overview[J]. IEEE Trans Geosci Remote Sens, 1997, 35(3): 675-686. DOI:10.1109/36.581987 |
| [20] |
武永利, 栾青, 田国珍. 基于6S模型的FY-3A/MERSI可见光到近红外波段大气校正[J]. 应用生态学报, 2011, 22(6): 1537-1542. |
| [21] |
NOAA NESDIS Center for Satellite Applications and Research. GOES-R advanced baseline imager (ABI) algorithm theoretical basis document for suspended matter/aerosol optical depth and aerosol size parameter Version3.0[R]. Washington, DC: NOAA, 2012.
|
| [22] |
CAPUZZO E, LYNAM C P, BARRY J, et al. A decline in primary production in the North Sea over 25 years, associated with reductions in zooplankton abundance and fish stock recruitment[J]. Glob Change Biol, 2018, 24(1): e352-e364. DOI:10.1111/gcb.2018.24.issue-1 |
| [23] |
国家卫星气象中心. 2018年第64期气象卫星台风监测报告[R].北京: 国家卫星气象中心, 2018.
|