CN112598608B - 一种基于目标区域的光学卫星快速融合产品制作方法 - Google Patents
一种基于目标区域的光学卫星快速融合产品制作方法 Download PDFInfo
- Publication number
- CN112598608B CN112598608B CN202011340149.6A CN202011340149A CN112598608B CN 112598608 B CN112598608 B CN 112598608B CN 202011340149 A CN202011340149 A CN 202011340149A CN 112598608 B CN112598608 B CN 112598608B
- Authority
- CN
- China
- Prior art keywords
- image
- imaging
- roi
- ccd
- virtual
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 230000004927 fusion Effects 0.000 title claims abstract description 54
- 238000000034 method Methods 0.000 title claims abstract description 47
- 238000004519 manufacturing process Methods 0.000 title claims abstract description 19
- 230000003287 optical effect Effects 0.000 title claims abstract description 19
- 238000003384 imaging method Methods 0.000 claims abstract description 125
- 238000012545 processing Methods 0.000 claims abstract description 32
- 238000007499 fusion processing Methods 0.000 claims abstract description 10
- 238000012937 correction Methods 0.000 claims abstract description 8
- 239000011159 matrix material Substances 0.000 claims description 15
- 239000000523 sample Substances 0.000 claims description 14
- 238000004364 calculation method Methods 0.000 claims description 12
- 230000008569 process Effects 0.000 claims description 12
- 239000003795 chemical substances by application Substances 0.000 claims description 9
- 230000002708 enhancing effect Effects 0.000 claims description 8
- 238000001914 filtration Methods 0.000 claims description 6
- 230000000007 visual effect Effects 0.000 claims description 6
- 230000004075 alteration Effects 0.000 claims description 5
- 238000012892 rational function Methods 0.000 claims description 5
- 230000003595 spectral effect Effects 0.000 claims description 5
- 230000008859 change Effects 0.000 claims description 3
- 238000013461 design Methods 0.000 claims description 3
- 238000009826 distribution Methods 0.000 claims description 3
- 230000000694 effects Effects 0.000 claims description 3
- 230000002401 inhibitory effect Effects 0.000 claims description 3
- 238000013507 mapping Methods 0.000 claims description 3
- 238000005457 optimization Methods 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 230000017105 transposition Effects 0.000 claims description 3
- 238000004458 analytical method Methods 0.000 claims description 2
- 230000015572 biosynthetic process Effects 0.000 claims description 2
- 238000010276 construction Methods 0.000 claims description 2
- 238000001228 spectrum Methods 0.000 claims description 2
- 238000001444 catalytic combustion detection Methods 0.000 claims 11
- 230000005764 inhibitory process Effects 0.000 claims 1
- 238000005516 engineering process Methods 0.000 abstract description 4
- 238000010191 image analysis Methods 0.000 abstract description 2
- 230000004044 response Effects 0.000 abstract description 2
- 230000005540 biological transmission Effects 0.000 description 3
- 238000011161 development Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000000701 chemical imaging Methods 0.000 description 1
- 238000005520 cutting process Methods 0.000 description 1
- 238000013500 data storage Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000010354 integration Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000008520 organization Effects 0.000 description 1
- 230000036544 posture Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000011002 quantification Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000003860 storage Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000001629 suppression Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J3/00—Spectrometry; Spectrophotometry; Monochromators; Measuring colours
- G01J3/28—Investigating the spectrum
- G01J3/2823—Imaging spectrometer
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/40—Image enhancement or restoration using histogram techniques
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T7/00—Image analysis
- G06T7/70—Determining position or orientation of objects or cameras
- G06T7/73—Determining position or orientation of objects or cameras using feature-based methods
- G06T7/75—Determining position or orientation of objects or cameras using feature-based methods involving models
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01J—MEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
- G01J3/00—Spectrometry; Spectrophotometry; Monochromators; Measuring colours
- G01J3/28—Investigating the spectrum
- G01J3/2823—Imaging spectrometer
- G01J2003/2826—Multispectral imaging, e.g. filter imaging
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
- G06T2207/10036—Multispectral image; Hyperspectral image
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
- G06T2207/10041—Panchromatic image
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Spectroscopy & Molecular Physics (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种基于目标区域的光学卫星快速融合产品制作方法。首采用一种动态搜索方式确定目标兴趣区域（ROI）的成像时间窗口；将ROI四角点的经纬度坐标反投至每一片CCD影像上，从而确定ROI成像影像范围；然后对ROI对应的原始数据范围进行影像解析处理，并对ROI原始全色和多光谱影像进行全谱段一体化虚拟重成像处理，输出面向融合的全色多光谱同分辨率影像产品；之后采用一种高色彩保真的全色多光谱融合技术对ROI全色多光谱进行融合处理；最后对ROI融合影像进行间接几何纠正处理，生成带地理坐标的目标融合产品。本发明具有数据量小、计算量小、任务响应迅速的优势，提升面向应急任务的信息获取效率。
Description
技术领域
本发明属于遥感卫星影像处理领域，涉及一种基于目标区域的光学卫星快速融合产品制作方案。
背景技术
近年来，遥感对地观测技术飞速发展，已步入数据获取手段多样化和信息海量化的新时代。我国高分辨率对地观测系统重大专项现已全面启动，到2020年，我国将新研制数十颗高分辨率遥感卫星，实现30～50颗遥感卫星同时在轨工作，保守估计每天获取的数据量将达到数百个TB。为了不断提升遥感卫星应用效能，高分辨率、宽覆盖、多光谱成像是当今和未来遥感对地观测发展的主要目标之一，其显著特点是获取的数据量呈几何级数增长。
然而，现有地面运行系统处理模式和卫星数据获取能力不匹配，严重制约了卫星数据的处理和应用保障效率。以传统的目标融合产品生产服务为例，现有运行系统仍按标准景方式组织产品服务，即：检索到目标所在景后，分别生产全色和多光谱标准景1级产品，并依次进行全色和多光谱1级产品的定向处理、精配准处理和融合处理，再对标准景融合产品进行几何纠正和切片，从而获取目标区域融合产品。随着卫星空间分辨率的不断提高，标准景产品数据量通常可以达到几十GB，其大侧摆大俯仰条件下的融合产品甚至达到上百GB，受计算资源、存储设备IO、网络带宽、终端单机性能等方面的限制，数据量过大所带来的问题日益凸显，从标准景产品订单下达、标准景融合产品生产、裁切出感兴趣目标区域甚至需要几个小时，已经严重影响了高分辨率卫星日常应用和保障的时效性。因此，必须针对现有地面处理系统中目标融合产品核心处理环节数据组织冗余、处理模型灵活性较差等问题，开展基于目标区域的融合产品快速制作方法的研究，提升面向应急任务的信息获取效率。
发明内容
本发明的目的在于提供一种基于目标区域的光学卫星快速融合产品制作方法，以解决上述背景技术中提出的问题。
为实现上述目的本发明采用以下技术方案：
一种基于目标区域的光学卫星快速融合产品制作方法，包括以下步骤：
步骤1，解析卫星下传原始数据中的辅助数据；
步骤2，确定ROI成像时间窗口；
步骤3，确定ROI成像影像范围；
步骤4，ROI影像解析处理；
步骤5，ROI全谱段一体化虚拟重成像；
步骤6，ROI全色多光谱融合处理；
步骤7，目标融合影像的几何纠正处理。
作为本发明进一步的方案：步骤1所述解析卫星下传原始数据中的辅助数据包括成像时间数据、姿态数据、轨道数据，为ROI目标融合产品快速处理提供几何参数。
作为本发明进一步的方案：步骤2所述确定R0I成像时间窗口的方法如下：
2.1在数据解析过程中，从T0时刻开始，以一定时间间隔Δt0，对当前成像区域位置进行实时计算，即计算推扫成像行两端的经纬度坐标，判断ROI目标点经纬度是否落于计算的经纬度范围内，从而获取ROI的概略时间窗口；步骤如下：
2.11)建立当前成像行的成像时间T0的严格几何成像模型，根据当前成像行的成像时间T0确定轨道参数和姿态参数，再联合相机参数即可建立(相机参数通常会在发射前进行严格的实验室检校，但卫星发射和运行空间环境极其复杂，会导致星上载荷结构发生变化，需进行相机参数的在轨几何定标，才能保证ROI区域的精准定位，此处用到的相机参数即为几何定标后的精确参数)，如下式：
其中
式中(Xg,Yg,Zg)与(Xgps,Ygps,Zgps)分别表示像点对应的物方点及GPS(全球定位系统)天线相位中心在WGS84坐标系下的坐标，后者由卫星上搭载的GPS获取，按成像时间内插得到的轨道参数；
分别代表WGS84坐标系到J2000坐标系的旋转矩阵、J2000坐标系到卫星本体坐标系的旋转矩阵、卫星本体坐标系到相机坐标系的旋转矩阵，
由星敏、陀螺通过组合定姿得到，
由相机参数文件得到；(BX,BY,BZ)body代表从传感器投影中心到GPS天线相位中心的偏心矢量在卫星本体坐标系下的坐标，则在卫星发射前由实验室检校；λ为比例系数；(ψx(s),ψy(s))代表探元s在相机坐标系下的指向角，s代表CCD中的探元号，axi，ayi(i＝0,1,2,3)为多项式系数(相机参数文件中的内定标参数)；
2.12)计算当前行首尾端点的地理坐标，得到点p0，q0；
2.13)T1＝T0+Δt0时刻，重复上述2.11)和2.12)步骤，得到点p1，q1；
2.14)判断ROI目标点经纬度是否落于矩形(p0,q0,p1,q1)内，如不是，则Δt0时候后继续重复上述计算；如是，则停止计算，得到ROI的概略时间窗口(Ti,Ti+1)；
2.2在2.1中确定的概略时间窗口(Ti,Ti+1)内，以Ti时刻为起点，进一步缩小时间间隔为Δt1，重复2.1中步骤，并计算推扫成像行两端的经纬度坐标，判断目标点经纬度是否落于两次计算经纬度范围内，从而获取更加精确的ROI时间窗口；
2.3如此迭代，采用更小的时间间隔Δti，直至获取的时间窗口能完全准确的覆盖ROI的成像范围，为了提高搜索的可靠性，最小时间间隔要能覆盖ROI的成像范围。
作为本发明进一步的方案：步骤3所述确定R0I成像影像范围方法如下：
3.1依据步骤2得到的ROI成像时间窗口即可确定轨道参数和姿态参数，再联合相机参数即可按式(1)建立在该时间窗口内成像的每一片CCD影像的严格几何成像模型；
3.2利用严格几何成像模型计算得到虚拟控制点，通过地形独立法将严格几何成像模型转化为有理函数模型RFM(有理函数模型)，以解决严格几何成像模型坐标反算难以实现的问题；
3.3利用RFM将ROI四角点的经纬度坐标反投至每一片CCD影像空间，得到ROI在图像中的坐标，判断ROI四角点坐标是否在影像成像范围内，从而确定ROI的成像CCD片号和具体的成像起始行列号。
作为本发明进一步的方案：步骤5所述ROI全谱段一体化虚拟重成像方法如下：
5.1无畸变虚拟CCD的设计，高分辨率光学遥感卫星为了获取较大的幅宽，全色和多光谱相机均会在焦平面上安置多个CCD，因此，在重成像过程中，采用一种虚拟的CCD代替原有的多片CCD进行成像，虚拟CCD是一个覆盖整个成像视场的无畸变的理想CCD线阵，且虚拟CCD的主点、主距及探元尺寸均与原始CCD一致，其像素个数是所有分片CCD像素个数总和去除重叠像素个数，同时，为了避免高程误差带来的重成像误差，虚拟CCD的位置通常应尽可能与真实CCD位置接近，保证原始成像视角与虚拟成像视角的差异最小；
5.2虚拟CCD成像模型构建，步骤如下：
5.21)建立虚拟CCD的稳态重成像模型，与原始分片CCD的严格几何成像模型基本相同，不同之处在于该模型是在稳定成像状态下建立，即虚拟扫描行的姿态角则是在滤波处理后的姿态角中内插获取；
5.22)生成RFM，利用虚拟CCD建立的稳态重成像模型，采用与地形无关的方法，生成虚拟CCD的RFM，作为虚拟CCD的几何成像模型；
5.3原始影像与虚拟影像的坐标映射，步骤如下：
5.31)对于原始影像像点(s,l)，通过探元号s确定探元指向角，通过行号l确定成像时间，再通过成像时间内插确定成像姿态和位置，建立该像点成像模型，并通过前方交会到高程面H上，获取对应物点的坐标(B,L,H)，高程H可以从参考DEM中内插获取或者采用平均高程；
5.32)利用虚拟影像的RFM反算物点(B,L,H)在虚拟影像上的位置(s’,l’)；
5.33)由此可获得虚拟影像上(s’,l’)的灰度，即原始影像像点(s,l)的灰度值；
5.34)重复上述流程，获取每一个虚拟像点的灰度值；
5.35)对ROI全色和多光谱的各个波段均进行上述处理，获得空间分辨率一致的全色和多光谱影像。
作为本发明进一步的方案：步骤6所述ROI全色多光谱融合处理方法如下：
6.1高斯滤波参数自适应优化，寻找一个最优的高斯滤波参数σ，使得模糊处理后得到的低分辨率全色影像Pan′ds与多光谱影像MS的清晰度最为相似，步骤如下：
6.11)计算全色及多光谱影像各波段的平均梯度与均值，平均梯度定义为：
6.12)以全色影像的均值为标准，调整多光谱各波段的平均梯度，将多光谱各波段乘上一个均值调整系数μ，使其均值与全色相同，均值调整系数μ定义为：
其中
是全色影像的均值，
是多光谱第i波段的均值；
由方差与平均梯度的定义可知，多光谱各波段调整后的平均梯度为：
AGm＝μ*AG (4)；
其中AGm是经过均值调整后的平均梯度；
6.13)拟合高斯参数σ值与高斯滤波得到Pan′ds平均梯度间的函数关系，从1开始，直到0.5，每隔0.1设置一次σ值，以不同的σ值进行高斯滤波得到Pan′ds并计算其平均梯度，再对得到的采样数据进行二次多项式拟合，即可得到σ与Pan′ds平均梯度的函数关系；
6.14)将6.12)中计算的多光谱各波段调整后的平均梯度代入6.13)中得到的拟合方程中，计算出最优参数；
6.2色彩偏离的极值抑制，像素级融合中，融合规则不会考虑图像的空间性与整体性，使得图像的每个像素值都会产生变化，在地物边缘处，由于本身灰度变化剧烈，融合前后的像素值差异较大，目前的遥感影像多是以10bits或12bits存储的，而图像的显示是在8bits的显示器上显示，在转位的过程中会进行图像整体灰度的拉伸，当融合后某些像素值较真实值偏差较多时，会干扰图像整体的灰度直方图，从而引起图像整体的色差，通过原始多光谱影像计算图像各行各列的极值分布，从而得到一个极值矩阵；将极值矩阵引入到融合过程中，限制像素值变化范围，以达到抑制融合色差的目的；
6.3空间信息增强调整；融合影像与原始多光谱影像的各波段相关系数越接近，其融合前后的色差就越小，以相关系数作为调整依据，以某一波段作为基准，以空间信息增强作为手段调整其他波段的清晰度，相关系数定义如下：
其中ffus是融合图像，g是多光谱影像，
和
是图像相应的均值；空间信息增强可以由以下公式进行描述：
通过调整μ1和μ2来达到空间信息增强的目的，在空间信息增强的过程中，波段的相关系数会逐渐降低，当各波段的相关系数趋于一致时，图像的色彩保持效果及图像清晰度将达到最优化。
作为本发明进一步的方案：步骤7所述目标融合影像的几何纠正处理方法，对于融合后的目标影像，基于全球DEM数据，采用分块间接法几何纠正，生成带地理坐标的2级目标融合产品。
与现有技术相比，本发明的有益效果是：本发明基于目标区域提出一种光学卫星快速融合产品制作方法，直接从数据预处理阶段进行目标的快速生产，输出目标区域面向融合的5谱段同分辨率影像产品，融合处理直接基于目标区域的全色谱段和4个多光谱谱段进行。该处理模式可以避免大量不必要数据处理和数据存储，直接面向任务获取需要的信息，相比于地面常规处理，具有数据量小、计算量小、任务响应迅速的优势。
附图说明
图1为本发明实施例的基于目标区域的光学卫星快速融合产品制作方法流程图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的阐述。
一种基于目标区域的光学卫星快速融合产品制作方法，包括以下步骤：
步骤1，解析光学卫星下传原始数据中的辅助数据，包括成像时间数据、姿态数据、轨道数据，为ROI目标融合产品快速处理提供几何参数。
步骤2，确定ROI成像时间窗口。对于光学线阵推扫成像，只要确定了成像时间即可确定姿态数据、轨道数据范围，同时也可以得到成像行范围。采用一种动态搜索方式确定成像时间窗口，实施例中包含以下步骤：
2.1在数据解析过程中，从T0时刻开始，以一定时间间隔Δt0，对当前成像区域位置进行实时计算，即计算推扫成像行两端的经纬度坐标，判断ROI目标点经纬度是否落于计算的经纬度范围内，从而获取ROI的概略时间窗口。步骤如下：
1)建立当前成像行的成像时间T0的严格几何成像模型。根据当前成像行的成像时间T0确定轨道参数和姿态参数，再联合相机参数即可建立(相机参数通常会在发射前进行严格的实验室检校，但卫星发射和运行空间环境极其复杂，会导致星上载荷结构发生变化，需进行相机参数的在轨几何定标，才能保证ROI区域的精准定位，此处用到的相机参数即为几何定标后的精确参数)，如下式。
其中
式中(Xg,Yg,Zg)与(Xgps,Ygps,Zgps)分别表示像点对应的物方点及GPS(全球定位系统)天线相位中心在WGS84坐标系下的坐标，后者由卫星上搭载的GPS获取，按成像时间内插得到的轨道参数；
分别代表WGS84坐标系到J2000坐标系的旋转矩阵、J2000坐标系到卫星本体坐标系的旋转矩阵、卫星本体坐标系到相机坐标系的旋转矩阵，
由星敏、陀螺通过组合定姿得到，
由相机参数文件得到；(BX,BY,BZ)body代表从传感器投影中心到GPS天线相位中心的偏心矢量在卫星本体坐标系下的坐标，则在卫星发射前由实验室检校；λ为比例系数；(ψx(s),ψy(s))代表探元s在相机坐标系下的指向角，s代表CCD中的探元号，axi，ayi(i＝0,1,2,3)为多项式系数(相机参数文件中的内定标参数)。
2)计算当前行首尾端点的地理坐标，得到点p0，q0。
3)T1＝T0+Δt0时刻，重复上述1)和2)步骤，得到点p1，q1。
4)判断ROI目标点经纬度是否落于矩形(p0,q0,p1,q1)内，如不是，则Δt0后继续重复上述计算；如是，则停止计算，得到ROI的概略时间窗口(Ti,Ti+1)；
2.2在2.1中确定的概略时间窗口(Ti,Ti+1)内，以Ti时刻为起点，进一步缩小时间间隔为Δt1，重复2.1中步骤，并计算推扫成像行两端的经纬度坐标，判断目标点经纬度是否落于两次计算经纬度范围内，从而获取更加精确的ROI时间窗口；
2.3如此迭代，采用更小的时间间隔Δti，直至获取的时间窗口能完全准确的覆盖ROI的成像范围。为了提高搜索的可靠性，最小时间间隔要能覆盖ROI的成像范围。
步骤3，确定ROI成像影像范围。除了获取成像行范围外，还需获取成像列范围，对于具有多片TDI CCD(时间延时积分电荷耦合器件)的相机，还需要确定ROI对应于哪一片CCD影像以及准确列范围。实施例中包含以下步骤：
3.1依据步骤2得到的ROI成像时间窗口即可确定轨道参数和姿态参数，再联合相机参数即可按式(1)建立在该时间窗口内成像的每一片CCD影像的严格几何成像模型。
3.2利用严格几何成像模型计算得到虚拟控制点，通过地形独立法将严格几何成像模型转化为有理函数模型RFM(有理函数模型)，以解决严格几何成像模型坐标反算难以实现的问题。
3.3利用RFM将ROI四角点的经纬度坐标反投至每一片CCD影像空间，得到ROI在图像中的坐标，判断ROI四角点坐标是否在影像成像范围内，从而确定ROI的成像CCD片号和具体的成像起始行列号。
步骤4，ROI影像解析处理。对ROI对应的原始数据范围进行影像解析处理，得到ROI原始的全色和多光谱图像。
步骤5，ROI全谱段一体化虚拟重成像。为了保证目标融合产品处理的精度与效率，基于空间几何约束关系，对全色与多光谱影像建立相同的像方与物方的对应关系，进行全谱段一体化虚拟重成像，生成相同分辨率的全谱段虚拟扫描景。实施例包含以下步骤：
5.1无畸变虚拟CCD的设计。高分辨率光学遥感卫星为了获取较大的幅宽，全色和多光谱相机均会在焦平面上安置多个CCD。针对这一问题，在重成像过程中，采用一种虚拟的CCD代替原有的多片CCD进行成像。虚拟CCD是一个覆盖整个成像视场的无畸变的理想CCD线阵，且虚拟CCD的主点、主距及探元尺寸均与原始CCD一致，其像素个数是所有分片CCD像素个数总和去除重叠像素个数。同时，为了避免高程误差带来的重成像误差，虚拟CCD的位置通常应尽可能与真实CCD位置接近，保证原始成像视角与虚拟成像视角的差异最小。
5.2虚拟CCD成像模型构建。步骤如下：
1)建立虚拟CCD的稳态重成像模型，与原始分片CCD的严格几何成像模型基本相同，不同之处在于该模型是在稳定成像状态下建立的，即虚拟扫描行的姿态角则是在滤波处理后的姿态角中内插获取。
2)生成RFM。利用虚拟CCD建立的稳态重成像模型，采用与地形无关的方法，生成虚拟CCD的RFM，作为虚拟CCD的几何成像模型。
5.3原始影像与虚拟影像的坐标映射。步骤如下：
1)对于原始影像像点(s,l)，通过探元号s确定探元指向角，通过行号l确定成像时间，再通过成像时间内插确定成像姿态和位置，建立该像点成像模型，并通过前方交会到高程面H上，获取对应物点的坐标(B,L,H)。高程H可以从参考DEM中内插获取或者采用平均高程。
2)利用虚拟影像的RFM反算物点(B,L,H)在虚拟影像上的位置(s’,l’)。
3)由此可获得虚拟影像上(s’,l’)的灰度，即原始影像像点(s,l)的灰度值。
4)重复上述流程，获取每一个虚拟像点的灰度值。
5)对ROI全色和多光谱的各个波段均进行上述处理，获得空间分辨率一致的全色和多光谱影像。
步骤6，ROI全色多光谱融合处理。为保证全色多光谱融合精度，采用一种高色彩保真的全色多光谱融合技术，实施例包含以下步骤：
6.1高斯滤波参数自适应优化。寻找一个最优的高斯滤波参数σ，使得模糊处理后得到的低分辨率全色影像Pan′ds与多光谱影像MS的清晰度最为相似，步骤如下：
1)计算全色及多光谱影像各波段的平均梯度与均值，平均梯度定义为：
2)以全色影像的均值为标准，调整多光谱各波段的平均梯度。将多光谱各波段乘上一个均值调整系数μ，使其均值与全色相同，均值调整系数μ定义为：
其中
是全色影像的均值，
是多光谱第i波段的均值。
由方差与平均梯度的定义可知，多光谱各波段调整后的平均梯度为：
AGm＝μ*AG (4)；
其中AGm是经过均值调整后的平均梯度。
3)拟合高斯参数σ值与高斯滤波得到Pan′ds平均梯度间的函数关系。从1开始，直到0.5，每隔0.1设置一次σ值，以不同的σ值进行高斯滤波得到Pan′ds并计算其平均梯度。再对得到的采样数据进行二次多项式拟合，即可得到σ与Pan′ds平均梯度的函数关系。
4)将2)中计算的多光谱各波段调整后的平均梯度代入3)中得到的拟合方程中，计算出最优参数。
6.2色彩偏离的极值抑制。像素级融合中，融合规则不会考虑图像的空间性与整体性，使得图像的每个像素值都会产生变化。在地物边缘处，由于本身灰度变化剧烈，融合前后的像素值差异较大。目前的遥感影像多是以10bits或12bits存储的，而图像的显示是在8bits的显示器上显示，在转位的过程中会进行图像整体灰度的拉伸。当融合后某些像素值较真实值偏差较多时，会干扰图像整体的灰度直方图，从而引起图像整体的色差。通过原始多光谱影像计算图像各行各列的极值分布，从而得到一个极值矩阵；将极值矩阵引入到融合过程中，限制像素值变化范围，以达到抑制融合色差的目的。
6.3空间信息增强调整。融合影像与原始多光谱影像的各波段相关系数越接近，其融合前后的色差就越小。以相关系数作为调整依据，以某一波段作为基准，以空间信息增强作为手段调整其他波段的清晰度。相关系数定义如下：
其中ffus是融合图像，g是多光谱影像，
和
是图像相应的均值；
空间信息增强可以由以下公式进行描述：
通过调整μ1和μ2来达到空间信息增强的目的，在空间信息增强的过程中，波段的相关系数会逐渐降低，当各波段的相关系数趋于一致时，图像的色彩保持效果及图像清晰度将达到最优化。
步骤7，目标融合影像的几何纠正处理。实施例中，对于融合后的目标影像，基于全球DEM数据，采用分块间接法几何纠正，生成带地理坐标的2级目标融合产品。
以上所述为本发明较佳实施例，对于本领域的普通技术人员而言，根据本发明的教导，在不脱离本发明的原理与精神的情况下，对实施方式所进行的改变、修改、替换和变型仍落入本发明的保护范围之内。
Claims (6)
1.一种基于目标区域的光学卫星快速融合产品制作方法，其特征在于，包括以下步骤：
步骤1，解析卫星下传原始数据中的辅助数据；
步骤2，确定ROI成像时间窗口；
步骤3，确定ROI成像影像范围；
步骤4，ROI影像解析处理；
步骤5，ROI全谱段一体化虚拟重成像；
步骤6，ROI全色多光谱融合处理；
所述ROI全色多光谱融合处理方法如下：
6.1高斯滤波参数自适应优化，寻找一个最优的高斯滤波参数σ，使得模糊处理后得到的低分辨率全色影像Pan′ds与多光谱影像MS的清晰度最为相似，步骤如下：
6.11)计算全色及多光谱影像各波段的平均梯度与均值，平均梯度定义为：
6.12)以全色影像的均值为标准，调整多光谱各波段的平均梯度，将多光谱各波段乘上一个均值调整系数μ，使其均值与全色相同，均值调整系数μ定义为：
其中
是全色影像的均值，
是多光谱第i波段的均值；
由方差与平均梯度的定义可知，多光谱各波段调整后的平均梯度为：
AGm＝μ*AG (4)；
其中AGm是经过均值调整后的平均梯度；
6.13)拟合高斯参数σ值与高斯滤波得到Pan′ds平均梯度间的函数关系，从1开始，直到0.5，每隔0.1设置一次σ值，以不同的σ值进行高斯滤波得到Pan′ds并计算其平均梯度，再对得到的采样数据进行二次多项式拟合，即可得到σ与Pan′ds平均梯度的函数关系；
6.14)将6.12)中计算的多光谱各波段调整后的平均梯度代入6.13)中得到的拟合方程中，计算出最优参数；
6.2色彩偏离的极值抑制，像素级融合中，融合规则不会考虑图像的空间性与整体性，使得图像的每个像素值都会产生变化，在地物边缘处，由于本身灰度变化剧烈，融合前后的像素值差异较大，目前的遥感影像是以10bits或12bits存储的，而图像的显示是在8bits的显示器上显示，在转位的过程中会进行图像整体灰度的拉伸，当融合后某些像素值较真实值偏差较多时，会干扰图像整体的灰度直方图，从而引起图像整体的色差，通过原始多光谱影像计算图像各行各列的极值分布，从而得到一个极值矩阵；将极值矩阵引入到融合过程中，限制像素值变化范围，以达到抑制融合色差的目的；
6.3空间信息增强调整；融合影像与原始多光谱影像的各波段相关系数越接近，其融合前后的色差就越小，以相关系数作为调整依据，以某一波段作为基准，以空间信息增强作为手段调整其他波段的清晰度，相关系数定义如下：
其中ffus是融合图像，g是多光谱影像，
和
是图像相应的均值；
空间信息增强可以由以下公式进行描述：
通过调整μ1和μ2来达到空间信息增强的目的，在空间信息增强的过程中，波段的相关系数会逐渐降低，当各波段的相关系数趋于一致时，图像的色彩保持效果及图像清晰度将达到最优化；
步骤7，目标融合影像的几何纠正处理。
2.如权利要求1所述一种基于目标区域的光学卫星快速融合产品制作方法，其特征在于：步骤1所述解析卫星下传原始数据中的辅助数据包括成像时间数据、姿态数据、轨道数据，为ROI目标融合产品快速处理提供几何参数。
3.如权利要求2所述一种基于目标区域的光学卫星快速融合产品制作方法，其特征在于：步骤2所述确定R0I成像时间窗口的方法如下：
2.1在数据解析过程中，从T0时刻开始，以一定时间间隔Δt0，对当前成像区域位置进行实时计算，即计算推扫成像行两端的经纬度坐标，判断ROI目标点经纬度是否落于计算的经纬度范围内，从而获取ROI的概略时间窗口；步骤如下：
2.11)建立当前成像行的成像时间T0的严格几何成像模型，根据当前成像行的成像时间T0确定轨道参数和姿态参数，再联合相机参数即可建立，如下式：
其中
式中(Xg,Yg,Zg)与(Xgps,Ygps,Zgps)分别表示像点对应的物方点及GPS天线相位中心在WGS84坐标系下的坐标，后者由卫星上搭载的GPS获取，按成像时间内插得到的轨道参数；
分别代表WGS84坐标系到J2000坐标系的旋转矩阵、J2000坐标系到卫星本体坐标系的旋转矩阵、卫星本体坐标系到相机坐标系的旋转矩阵，
由星敏、陀螺通过组合定姿得到，
由相机参数文件得到；(BX,BY,BZ)body代表从传感器投影中心到GPS天线相位中心的偏心矢量在卫星本体坐标系下的坐标，则在卫星发射前由实验室检校；λ为比例系数；(ψx(s),ψy(s))代表探元s在相机坐标系下的指向角，s代表CCD中的探元号，axi，ayi为多项式系数，其中i＝0,1,2,3；
2.12)计算当前行首尾端点的地理坐标，得到点p0，q0；
2.13)T1＝T0+Δt0时刻，重复上述2.11)和2.12)步骤，得到点p1，q1；
2.14)判断ROI目标点经纬度是否落于矩形(p0,q0,p1,q1)内，如不是，则Δt0时候后继续重复上述计算；如是，则停止计算，得到ROI的概略时间窗口(Ti,Ti+1)；
2.2在2.1中确定的概略时间窗口(Ti,Ti+1)内，以Ti时刻为起点，进一步缩小时间间隔为Δt1，重复2.1中步骤，并计算推扫成像行两端的经纬度坐标，判断目标点经纬度是否落于两次计算经纬度范围内，从而获取更加精确的ROI时间窗口；
2.3如此迭代，采用更小的时间间隔Δti，直至获取的时间窗口能完全准确的覆盖ROI的成像范围，为了提高搜索的可靠性，最小时间间隔要能覆盖ROI的成像范围。
4.如权利要求3所述一种基于目标区域的光学卫星快速融合产品制作方法，其特征在于：步骤3所述确定R0I成像影像范围方法如下：
3.1依据步骤2得到的ROI成像时间窗口即可确定轨道参数和姿态参数，再联合相机参数即可按式(1)建立在该时间窗口内成像的每一片CCD影像的严格几何成像模型；
3.2利用严格几何成像模型计算得到虚拟控制点，通过地形独立法将严格几何成像模型转化为有理函数模型RFM，以解决严格几何成像模型坐标反算难以实现的问题；
3.3利用RFM将ROI四角点的经纬度坐标反投至每一片CCD影像空间，得到ROI在图像中的坐标，判断ROI四角点坐标是否在影像成像范围内，从而确定ROI的成像CCD片号和具体的成像起始行列号。
5.如权利要求4所述一种基于目标区域的光学卫星快速融合产品制作方法，其特征在于：步骤5所述ROI全谱段一体化虚拟重成像方法如下：
5.1无畸变虚拟CCD的设计，高分辨率光学遥感卫星为了获取较大的幅宽，全色和多光谱相机均会在焦平面上安置多个CCD，因此在重成像过程中，采用一种虚拟的CCD代替原有的多片CCD进行成像，虚拟CCD是一个覆盖整个成像视场的无畸变的理想CCD线阵，且虚拟CCD的主点、主距及探元尺寸均与原始CCD一致，其像素个数是所有分片CCD像素个数总和去除重叠像素个数，同时，为了避免高程误差带来的重成像误差，虚拟CCD的位置与真实CCD位置接近，保证原始成像视角与虚拟成像视角的差异最小；
5.2虚拟CCD成像模型构建，步骤如下：
5.21)建立虚拟CCD的稳态重成像模型，与原始分片CCD的严格几何成像模型基本相同，不同之处在于该模型是在稳定成像状态下建立，即虚拟扫描行的姿态角则是在滤波处理后的姿态角中内插获取；
5.22)生成RFM，利用虚拟CCD建立的稳态重成像模型，采用与地形无关的方法，生成虚拟CCD的RFM，作为虚拟CCD的几何成像模型；
5.3原始影像与虚拟影像的坐标映射，步骤如下：
5.31)对于原始影像像点(s,l)，通过探元号s确定探元指向角，通过行号l确定成像时间，再通过成像时间内插确定成像姿态和位置，建立该像点成像模型，并通过前方交会到高程面H上，获取对应物点的坐标(B,L,H)，高程H可以从参考DEM中内插获取或者采用平均高程；
5.32)利用虚拟影像的RFM反算物点(B,L,H)在虚拟影像上的位置(s’,l’)；
5.33)由此可获得虚拟影像上(s’,l’)的灰度，即原始影像像点(s,l)的灰度值；
5.34)重复上述5.31)-5.33)流程，获取每一个虚拟像点的灰度值；
5.35)对ROI全色和多光谱的各个波段均进行上述5.31)-5.34)处理，获得空间分辨率一致的全色和多光谱影像。
6.如权利要求5所述一种基于目标区域的光学卫星快速融合产品制作方法，其特征在于：步骤7所述目标融合影像的几何纠正处理方法，对于融合后的目标影像，基于全球DEM数据，采用分块间接法几何纠正，生成带地理坐标的2级目标融合产品。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011340149.6A CN112598608B (zh) | 2020-11-25 | 2020-11-25 | 一种基于目标区域的光学卫星快速融合产品制作方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202011340149.6A CN112598608B (zh) | 2020-11-25 | 2020-11-25 | 一种基于目标区域的光学卫星快速融合产品制作方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112598608A CN112598608A (zh) | 2021-04-02 |
| CN112598608B true CN112598608B (zh) | 2022-07-19 |
Family
ID=75184568
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202011340149.6A Active CN112598608B (zh) | 2020-11-25 | 2020-11-25 | 一种基于目标区域的光学卫星快速融合产品制作方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112598608B (zh) |
Families Citing this family (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113420982A (zh) * | 2021-06-23 | 2021-09-21 | 中国人民解放军63921部队 | 航天遥感测绘融合应用方法及系统 |
| CN113393499B (zh) * | 2021-07-12 | 2022-02-01 | 自然资源部国土卫星遥感应用中心 | 一种高分七号卫星全色影像和多光谱影像自动配准方法 |
| CN114332085B (zh) * | 2022-03-11 | 2022-06-24 | 西安中科西光航天科技有限公司 | 一种光学卫星遥感影像检测方法 |
| CN114723637B (zh) * | 2022-04-27 | 2024-06-18 | 上海复瞰科技有限公司 | 一种色差调整方法及系统 |
| CN114972545B (zh) * | 2022-08-01 | 2022-11-04 | 潍坊绘圆地理信息有限公司 | 一种高光谱卫星的在轨数据快速预处理方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104050643A (zh) * | 2014-06-18 | 2014-09-17 | 华中师范大学 | 遥感影像几何与辐射一体化相对校正方法及系统 |
| CN108960483A (zh) * | 2017-05-24 | 2018-12-07 | 电视广播有限公司 | 卫星调度方法、处理系统以及软件程序产品 |
| CN111612693A (zh) * | 2020-05-19 | 2020-09-01 | 中国科学院微小卫星创新研究院 | 一种旋转大幅宽光学卫星传感器校正方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US9554063B2 (en) * | 2015-06-12 | 2017-01-24 | Google Inc. | Using infrared images of a monitored scene to identify windows |
-
2020
- 2020-11-25 CN CN202011340149.6A patent/CN112598608B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104050643A (zh) * | 2014-06-18 | 2014-09-17 | 华中师范大学 | 遥感影像几何与辐射一体化相对校正方法及系统 |
| CN108960483A (zh) * | 2017-05-24 | 2018-12-07 | 电视广播有限公司 | 卫星调度方法、处理系统以及软件程序产品 |
| CN111612693A (zh) * | 2020-05-19 | 2020-09-01 | 中国科学院微小卫星创新研究院 | 一种旋转大幅宽光学卫星传感器校正方法 |
Non-Patent Citations (3)
| Title |
|---|
| Geometric stitching of a HaiYang-1C ultra violet imager with a distorted virtual camera;Jinshan Cao.et.;《Optics express》;20200424;第28卷(第9期);第14109-14126页 * |
| ZY-3卫星全色与多光谱影像融合方法比较;李霖等;《农业工程学报》;20140831;第30卷(第16期);第157-165页 * |
| 光学线阵相机在轨几何定标平台设计与实现;常学立等;《航天返回与遥感》;20150430;第36卷(第2期);第24-31页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112598608A (zh) | 2021-04-02 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN112598608B (zh) | 一种基于目标区域的光学卫星快速融合产品制作方法 | |
| Grodecki et al. | IKONOS geometric accuracy | |
| Teo | Bias compensation in a rigorous sensor model and rational function model for high-resolution satellite images | |
| Teo et al. | DEM-aided block adjustment for satellite images with weak convergence geometry | |
| CN112419380B (zh) | 一种基于云掩膜的静止轨道卫星序列影像高精度配准方法 | |
| CN109696182A (zh) | 一种星载推扫式光学传感器内方位元素定标方法 | |
| CN111912430B (zh) | 高轨光学卫星的在轨几何定标方法、装置、设备及介质 | |
| US20030086604A1 (en) | Three-dimensional database generating system and method for generating three-dimensional database | |
| CN106887016B (zh) | 一种gf-4卫星序列图像自动相对配准方法 | |
| CN114972545B (zh) | 一种高光谱卫星的在轨数据快速预处理方法 | |
| CN111473802A (zh) | 一种基于线阵推扫的光学传感器内方位元素定标方法 | |
| CN113514829A (zh) | 面向InSAR的初始DSM的区域网平差方法 | |
| CN114241064B (zh) | 一种遥感卫星内外方位元素实时几何定标方法 | |
| CN114708211B (zh) | 一种卫星遥感成像的优化处理方法及系统 | |
| Re et al. | Performance evaluation of 3DPD, the photogrammetric pipeline for the CaSSIS stereo images | |
| CN110516588B (zh) | 一种遥感卫星系统 | |
| Hu et al. | Planetary3D: A photogrammetric tool for 3D topographic mapping of planetary bodies | |
| Crespi et al. | GeoEye-1: Analysis of radiometric and geometric capability | |
| Seo et al. | Kompsat-2 direct sensor modeling and geometric calibration/validation | |
| Reinartz et al. | Accuracy analysis for DEM and orthoimages derived from SPOT HRS stereo data without using GCP | |
| Reulke et al. | Improvement of spatial resolution with staggered arrays as used in the airborne optical sensor ADS40 | |
| Tonolo et al. | Georeferencing of EROS-A1 high resolution images with rigorous and rational function model | |
| Schläpfer et al. | Parametric geocoding of AVIRIS data using a ground control point derived flightpath | |
| Jacobsen | Comparison of image orientation by IKONOS, QuickBird and OrbView-3 | |
| Zhu et al. | Roi-Orientated Sensor Correction Based on Virtual Steady Reimaging Model for Wide Swath High Resolution Optical Satellite Imagery |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |