CN105427267B - 获取刃边法测量光学遥感载荷在轨mtf精度的方法 - Google Patents
获取刃边法测量光学遥感载荷在轨mtf精度的方法 Download PDFInfo
- Publication number
- CN105427267B CN105427267B CN201510711683.6A CN201510711683A CN105427267B CN 105427267 B CN105427267 B CN 105427267B CN 201510711683 A CN201510711683 A CN 201510711683A CN 105427267 B CN105427267 B CN 105427267B
- Authority
- CN
- China
- Prior art keywords
- mrow
- mtf
- sword
- sub
- reference picture
- 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
- 238000005259 measurement Methods 0.000 title claims abstract description 73
- 230000003287 optical effect Effects 0.000 title claims abstract description 69
- 238000000034 method Methods 0.000 title claims abstract description 67
- 238000012360 testing method Methods 0.000 claims abstract description 18
- 238000013507 mapping Methods 0.000 claims description 36
- 230000007850 degeneration Effects 0.000 claims description 29
- 239000011159 matrix material Substances 0.000 claims description 26
- 238000009792 diffusion process Methods 0.000 claims description 12
- 238000010606 normalization Methods 0.000 claims description 9
- 241000208340 Araliaceae Species 0.000 claims description 7
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 claims description 7
- 235000003140 Panax quinquefolius Nutrition 0.000 claims description 7
- 235000008434 ginseng Nutrition 0.000 claims description 7
- 238000012545 processing Methods 0.000 claims description 7
- 238000004422 calculation algorithm Methods 0.000 claims description 4
- 239000000284 extract Substances 0.000 claims description 3
- 238000006116 polymerization reaction Methods 0.000 claims description 2
- 101100258233 Caenorhabditis elegans sun-1 gene Proteins 0.000 claims 1
- 101100024583 Mus musculus Mtf1 gene Proteins 0.000 claims 1
- 240000002853 Nelumbo nucifera Species 0.000 claims 1
- 235000006508 Nelumbo nucifera Nutrition 0.000 claims 1
- 235000006510 Nelumbo pentapetala Nutrition 0.000 claims 1
- 235000013399 edible fruits Nutrition 0.000 claims 1
- 238000000605 extraction Methods 0.000 abstract description 6
- 238000004458 analytical method Methods 0.000 abstract description 3
- 238000003384 imaging method Methods 0.000 description 10
- 238000005070 sampling Methods 0.000 description 10
- 238000011156 evaluation Methods 0.000 description 7
- 230000008569 process Effects 0.000 description 6
- 238000010586 diagram Methods 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 4
- 230000015556 catabolic process Effects 0.000 description 3
- 238000006731 degradation reaction Methods 0.000 description 3
- 238000013461 design Methods 0.000 description 3
- 230000000694 effects Effects 0.000 description 3
- 238000005457 optimization Methods 0.000 description 3
- 230000009466 transformation Effects 0.000 description 3
- 238000013316 zoning Methods 0.000 description 3
- 238000012935 Averaging Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000006872 improvement Effects 0.000 description 2
- VMXUWOKSQNHOCA-UKTHLTGXSA-N ranitidine Chemical compound [O-][N+](=O)\C=C(/NC)NCCSCC1=CC=C(CN(C)C)O1 VMXUWOKSQNHOCA-UKTHLTGXSA-N 0.000 description 2
- 230000009897 systematic effect Effects 0.000 description 2
- 238000012546 transfer Methods 0.000 description 2
- 238000004364 calculation method Methods 0.000 description 1
- 230000001186 cumulative effect Effects 0.000 description 1
- 238000003708 edge detection Methods 0.000 description 1
- 239000004744 fabric Substances 0.000 description 1
- 238000000556 factor analysis Methods 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012634 optical imaging Methods 0.000 description 1
- 238000011002 quantification Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 238000005316 response function Methods 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 238000010206 sensitivity analysis Methods 0.000 description 1
- 238000000926 separation method Methods 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
- 230000003313 weakening effect Effects 0.000 description 1
Classifications
-
- 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
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
-
- 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/20—Special algorithmic details
- G06T2207/20172—Image enhancement details
-
- 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/30—Subject of image; Context of image processing
- G06T2207/30168—Image quality inspection
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
Abstract
本发明提供了一种获取基于刃边法测量光学遥感载荷在轨MTF的测量精度的方法。该方法基于实际测试的MTF值设置MTF参考值,针对刃边法MTF测试特点,基于实际测试时载荷获取的刃边图像提取关键参数生成参考图像,利用刃边法对参考图像进行计算得到MTF计算值,然后将计算值与MTF参考值对比,获得绝对误差和相对误差,由此得到MTF的测量精度(或不确定度),从而解决了现有精度分析方法无法获得或估计刃边法实际测量光学遥感载荷在轨MTF测量精度的问题。
Description
技术领域
本发明涉及对地观测遥感技术领域,特别涉及一种获取基于刃边法测量光学遥感载荷在轨MTF的测量精度的方法。
背景技术
调制传递函数MTF(Modulation Transfer Function)是评价光学遥感载荷成像性能的一个重要指标,它表示了各个空间频率信号经过光学成像系统后调制度损失的百分比,反映成像系统在对目标成像过程中信号的扩散与削弱程度,是目前国际上普遍使用的评定光学遥感载荷成像性能以及空间分辨率的指标之一。
根据选用靶标的不同,目前测量光学遥感载荷在轨MTF的方法主要有刃边法、三线靶标法、点源/点阵法、脉冲法、辐射状靶标法。在这些方法中,刃边法因其能够获得不同空间频率的MTF曲线、能够更加完整地刻画光学遥感载荷空间响应特性、在实际应用中靶标布设与选取条件相对宽松(既可以是人工布设的靶标,也可以选取合乎要求的刃边类型的地物目标),因而是评估中高分辨率光学遥感载荷在轨MTF最普遍使用的方法。诸如SPOT5/6/7、IKONOS、QUICKBIRD、GEOEYE、ZY02/03-HR、GF1/2等高分辨率卫星的光学遥感载荷均采用刃边法进行在轨MTF评测。
为了提高刃边法评测载荷在轨MTF的精度,各卫星运行机构相继布设高质量的刃边靶标,基于刃边法的MTF在轨评测技术也日趋成熟。图1为国内外典型的高分辨率刃边靶标。其中,A为美国德克萨斯州Big Spring靶标;B为台湾澎湖靶标;法国Salon de Province靶标;D为加拿大Mt.Albert,Ontario靶标。
图2为现有技术采用刃边法评测光学遥感载荷在轨MTF的数据处理流程图。具体处理流程如下:
步骤S101,提取刃边子图像:从经过相对辐射校正的含有刃边靶标或刃边地物(以下统称为刃边靶标)的图像数据中提取不小于15像元×15像元的有效计算区域(刃边子图像);特别说明一下,根据刃边子图像中刃边的方向不同,可以分别测量得到载荷阵列方向(对应于所获取图像的x方向)的MTF和载荷运动方向(对应于所获取图像的y方向)的MTF,如图3所示,为了描述简单,后续描述如不特别说明,均针对载荷阵列方向的MTF,载荷运动方向的MTF计算与精度测量按照同样方法参照实现。
步骤S102,边缘探测:提取刃边子图像的每一行的灰度数据,采用一定的模型进行差值,利用微分法、质心法、Hough变换法、Fermi函数法等方法,确定该行数据灰度变换廓线的边缘位置到亚像元精度(刃边边缘的中心位置)。对每一行数据采用同样方法,获得刃边子图像每一行数据的亚像元精度的刃边边缘中心位置,然后通过直线拟合对探测到的边缘中心位置进行进一步调整,使探测到的各行边缘中心位置点在一条直线上,计算该拟合直线的斜率,可以得到刃边的倾角θ;
步骤S103,提取边扩散函数ESF(Edge Spread Function):利用步骤S102探测到的每一行的刃边边缘中心位置,对刃边子图像的各行数据进行配准、合并、插值处理,得到亚像元精度的ESF;
步骤S104,求取线扩散函数LSF(Line Spread Function):对步骤S103得到的ESF进行差分,并采用窗函数进行截断处理,得到宽度为-5像元~+5像元(或-10像元~+10像元)亚像元精度的LSF;
步骤S105,MTF求取:对步骤S104得到的LSF进行傅里叶变换,对变换结果取模值并做归一化处理,即可得到MTF序列曲线和Nyquist频率下的MTF值MTFNyquist(对宽度为-5像元~+5像元的LSF,其MTFNyquist为归一化MTF序列中的第6个值;对宽度为-10像元~+10像元的LSF,其MTFNyquist为归一化MTF序列中的第11个值)。
作为一种普遍使用的光学遥感载荷在轨MTF测量方法,现有的文献更多地涉及到MTF评测方法的改进和优化;也有部分文献采用仿真方法分析了亮暗区域对比度、靶标倾角、计算区域尺寸、噪声等因素对于MTF算法评测精度的影响,但是上述分析并没有与载荷实际成像条件结合起来。比如,文献中通常假设刃边靶标的亮暗区域的灰度值分别为200、100,假设噪声的标准差为0.7、1、1.5,靶标倾角为5°;此外,文献中通常假设刃边亮暗区域的噪声水平是相同的,并且往往没有考虑随机相位的影响。这些假设往往与载荷实际成像条件及获取图像并不一致,因而无法真实反映光学遥感载荷在轨实际成像条件下刃边法评测MTF的精度。
发明内容
(一)要解决的技术问题
本发明提出了一种基于载荷实际成像条件获取刃边法测量光学遥感载荷在轨MTF测量精度的方法,以解决现有精度分析方法无法获得或估计刃边法实际测量光学遥感载荷在轨MTF测量精度的问题。
(二)技术方案
本发明获取基于刃边法测量光学遥感载荷在轨MTF在第一方向的测量精度的方法包括:
步骤A:在光学遥感载荷获取的实测刃边靶标图像中截取尺寸为m行n列的多个刃边子图像,由该多个刃边子图像获取光学遥感载荷在第一方向的在轨MTF参考值和刃边倾角参考值;
步骤B:由光学遥感载荷在第一方向的在轨MTF参考值,获取光学遥感载荷在该第一方向的扩散尺度,构建亚像元尺度下的参考图像退化模板psf(x,y);
步骤C:由光学遥感载荷获取的实测刃边靶标图像提取被刃边分割的两个均匀区域的灰度均值和随机噪声标准差;
步骤D:根据两个均匀区域的随机噪声标准差,获得用于获取后续步骤中亚像元水平下刃边噪声参考图像的噪声参数:刃边噪声参考图像中两个区域的噪声标准差;
步骤E:利用刃边靶标参数,构建亚像元水平的第一方向的理想刃边参考图像I(x,y),其中,所述刃边靶标参数包括:刃边子图像尺寸、靶标图像被刃边分割的两个均匀区域的灰度均值、光学遥感载荷在第一方向的刃边倾角参考值;
步骤F:将所述理想刃边参考图像I(x,y)与参考图像退化模板psf(x,y)进行二维卷积,截取卷积结果矩阵中有效数据,得到亚像元水平的第一方向的无噪声刃边退化图像Ipsf(x,y);
步骤G:对所述刃边退化图像Ipsf(x,y)中的数据进行数据聚合,生成整像元水平的第一方向的无噪声刃边退化参考图像g(x,y);
步骤H:利用所述刃边噪声参考图像中两个区域的噪声标准差,构建K幅亚像元水平的第一方向的带噪声刃边参考图像Ik(x,y),k=1,2,…K;
步骤I:将带噪声刃边参考图像Ik(x,y)k=1,2,…K,分别与参考图像退化模板psf(x,y)进行二维卷积,截取卷积结果矩阵中的有效数据得到K幅亚像元水平的第一方向的带噪声刃边退化图像Ipsf-k(x,y),k=1,2,…K;
步骤J:对所述带噪声刃边退化图像Ipsf-k(x,y),k=1,2,…K,进行数据聚合处理,生成K幅整像元水平的阵列方向带噪声退化参考图像gk(x,y),k=1,2,…K;以及
步骤K:利用刃边法MTF测试算法,对所述整像元水平阵列无噪声刃边退化图像g(x,y),以及所述K幅整像元水平的阵列方向带噪声刃边退化参考图像gk(x,y),k=1,2,…K,分别进行第一方向的MTF计算,将Nyquist频率处的归一化MTF计算值与第一方向的在轨MTF参考值进行对比,得到光学遥感载荷在轨MTF精度;
其中,所述第一方向为载荷阵列方向和载荷运动方向其中之一。
(三)有益效果
从上述技术方案可以看出,本发明法获取刃边法测量光学遥感载荷在轨MTF测量精度的方法具有以下有益效果:
(1)基于载荷实际成像条件下所获取的图像,提取刃边法在轨MTF评测中的关键参数(刃边亮暗区域灰度、刃边靶标倾角、MTF计算有效区域尺寸、噪声等),并考虑刃边两侧区域随机噪声水平不同以及随机采样相位的影响,得到刃边法测量光学遥感载荷在轨MTF的测量精度,该测量精度与MTF在轨实际测量密切相关,可为MTF在轨实际测量结果提供更可靠的定量化精度和置信度量;
(2)基于实测条件获得各种MTF评测方法的测量精度,便于更有效对比不同MTF评测方法的优劣。
附图说明
图1为国内外典型的高分辨率刃边靶标;
图2为现有技术采用刃边法评测光学遥感载荷在轨MTF的数据处理流程图;
图3为现有技术从刃边靶标图像中提取的刃边子图像(MTF计算有效区域)的示意图;
图4为根据本发明实施例获取刃边法评测光学遥感载荷在轨MTF精度的方法流程图;
图5A为测量载荷阵列方向MTF的刃边子图像的刃边及倾角示意图。
图5B为测量载荷飞行方向MTF的刃边子图像的刃边及倾角示意图;
图6为按照步骤B生成的退化模板示意图(σ1=σ2=0.61244);
图7为实施步骤E时参考图像P、原点O、倾角为θx-ref的刃边直线的示意图;
图8为从某星载光学载荷获取的阵列方向刃边靶标图像(方框内区域为用于阵列方向MTF计算的刃边子图像);
图9为基于图8进行MTF计算及所获得图像参数,生成用于获取其MTF精度所需要的刃边参考图像的过程示意图。
具体实施方式
本发明基于实际测试的MTF值设置MTF参考值,针对刃边法MTF测试特点,基于实际测试时载荷获取的刃边图像提取关键参数生成参考图像,利用刃边法对参考图像进行计算得到MTF计算值,然后将计算值与MTF参考值对比,获得绝对误差和相对误差,由此得到MTF的测量精度。
在本发明的一个示例性实施例中,提供了一种获取刃边法测量光学遥感载荷在轨MTF在阵列方向的测量精度的方法。需要说明的是,本实施例是以载荷阵列方向的MTF,即MTFx为例进行说明。关于载荷运动方向的MTF,即MTFy可按照同样方法参照实现。
图4为根据本发明实施例获取刃边法评测光学遥感载荷在轨MTF精度的方法流程图。如图4所示,本实施例获取光学遥感载荷在轨MTF精度的方法包括:
步骤A:由光学遥感载荷获取实测刃边靶标图像,在该实测刃边靶标图像中截取的尺寸为m×n的多个刃边子图像在归一化频率下的MTF值MTFx-Nyquist和刃边倾角θx经平均获取光学遥感载荷在阵列方向的在轨MTF参考值MTFx-ref和刃边倾角θx-ref;
MTFx-ref采用多次实际测试过程中测量得到的MTF值的平均值:即对测试中载荷获取的刃边靶标图像,采用刃边法计算载荷的归一化Nyqiust频率(f=0.5)下的MTF值(MTFx-Nyquist)作为MTFx-ref。由于受噪声和计算区域等不确定因素影响,每次MTFx-Nyquist计算的结果会有差别,MTFx-ref取多次计算的平均值,步骤A又可以分为以下子步骤:
子步骤A1:在光学遥感载荷获取的实测刃边靶标图像中截取满足MTF计算要求的多个刃边子图像D;
其中,该多个刃边子图像D具有相同尺寸m行×n列(或接近相同尺寸),后续将基于该尺寸m行×n列构建参考图像。多个刃边子图像D在刃边靶标图像的位置并没有特别的限定。优选地,在选取刃边子图像D时保证刃边位于刃边子图像D中间位置效果会更好一些。
子步骤A2:对于多个刃边子图像中的每一个,采用刃边法计算其在载荷阵列方向的归一化Nyquist频率下的MTF,即MTFx-Nyquist,并在计算过程中获得刃边倾角θx;
其中,θx为刃边子图像中的刃边与图像列方向的夹角,范围为-90°~90°。图5A显示了测量载荷阵列方向MTF的4种刃边子图像中的刃边以及相应的刃边倾角θx。其中,图(a)和图(c)中由亮暗区域所形成的刃边向图像列方向偏左的方向倾斜,θx的取值范围为-90°~0°;图(b)和图(d)中由亮暗区域所形成的刃边向图像列方向偏右的方向倾斜,θx的取值范围为0°~90°。
对于测量载荷飞行方向MTF的刃边子图像,刃边倾角θy为子图像中的刃边与图像行方向的夹角,范围为-90°~90°。图5B显示了测量载荷飞行方向MTF的4种刃边子图像中的刃边以及相应的刃边倾角θy。其中,图(a)和图(c)中由亮暗区域所形成的刃边向图像行方向偏上的方向倾斜,θy的取值范围为0°~90°;图(b)和图(d)中由亮暗区域所形成的刃边向图像列方向偏下的方向倾斜,θy的取值范围为-90°~0°。
子步骤A3:对多个刃边子图像在载荷阵列方向的归一化Nyquist频率下的MTF,即MTFx-Nyquist和刃边倾角θx分别求平均,获取光学遥感载荷在载荷阵列方向的在轨MTF参考值MTFx-ref和刃边倾角参考值θx-ref。
步骤B:由光学遥感载荷在阵列方向的在轨MTF参考值MTFx-ref,获取光学遥感载荷在该方向的扩散尺度σx,构建亚像元尺度(水平)下的参考图像退化模板psf(x,y);
光学遥感载荷系统响应函数PSF通常采用高斯模型表示,结合光电系统数字采样效应,表示为:
其中:
其中,x,y分别表示光学遥感载荷在阵列方向和运动方向的坐标,σx、σy分别表示载荷系统函数PSF在阵列方向和运动方向的扩散尺度。则MTF参考值MTFref与该光学遥感载荷系统函数PSF的扩散尺度σ之间存在如下关系:
利用公式(2),已知MTF参考值,即可得到光学遥感载荷系统函数PSF的关键参数——扩散尺度。比如,当MTFx-ref=0.1时,经由式2-1计算,对应的系统函数PSF的扩散尺度σx为0.61244。
由公式(1-1)可知,光学遥感载荷的系统响应函数PSF在阵列方向和运动方向是可以分离的,当光学遥感载荷在运动方向的在轨MTF参考值MTFy-ref可以测量得到时,可以按照以上同样方法获得载荷运动方向的扩散尺度σy;如果无法测量MTFy-ref时,可以令MIFy-ref=MTFx-ref,从而有σy=σx,将不会影响本发明在阵列方向MTF的测量精度。
基于系统函数PSF的扩散尺度σx、σy,可以生成亚像元水平的参考图像归一化退化模板为:
本实施例中,参考图像退化模板的尺寸选为整像元水平下的10像元×10像元(-5像元≤x≤4.9像元,-5像元≤y≤4.9像元),两数据点之间的间隔为0.1像元。图6为按照步骤B生成的亚像元水平参考图像退化模板psf(x,y)的示意图,图中,psf(x,y)的尺寸选为10像元×10像元(-5像元≤x≤4.9像元,-5像元≤y≤4.9像元),扩散尺度σx=σy=0.61244,该模板对应的MTFx-rey=MTFy-ref=0.1。
需要说明的是,本发明并不以上述实施例为限,本领域技术人员可以合理设置参考图像退化模板的尺寸以及数据点之间的间隔。但需要注意的是,参考图像退化模板psf(x,y)中数据点之间的间隔应当介于0.01~0.1像元之间。
步骤C:由光学遥感载荷获取的实测刃边靶标图像提取被刃边分割的两个均匀区域的灰度均值和随机噪声标准差,其中,设该刃边靶标图像中刃边左侧均匀区域ZL的平均灰度为μL,随机噪声标准差σnoise-L;刃边右侧均匀区域ZR的平均灰度为μR,ZR的随机噪声标准差σnoise-R;
本步骤中,上述实测刃边靶标图像中被刃边分割的两个均匀区域分为以下两种情况:当测量载荷阵列方向的MTF时,实测刃边靶标图像被刃边分割成左右两个均匀区域,刃边左侧为亮区域,右侧为暗区域,或者刃边左侧为暗区域,右侧为亮区域;当测量飞行方向的MTF时,图像被刃边分割成上下两个均匀区域,刃边上部为亮区域,下部为暗区域,或者刃边上部为暗区域,下部为亮区域。
本实施例中,按照阵列方向MTF测量进行描述。具体的方法是:基于光学遥感载荷获取的实测刃边靶标图像,在靶标图像中的刃边两侧的均匀区域中分别选取一个尽可能大的矩形区域(ZL和ZR),为避免邻近区域的干扰,选取的区域应当满足:与图像中的刃边中心位置距离大于3个像元,与靶标的各个边缘的距离大于2个像元。分别计算这两个所选取矩形区域的灰度均值和标准差,得到该图像中刃边左右两侧均匀区域的灰度均值μL、μR和区域ZL、ZR的随机噪声标准差σnoise-L、σnoise-R。该过程对于本领域技术人员是公知的,在相关的技术文件中也可以得到,此处不再进行赘述。
步骤D:根据刃边左侧均匀区域ZL(mZL行×nZL列)的随机噪声标准差σnoise-L和刃边右侧均匀区域ZR(mZR行×nZR列)的随机噪声标准差σnoise-R,获得用于构建后续步骤中亚像元水平下刃边噪声参考图像的噪声参数:刃边噪声参考图像中刃边左右两个区域的噪声标准差σnL和σnR;
其中,获取刃边噪声参考图像刃边左区域的噪声标准差σnL可以分为以下子步骤:
子步骤D1:构建一个均值为0、标准差σINZL的初始值为σnoise-L、尺寸为[(mZL+10)×10行]×[(nZL+10)×10列]的亚像元水平的噪声矩阵INZL(x,y),该噪声矩阵INZL(x,y)对应的数据点间距为0.1像元;
子步骤D2:令噪声矩阵INZL(x,y)与步骤B得到的参考图像退化模板Psf(x,y)进行卷积:INZL(x,y)*psf(x,y),截取卷积结果矩阵中最中间的有效数据((mZL×10+1)行×(nZL×10+1)列)的数据(在Matlab中,使用卷积函数conv2的‘valid’参数)IPNZL;
子步骤D3:对矩阵IPNZL中的数据每10×10个数据为一组进行聚合,得到一个mZL行×nZL列的矩阵PNZL,PNZL数据点间距为1像元;
子步骤D4:计算PNZL的噪声标准差σPNZL;
子步骤D5:对比σPNZL和σnoise-L,如果σPNZL不等于σnoise-L,则调整子步骤D1中噪声矩阵INZL(x,y)的标准差σINZL,然后重复子步骤D2~D5,直到σPNZL=σnoise-L。
则此时步骤D1中噪声矩阵INZL(x,y)的标准差σINZL,即为后续构建亚像元水平下刃边参考图像刃边左区域噪声标准差σnL。
采用同样的方法,获取刃边参考图像刃边右区域的噪声标准差σnR。
步骤E:利用步骤A、B、C、D中提取的刃边靶标参数,包括:刃边子图像尺寸,即行数m和列数n;靶标图像中被刃边分割的左右两侧均匀区域的灰度均值μL、μR;光学遥感载荷在载荷阵列方向的刃边倾角参考值θx-ref,构建亚像元水平的阵列方向的理想刃边参考图像I(x,y);
步骤E又可以分为以下子步骤:
子步骤E1:基于步骤A提取的刃边子图像D的尺寸(m行n列,数据点间隔为1像元,构建一个(m+1+10)×10行,(n+1+10)×10列的亚像元水平的参考图像P,参考图像P的数据点(像素点)间距为0.1像元。
子步骤E2,在参考图像P中,令与图像P的列方向(如图7所示的y轴方向)成θx-ref角度的一条直线l,从参考图像P中第 像素的左上角点(以下简称O点)斜穿而过,则根据该直线l穿越P中某一像素时,分割该像素为左右两个部分,设参考图像P中各个像素的面积为1,根据P中各个像素被直线l分割的左右两部分的面积为权重,对参考图像P进行灰度赋值:
(1)对于P中没有被直线l穿过的像素p(x,y),如果p(x,y)位于直线l左侧,则p的灰度I(x,y)赋值为:
I(x,y)=μL
(2)对于参考图像P中没有被直线l穿过的像素p(x,y),如果p(x,y)位于直线l右侧,则该p(x,y)的灰度I(x,y)赋值为:
I(x,y)=μR
(3)对于被直线l穿过的像素p(x,y),则计算直线l穿过p时将p分割成为左右两部分,所占面积分别为S1、S2(S1+S2=1),则该像素p(x,y)的灰度I(x,y)赋值为:
I(x,y)=S1·μL+S2·μR
图7显示了实施步骤E时构建的参考图像P、原点O、倾角为θx-ref的刃边直线的示意图。图中,O点四周的四个像素分别为O1、O2、O3、O4,O1没有被直线l穿过,且位于l左侧,所以O1的灰度值为μL;O2被直线l穿过,设O2面积为1,直线l分割O2左右两部分的面积分别为:tgθx-ref和1-tgθx-ref,则O2的灰度值为μL·tgθx-ref+μR·(1-tgθx-ref);O3同样被直线l穿过,直线l分割O3左右两部分的面积分别为:1-tgθx-ref和tgθx-ref,则O3的灰度值为μL·(1-tgθx-ref)+μR·tgθx-ref;O4没有被直线l穿过,且位于l右侧,所以O4的灰度值为μR。
由此,构建出一幅亚像元水平的载荷阵列方向理想刃边参考图像I(x,y)。
步骤F:将步骤E得到的亚像元水平的阵列方向的理想刃边参考图像I(x,y)与步骤B得到的参考图像退化模板psf(x,y)进行二维卷积,截取卷积结果矩阵中有效数据(在Matlab中,使用卷积函数conv2的‘valid’参数),得到(m+1)×10+1行、(n+1)×10+1列的亚像元水平的阵列方向的无噪声刃边退化图像Ipsf(x,y);
步骤G:对步骤F得到的亚像元水平的阵列方向的刃边退化图像Ipsf(x,y),以Ipsf(x,y)中位于行、列的像素为起始位置,在图像Ipsf(x,y)中上下左右每10×10个像素为1组进行数据聚合,生成整像元水平的阵列方向的无噪声刃边退化参考图像g(x,y);
步骤H:利用步骤D得到的亚像元水平下刃边参考图像刃边左右两侧区域的随机噪声标准差σnL、σnL,构建K幅(K通常取100)亚像元水平的阵列方向的带噪声刃边参考图像Ik(x,y),k=1,2,…K;
步骤H又可分为以下子步骤:
子步骤H1:按照子步骤E1所述方法,再次构建K幅(m+1+10)×10行、In+1+10)×10列的亚像元水平的参考图像Pk,k=1,2,…K,参考图像Pk的数据点(像素点)间距为0.1像元。
子步骤H2:分别构建K个均值为0、标准差σnL、尺寸为(m+1+10)×10行、(n+1+10)×10列的噪声矩阵NPLk,k=1,2,…K,和K个均值为0、标准差σnR、尺寸为(m+1+10)×10行、(n+1+10)×10列的噪声矩阵NPRk,k=1,2,…K,噪声矩阵NPLk和NPRk的数据点间距为0.1像元;
子步骤H3,按照类似子步骤E2的方法,在参考图像Pk,k=1,2,…K中,令与参考图像Pk的列方向(如图7所示的y轴方向)成θx-ref角度的一条直线l从参考图像Pk中第像素的左上角点(以下简称O点)斜穿而过,则根据该直线l穿越Pk中某一像素时,分割该像素为左右两个部分,设参考图像Pk中各个像素的面积为1,根据参考图像Pk中各个像素被直线l分割的左右两部分的面积为权重,对参考图像Pk进行灰度赋值:
(1)对于Pk中没有被直线l穿过的像素pk(x,y),如果pk(x,y)位于直线l左侧,则pk(x,y)的灰度Ik(x,y)赋值为:
Ik(x,y)=μL+NPLk(x,y);
(2)如果pk(x,y)位于直线l右侧,则该Pk(x,y)的灰度Ik(x,y)赋值为:
Ik(x,y)=μR+NPRk(x,y)。
(3)对于被直线l穿过的像素pk(x,y),则计算直线l穿过pk时将pk分割成为左右两部分,所占面积分别为S1k、S2k(S1k+S2k=1),则该像素pk(x,y)的灰度Ik(x,y)赋值为:
Ik(x,y)=S1k·[μL+NPLk(x,y)]+S2k·[μR+NPRk(x,y)]
由此,构建出K幅亚像元水平的阵列方向带噪声刃边参考图像Ik(x,y),k=1,2,…K。
步骤I:将步骤H得到的亚像元水平的阵列方向带噪声刃边参考图像Ik(x,y),k=1,2,…K,分别与步骤B得到的参考图像退化模板psf(x,y)进行二维卷积,截取卷积结果矩阵中最中间的有效数据(在Matlab中,使用卷积函数conv2的‘valid’参数),得到K幅(m+1)×10+1行、(n+1)×10+1列的亚像元水平的阵列方向带噪声刃边退化图像Ipsf-k(x,y),k=1,2,…K;
步骤J:对步骤I得到的亚像元水平的阵列方向的带噪声刃边退化图像Ipsf-k(x,y),k=1,2,…K,进行数据聚合处理,生成K幅整像元水平的阵列方向带噪声退化参考图像gk(x,y),k=1,2,…K;
需要注意的是,采样相位是MTF测量中另一个必须考虑的随机因素。载荷在对地物成像并以像元为单位进行数字化采样时,在不同的采样相位下,获取的图像灰度值会不同,并最终影响MTF的计算结果,因而在步骤J生成整像元水平的带噪声退化参考图像gk(x,y)k=1,2,…K时,需要将随机采样相位因素引入。步骤K又可以分为以下子步骤:
子步骤J1:在(-5、-4、-3、-2、-1、0、1、2、3、4、5)这11个数中每次随机选出2个数,共随机选出K对数值(Sxk,Syk)k=1,2,…K,作为随机相位的偏移值;
子步骤J2:对第k幅亚像元水平的阵列方向带噪声刃边退化图像Ipsf-k(x,y),k=1,2,…K,以Ipsf-k(x,y)中第行、列的像素为聚合起始位置,在图像Ipsf-k(x,y)中上下左右每10×10个像素为1组进行聚合,即对该10×10个像素组内所有像素数据点的灰度数值累加求平均并四舍五入取整,最终得到K幅带有随机噪声和随机相位的、间距为1像元的整像元水平的阵列方向带噪声退化刃边参考图像gk(x,y),k=1,2,…K。
步骤G得到的整像元水平的阵列方向无噪声刃边退化图像,没有包含随机噪声和随机采样相位造成的随机误差,可以用于分析光学遥感载荷在轨MTF阵列方向的MTF系统误差。
步骤J得到的整像元水平的阵列方向带噪声刃边退化参考图像gk(x,y),k=1,2,…K,综合了MTF实测过程中刃边靶标亮暗反射区域不同辐亮度反射、刃边靶标倾角、随机噪声、CCD相机随机采样相位、CCD相机成像时的模糊退化效应、MTF计算图像区域尺寸等多重关键影响因素,因而使用MTF刃边算法针对gk(x,y),k=1,2,…K,计算得到的MTF值,最能够反映实际测量的情况。
图8和图9显示了按照上述步骤A~步骤J生成用于获取实测某星载光学载荷阵列方向MTF的测量精度时,所生成的一幅带噪声刃边参考图像的过程。图8显示了某星载光学载荷获取的刃边靶标图像,图中方框内的区域为用于测量阵列方向MTF的刃边子图像。通过该图像,按照步骤A得到MTFx-ref=0.070,θ=8.75°,子图像尺寸m=40,n=40;按照步骤C,得到μL=804.25,μR=191.68,σnoise-A=8.79,σnoise_B=6.16;按照步骤D,得到σnL=130,σnR=190。图9为基于图8进行MTF测量及所获得图像参数,生成用于获取其MTF精度所需要的刃边参考图像的过程示意图。其中,图9a为按照步骤E基于上述参数(μL=804.25,μR=191.68,θx-ref=8.75°,m=40,n=40)构建的亚像元水平阵列方向理想刃边参考图像I(x,y),图像尺寸为510×510;图9b为按照步骤H得到的1幅亚像元水平下有噪声参考图像I1(x,y),图像尺寸为510×510;图9c为按照步骤I得到的1幅亚像元水平下有噪声退化参考图像,图像尺寸为411×411;图9d为按照步骤J,经过数据聚合并考虑随机相位得到的整像元水平下有噪声退化参考图像,图像尺寸为40×40。
载荷飞行方向的无噪声和带噪声刃边退化参考图像可以参照同样的方法得到。
步骤K:利用刃边法MTF测试算法,对步骤G得到的整像元水平阵列无噪声刃边退化图像g(x,y),以及步骤J生成的K幅整像元水平的阵列方向带噪声刃边退化参考图像gk(x,y),k=1,2,…K,分别进行阵列方向的MTF计算,然后将Nyquist频率处的MTF计算值MTFcal与步骤A中的阵列方向的MTF参考值MTFx-ref进行对比,得到光学遥感载荷在轨MTF精度。
本步骤中,由于载荷阵列方向和飞行方向精度计算方法相同,为描述方便,以下将MTF计算值(MTFx-cal、MTFy-cal)统一表示为MTFcal,MTF参考值(MTFx-ref、MTFy-ref)统一表示为MTFref。
该获取MTF测量精度的过程又可以分为以下两种情况:
一、不考虑随机因素的情况
如果不考虑测量中的随机因素(噪声和采样相位影响),基于步骤G得到的无噪声退化参考图像g(x,y),和步骤A得到MTF参考值MTFref,得到MTF测量的系统误差。
即基于步骤G得到的退化参考图像g(x,y),按照图1所示的刃边法MTF数据处理流程计算MTF值得到MTFcal,然后以步骤A得到的MTF参考值MTFref为基准,按照公式(13)计算相对误差,即为系统误差。
二、考虑随机因素的情况
考虑测量中的随机噪声和采样相位等随机因素,基于步骤J得到的K幅带噪声刃边参考图像gk(x,y),和步骤A得到MTF参考值MTFref,采用统计方法计算MTF整体测量精度。
即对步骤J生成的K幅考虑随机噪声和采样相位的刃边随机参考图像gk(x,y),k=1,2,…K,分别利用图1所示的刃边法计算MTF,并提取出Nyquist频率下的MTF值(MTFk-cal,k=1,2,…,K,K常取100或1000),与步骤A的MTFref进行比较,按照下式得到:相对误差、平均相对误差和误差标准差。
相对误差:
平均相对误差:
均方根误差:
通常情况下,采用为在当前MTF测量下的平均测量误差,以σΔ作为MTF测量精度的1σ精度。
由于在MTF测量中存在随机因素,通常还会按照下述方式来表示MTF测量的置信度:以σΔ作为MTF测量精度的1σ精度,则MTF测量精度在区间的置信度为66.7%;以2σΔ作为MTF测量精度的2σ精度,则表示测量精度在区间的置信度为95.45%。
另外,由于通常有指标要求MTF测量的精度要满足优于10%,针对这种情况,可计算相对误差的绝对值|Δn|,统计|Δn|小于10%的个数N1,计算则可认为MTF测量精度优于10%的置信度为
以上方法除了可以获取光学遥感载荷在轨MTF测量精度,还可以通过在步骤E调整MTF测量过程中的关键参数,如靶标对比度A、B、靶标尺寸、靶标倾角等,以及在步骤H和步骤J调整噪声标准差和随机相位,然后计算系统误差和整体精度,进行关键误差因素的敏感性分析,以实现算法的优化和靶标设计与布局的优化。
至此,已经结合附图对本实施例进行了详细描述。依据以上描述,本领域技术人员应当对本发明获取刃边法测量光学遥感载荷在轨MTF精度的方法有了清楚的认识。
需要说明的是,在附图或说明书正文中,未绘示或描述的实现方式,均为所属技术领域中普通技术人员所知的形式,并未进行详细说明。此外,上述对各元件和方法的定义并不仅限于实施例中提到的各种方式,本领域普通技术人员可对其进行简单地更改或替换,例如:
(1)带噪声刃边参考图像Ik(x,y)的数目可以根据需要进行调整,可以取50、100或者200,本领域技术人员可以根据需要进行选择;
(2)本文可提供包含特定值的参数的示范,但这些参数无需确切等于相应的值,而是可在可接受的误差容限或设计约束内近似于相应值;
(3)实施例中提到的方向用语,例如“上”、“下”、“前”、“后”、“左”、“右”等,仅是参考附图的方向,并非用来限制本发明的保护范围。
综上所述,本发明提供一种获取刃边法测量光学遥感载荷在轨MTF精度的方法。该方法针对现有刃边法评估光学MTF方法的技术特点和实际观测参数,提出了系统性获取其MTF测量精度(或置信度)的方法,该方法还可以用于关键技术参数的敏感性分析,从而可为刃边的靶标设计与布设以及实施刃边法的MTF在轨评估提供理论指导和依据。
以上所述的具体实施例,对本发明的目的、技术方案和有益效果进行了进一步详细说明,所应理解的是,以上所述仅为本发明的具体实施例而已,并不用于限制本发明,凡在本发明的精神和原则之内,所做的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (14)
1.一种获取基于刃边法测量光学遥感载荷在轨MTF在第一方向的测量精度的方法,其特征在于,包括:
步骤A:在光学遥感载荷获取的实测刃边靶标图像中截取尺寸为m行n列的多个刃边子图像,由该多个刃边子图像获取光学遥感载荷在第一方向的在轨MTF参考值和刃边倾角参考值;
步骤B:由光学遥感载荷在第一方向的在轨MTF参考值,获取光学遥感载荷在该第一方向的扩散尺度,构建亚像元尺度下的参考图像退化模板psf(x,y);
步骤C:由光学遥感载荷获取的实测刃边靶标图像提取被刃边分割的两个均匀区域的灰度均值和随机噪声标准差;
步骤D:根据两个均匀区域的随机噪声标准差,获得用于获取后续步骤中亚像元水平下刃边噪声参考图像的噪声参数:刃边噪声参考图像中两个区域的噪声标准差;
步骤E:利用刃边靶标参数,构建亚像元水平的第一方向的理想刃边参考图像I(x,y),其中,所述刃边靶标参数包括:刃边子图像尺寸、靶标图像被刃边分割的两个均匀区域的灰度均值、光学遥感载荷在第一方向的刃边倾角参考值;
步骤F:将所述理想刃边参考图像I(x,y)与参考图像退化模板psf(x,y)进行二维卷积,截取卷积结果矩阵中有效数据,得到亚像元水平的第一方向的无噪声刃边退化图像Ipsf(x,y);
步骤G:对所述刃边退化图像Ipsf(x,y)中的数据进行数据聚合,生成整像元水平的第一方向的无噪声刃边退化参考图像g(x,y);
步骤H:利用所述刃边噪声参考图像中两个区域的噪声标准差,构建K幅亚像元水平的第一方向的带噪声刃边参考图像Ik(x,y),k=1,2,…K;
步骤I:将带噪声刃边参考图像Ik(x,y)k=1,2,…K,分别与参考图像退化模板psf(x,y)进行二维卷积,截取卷积结果矩阵中的有效数据得到K幅亚像元水平的第一方向的带噪声刃边退化图像Ipsf-k(x,y),k=1,2,…K;
步骤J:对所述带噪声刃边退化图像Ipsf-k(x,y),k=1,2,…K,进行数据聚合处理,生成K幅整像元水平的阵列方向带噪声退化参考图像gk(x,y),k=1,2,…K;以及
步骤K:利用刃边法MTF测试算法,对所述整像元水平的第一方向的无噪声刃边退化参考图像g(x,y),以及所述K幅整像元水平的阵列方向带噪声刃边退化参考图像gk(x,y),k=1,2,…K,分别进行第一方向的MTF计算,将Nyquist频率处的归一化MTF计算值与第一方向的在轨MTF参考值进行对比,得到光学遥感载荷在轨MTF精度;
其中,所述第一方向为载荷阵列方向和载荷运动方向其中之一。
2.根据权利要求1所述的方法,其特征在于,所述步骤A包括:
子步骤A1:在光学遥感载荷获取的实测刃边靶标图像中截取m行×n列的多个刃边子图像D,其中,在选取刃边子图像D时刃边位于刃边子图像D中间位置;
子步骤A2:对于多个刃边子图像中的每一个,采用刃边法计算其在第一方向的归一化Nyquist频率下的MTF和刃边倾角;
子步骤A3:对多个刃边子图像的在第一方向的Nyquist频率下的MTF和刃边倾角分别求平均,获取光学遥感载荷在第一方向的在轨MTF参考值和刃边倾角参考值。
3.根据权利要求1所述的方法,其特征在于,所述步骤B中,所述参考图像退化模板为:
<mrow>
<mi>p</mi>
<mi>s</mi>
<mi>f</mi>
<mrow>
<mo>(</mo>
<mi>x</mi>
<mo>,</mo>
<mi>y</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<msup>
<mi>e</mi>
<mrow>
<mo>-</mo>
<mfrac>
<msup>
<mi>x</mi>
<mn>2</mn>
</msup>
<mrow>
<mn>2</mn>
<msubsup>
<mi>&sigma;</mi>
<mi>x</mi>
<mn>2</mn>
</msubsup>
</mrow>
</mfrac>
</mrow>
</msup>
<msup>
<mi>e</mi>
<mrow>
<mo>-</mo>
<mfrac>
<msup>
<mi>y</mi>
<mn>2</mn>
</msup>
<mrow>
<mn>2</mn>
<msubsup>
<mi>&sigma;</mi>
<mi>y</mi>
<mn>2</mn>
</msubsup>
</mrow>
</mfrac>
</mrow>
</msup>
</mrow>
其中,x,y分别表示光学遥感载荷在阵列方向和运动方向的坐标,σx、σy分别表示载荷系统函数PSF在阵列方向和运动方向的扩散尺度,其满足:
<mrow>
<msub>
<mi>MTF</mi>
<mrow>
<mi>x</mi>
<mo>-</mo>
<mi>r</mi>
<mi>e</mi>
<mi>f</mi>
</mrow>
</msub>
<mo>=</mo>
<msup>
<mi>e</mi>
<mrow>
<mo>-</mo>
<mn>2</mn>
<msup>
<mi>&pi;</mi>
<mn>2</mn>
</msup>
<msup>
<msub>
<mi>&sigma;</mi>
<mi>x</mi>
</msub>
<mn>2</mn>
</msup>
<mo>/</mo>
<mn>4</mn>
</mrow>
</msup>
<mo>&times;</mo>
<mfrac>
<mn>2</mn>
<mi>&pi;</mi>
</mfrac>
</mrow>
<mrow>
<msub>
<mi>MTF</mi>
<mrow>
<mi>y</mi>
<mo>-</mo>
<mi>r</mi>
<mi>e</mi>
<mi>f</mi>
</mrow>
</msub>
<mo>=</mo>
<msup>
<mi>e</mi>
<mrow>
<mo>-</mo>
<mn>2</mn>
<msup>
<mi>&pi;</mi>
<mn>2</mn>
</msup>
<msup>
<msub>
<mi>&sigma;</mi>
<mi>y</mi>
</msub>
<mn>2</mn>
</msup>
<mo>/</mo>
<mn>4</mn>
</mrow>
</msup>
<mo>&times;</mo>
<mfrac>
<mn>2</mn>
<mi>&pi;</mi>
</mfrac>
</mrow>
其中,MTFx-ref为光学遥感载荷在载荷阵列方向的在轨MTF参考值,MTFy-ref为光学遥感载荷在载荷运动方向的在轨MTF参考值;
其中,所述MTFx-ref和MTFy-ref的其中之一通过步骤A得到,其中另一通过步骤A得到或令其与所述其中之一相等。
4.根据权利要求1所述的方法,其特征在于,所述步骤C中:
当第一方向为载荷阵列方向时,实测刃边靶标图像被刃边分割成左右两个均匀区域;
当第一方向为载荷运动方向时,实测刃边靶标图像被刃边分割成上下两个均匀区域。
5.根据权利要求4所述的方法,其特征在于,所述第一方向为载荷阵列方向,所述步骤D包括:
子步骤D1:构建一个均值为0、标准差σINZL的初始值为σnoise-L、尺寸为[(mZL+10)×10行]×[(nZL+10)×10列]的亚像元水平的噪声矩阵INZL(x,y),该噪声矩阵INZL(x,y)对应的数据点间距为0.1像元;
其中,σnoise-L为刃边靶标图像中刃边左侧均匀区域ZL的随机噪声标准差;mZL和nZL分别为刃边左侧均匀区域ZL的行数和列数;
子步骤D2:令噪声矩阵INZL(x,y)与参考图像退化模板psf(x,y)进行卷积,截取卷积结果矩阵中的有效数据((mZL×10+1)行×(nZL×10+1)列)的数据IPNZL;
子步骤D3:对矩阵IPNZL中的数据每10×10个数据为一组进行聚合,得到一个mZL行×nZL列的矩阵PNZL,PNZL数据点间距为1像元;
其中,mZL和nZL分别为刃边左侧均匀区域ZL的行数和列数;
子步骤D4:计算PNZL的噪声标准差σPNZL;
子步骤D5:对比σPNZL和σnoise-L,如果σPNZL不等于σnoise-L,则调整子步骤D1中噪声矩阵INZL(x,y)的标准差σINZL,然后重复子步骤D2~D5,直到σPNZL=σnoise-L,则将此时噪声矩阵INZL(x,y)的标准差σINZL作为后续构建亚像元水平下刃边参考图像刃边左区域噪声标准差σnL;
其中,按照子步骤D1~D5得到噪声矩阵INZR(x,y)的标准差σINZR将作为后续构建亚像元水平下刃边参考图像刃边右区域噪声标准差σnR。
6.根据权利要求4所述的方法,其特征在于,所述第一方向为载荷阵列方向,所述步骤E包括:
子步骤E1:基于刃边子图像的尺寸,构建一(m+1+10)×10行,(n+1+10)×10列的亚像元水平的参考图像P,参考图像P的数据点间距为0.1像元;
其中,刃边子图像的数据点间隔为1像元;
子步骤E2,在参考图像P中,令与图像P的列方向成θx-ref角度的一条直线l,从参考图像P中第像素的左上角点O点斜穿而过,则该直线l穿越P中某一像素时,分割该像素为左右两个部分,设参考图像P中各个像素的面积为1,根据P中各个像素被直线l分割的左右两部分的面积为权重,对参考图像P进行灰度赋值,灰度赋值后的参考图像P即为亚像元水平的载荷阵列方向的理想刃边参考图像I(x,y);
其中,θx-ref为光学遥感载荷在载荷阵列方向的刃边倾角参考值。
7.根据权利要求6所述的方法,其特征在于,所述子步骤E中对参考图像P进行灰度赋值包括:
(1)对于参考图像P中没有被直线l穿过的像素p(x,y),如果p(x,y)位于直线l左侧,则p的灰度I(x,y)赋值为:
I(x,y)=μL
(2)对于参考图像P中没有被直线l穿过的像素p(x,y),如果p(x,y)位于直线l右侧,则该p(x,y)的灰度I(x,y)赋值为:
I(x,y)=μR
(3)对于被直线l穿过的像素p(x,y),则计算直线l穿过p时将p分割成为左右两部分,所占面积分别为S1、S2,则该像素p(x,y)的灰度I(x,y)赋值为:
I(x,y)=S1·μL+S2·μR
其中,μL、μR分别为靶标图像中被刃边分割的左右两侧均匀区域的灰度均值。
8.根据权利要求1所述的方法,其特征在于,所述步骤F中无噪声刃边退化图像Ipsf(x,y)的尺寸为:(m+1)×10+1行、(n+1)×10+1列。
9.根据权利要求1所述的方法,其特征在于,所述步骤G中,对刃边退化图像Ipsf(x,y),以Ipsf(x,y)中位于行、列的像素为起始位置,在图像Ipsf(x,y)中上下左右每10×10个像素为1组进行数据聚合,生成整像元水平的第一方向的无噪声刃边退化参考图像g(x,y)。
10.根据权利要求1所述的方法,其特征在于,所述第一方向为载荷阵列方向,所述步骤H包括:
子步骤H1:构建K幅(m+1+10)×10行、(n+1+10)×10列的亚像元水平的参考图像Pk,k=1,2,…K,参考图像Pk的数据点间距为0.1像元;
子步骤H2:分别构建K个均值为0、标准差σnL、尺寸为(m+1+10)×10行、(n+1+10)×10列的噪声矩阵NPLk,k=1,2,…K,和K个均值为0、标准差σnR、尺寸为(m+1+10)×10行、(n+1+10)×10列的噪声矩阵NPRk,k=1,2,…K,噪声矩阵NPLk和NPRk,k=1,2,…K的数据点间距为0.1像元;
其中,σnL和σnR分别为亚像元水平下刃边参考图像刃边左区域和右区域的噪声标准差;
子步骤H3,在参考图像Pk,k=1,2,…K中,令与参考图像Pk的列方向成θx-ref角度的一条直线l从参考图像Pk中第像素的左上角点O点斜穿而过,则根据该直线l穿越Pk中某一像素时,分割该像素为左右两个部分,设参考图像Pk中各个像素的面积为1,根据Pk中各个像素被直线l分割的左右两部分的面积为权重,对参考图像Pk进行灰度赋值,生成K幅亚像元水平的阵列方向带随机噪声的刃边参考图像Ik(x,y),k=1,2,…K;
其中,所述θx-ref为光学遥感载荷在载荷阵列方向的刃边倾角参考值。
11.根据权利要求10所述的方法,其特征在于,所述子步骤H2中,对参考图像Pk,k=1,2,…K,进行灰度赋值包括:
(1)对于Pk中没有被直线l穿过的像素pk(x,y),如果pk(x,y)位于直线l左侧,则pk(x,y)的灰度Ik(x,y)赋值为:
Ik(x,y)=μL+NPLk(x,y);
(2)如果pk(x,y)位于直线l右侧,则该pk(x,y)的灰度Ik(x,y)赋值为:
Ik(x,y)=μR+NPRk(x,y);
(3)对于被直线l穿过的像素pk(x,y),则计算直线l穿过pk时将pk分割成为左右两部分,所占面积分别为S1k、S2k,则该像素pk(x,y)的灰度Ik(x,y)赋值为:
Ik(x,y)=S1k·[μL+NPLk(x,y)]+S2k·[μR+NPRk(x,y)]。
12.根据权利要求1所述的方法,其特征在于,所述步骤J包括:
子步骤J1:在(-5、-4、-3、-2、-1、0、1、2、3、4、5)这11个数中每次随机选出2个数,共随机选出K对数值(Sxk,Syk)k=1,2,…K,作为随机相位的偏移值;
子步骤J2:对第k幅亚像元水平的阵列方向带噪声刃边退化图像Ipsf-k(x,y),k=1,2,…K,以Ipsf-k(x,y)中第行、列的像素为聚合起始位置,在图像Ipsf-k(x,y)中上下左右每10×10个像素为1组进行聚合,得到K幅带有随机噪声和随机相位的、间距为1像元的整像元水平的阵列方向带噪声退化刃边参考图像gk(x,y),k=1,2,…K。
13.根据权利要求1至12中任一项所述的方法,其特征在于,在考虑随机因素时,以MTF的均方根误差作为光学遥感载荷在轨MTF的1σ精度,以MTF测量精度在区间作为光学遥感载荷在轨MTF的1σ置信度,其中:
系统的平均相对误差:
<mfenced open = "" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<mover>
<mi>&Delta;</mi>
<mo>&OverBar;</mo>
</mover>
<mo>=</mo>
<mi>a</mi>
<mi>v</mi>
<mi>g</mi>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<msub>
<mi>MTF</mi>
<mrow>
<mi>k</mi>
<mo>-</mo>
<mi>c</mi>
<mi>a</mi>
<mi>l</mi>
</mrow>
</msub>
<mo>-</mo>
<msub>
<mi>MTF</mi>
<mrow>
<mi>r</mi>
<mi>e</mi>
<mi>f</mi>
</mrow>
</msub>
</mrow>
<mrow>
<msub>
<mi>MTF</mi>
<mrow>
<mi>r</mi>
<mi>e</mi>
<mi>f</mi>
</mrow>
</msub>
</mrow>
</mfrac>
<mo>&times;</mo>
<mn>100</mn>
<mi>%</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mi>N</mi>
</mfrac>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>k</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>K</mi>
</munderover>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<msub>
<mi>MTF</mi>
<mrow>
<mi>k</mi>
<mo>-</mo>
<mi>c</mi>
<mi>a</mi>
<mi>l</mi>
</mrow>
</msub>
<mo>-</mo>
<msub>
<mi>MTF</mi>
<mrow>
<mi>r</mi>
<mi>e</mi>
<mi>f</mi>
</mrow>
</msub>
</mrow>
<mrow>
<msub>
<mi>MTF</mi>
<mrow>
<mi>r</mi>
<mi>e</mi>
<mi>f</mi>
</mrow>
</msub>
</mrow>
</mfrac>
<mo>)</mo>
</mrow>
<mo>&times;</mo>
<mn>100</mn>
<mi>%</mi>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
均方根误差:
<mfenced open = "" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mi>&sigma;</mi>
<mi>&Delta;</mi>
</msub>
<mo>=</mo>
<mi>s</mi>
<mi>t</mi>
<mi>d</mi>
<mrow>
<mo>(</mo>
<mfrac>
<mrow>
<msub>
<mi>MTF</mi>
<mrow>
<mi>k</mi>
<mo>-</mo>
<mi>c</mi>
<mi>a</mi>
<mi>l</mi>
</mrow>
</msub>
<mo>-</mo>
<msub>
<mi>MTF</mi>
<mrow>
<mi>r</mi>
<mi>e</mi>
<mi>f</mi>
</mrow>
</msub>
</mrow>
<mrow>
<msub>
<mi>MTF</mi>
<mrow>
<mi>r</mi>
<mi>e</mi>
<mi>f</mi>
</mrow>
</msub>
</mrow>
</mfrac>
<mo>&times;</mo>
<mn>100</mn>
<mi>%</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>=</mo>
<msqrt>
<mfrac>
<mrow>
<munderover>
<mo>&Sigma;</mo>
<mrow>
<mi>k</mi>
<mo>=</mo>
<mn>1</mn>
</mrow>
<mi>K</mi>
</munderover>
<msup>
<mrow>
<mo>(</mo>
<msub>
<mi>MTF</mi>
<mrow>
<mi>k</mi>
<mo>-</mo>
<mi>c</mi>
<mi>a</mi>
<mi>l</mi>
</mrow>
</msub>
<mo>-</mo>
<mover>
<mi>&Delta;</mi>
<mo>&OverBar;</mo>
</mover>
<mo>)</mo>
</mrow>
<mn>2</mn>
</msup>
</mrow>
<mrow>
<mi>K</mi>
<mo>-</mo>
<mn>1</mn>
</mrow>
</mfrac>
</msqrt>
<mo>&times;</mo>
<mn>100</mn>
<mi>%</mi>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
其中,MTFk-cal为刃边随机参考图像gk(x,y)利用刃边法计算得到的MTF值,并提取的Nyquist频率下的MTF值,MTFref为光学遥感载荷在第一方向的在轨MTF参考值。
14.根据权利要求1至12中任一项所述的方法,所述步骤K中,在不考虑随机因素时,系统误差为:
<mrow>
<mi>&Delta;</mi>
<mo>=</mo>
<mfrac>
<mrow>
<msub>
<mi>MTF</mi>
<mrow>
<mi>c</mi>
<mi>a</mi>
<mi>l</mi>
</mrow>
</msub>
<mo>-</mo>
<msub>
<mi>MTF</mi>
<mrow>
<mi>r</mi>
<mi>e</mi>
<mi>f</mi>
</mrow>
</msub>
</mrow>
<mrow>
<msub>
<mi>MTF</mi>
<mrow>
<mi>r</mi>
<mi>e</mi>
<mi>f</mi>
</mrow>
</msub>
</mrow>
</mfrac>
<mo>&times;</mo>
<mn>100</mn>
<mi>%</mi>
<mo>,</mo>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>=</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>,</mo>
<mo>...</mo>
<mo>,</mo>
<mi>K</mi>
<mo>)</mo>
</mrow>
</mrow>
其中,MTFcal为基于无噪声退化参考图像g(x,y)利用刃边法计算得到的MTF值,MTFref为光学遥感载荷在第一方向的在轨MTF参考值。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510711683.6A CN105427267B (zh) | 2015-10-28 | 2015-10-28 | 获取刃边法测量光学遥感载荷在轨mtf精度的方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201510711683.6A CN105427267B (zh) | 2015-10-28 | 2015-10-28 | 获取刃边法测量光学遥感载荷在轨mtf精度的方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN105427267A CN105427267A (zh) | 2016-03-23 |
CN105427267B true CN105427267B (zh) | 2018-03-23 |
Family
ID=55505448
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201510711683.6A Active CN105427267B (zh) | 2015-10-28 | 2015-10-28 | 获取刃边法测量光学遥感载荷在轨mtf精度的方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN105427267B (zh) |
Families Citing this family (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105976317B (zh) * | 2016-04-28 | 2019-03-19 | 中国科学院遥感与数字地球研究所 | 一种图像空间退化模拟方法及系统 |
CN117152069B (zh) * | 2023-08-17 | 2024-08-20 | 哈尔滨工业大学 | 一种空间目标高精度指向与成像评估试验装置及其试验方法 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104298844A (zh) * | 2014-05-23 | 2015-01-21 | 中国科学院光电研究院 | 获取点阵法测量光学遥感载荷在轨mtf测量精度的方法 |
CN104434150A (zh) * | 2013-09-18 | 2015-03-25 | 中国科学院深圳先进技术研究院 | 数字x线成像系统的二维调制传递函数测量方法及系统 |
CN104833480A (zh) * | 2015-04-27 | 2015-08-12 | 北京空间机电研究所 | 一种测试刃边倾角对调制传递函数影响的装置及方法 |
-
2015
- 2015-10-28 CN CN201510711683.6A patent/CN105427267B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104434150A (zh) * | 2013-09-18 | 2015-03-25 | 中国科学院深圳先进技术研究院 | 数字x线成像系统的二维调制传递函数测量方法及系统 |
CN104298844A (zh) * | 2014-05-23 | 2015-01-21 | 中国科学院光电研究院 | 获取点阵法测量光学遥感载荷在轨mtf测量精度的方法 |
CN104833480A (zh) * | 2015-04-27 | 2015-08-12 | 北京空间机电研究所 | 一种测试刃边倾角对调制传递函数影响的装置及方法 |
Non-Patent Citations (3)
Title |
---|
MTF assessment of high resoIution satellite images using ISo 12233 slanted-edge method;Hyundeok Hwang 等;《SPIE》;20081001;710905-1-710905-9 * |
一种优化的刃边法MTF在轨评估算法;徐航 等;《遥感信息》;20121215;第27卷(第6期);10-16 * |
高精度刃边法测量红外相机MTF的研究;胡涛 等;《红外》;20150710;第36卷(第7期);10-15 * |
Also Published As
Publication number | Publication date |
---|---|
CN105427267A (zh) | 2016-03-23 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Giannantonio et al. | CMB lensing tomography with the DES Science Verification galaxies | |
Guo et al. | FSDAF 2.0: Improving the performance of retrieving land cover changes and preserving spatial details | |
Wu et al. | An error-bound-regularized sparse coding for spatiotemporal reflectance fusion | |
CN103793907B (zh) | 一种水体信息的提取方法及装置 | |
Canty et al. | Automatic radiometric normalization of multitemporal satellite imagery | |
CN101465002B (zh) | 椭圆目标的亚像素边缘定位方法 | |
Gray et al. | STAGES: the Space Telescope A901/2 Galaxy Evolution Survey | |
CN101980293A (zh) | 一种基于刃边图像的高光谱遥感系统mtf检测方法 | |
Theiler et al. | Local coregistration adjustment for anomalous change detection | |
CN105930772A (zh) | 基于sar影像与光学遥感影像融合的城市不透水面提取方法 | |
GB2522365B (en) | Method for processing an image | |
Schlapfer et al. | Spatial PSF nonuniformity effects in airborne pushbroom imaging spectrometry data | |
CN103413292B (zh) | 基于约束最小二乘的高光谱图像非线性丰度估计方法 | |
Tang et al. | A multiple-point spatially weighted k-NN method for object-based classification | |
CN104615977A (zh) | 综合关键季相特征和模糊分类技术的冬小麦遥感识别方法 | |
CN106447653B (zh) | 基于空间约束的卡方变换的多时相遥感影像变化检测方法 | |
CN110388986B (zh) | 基于tasi数据的土地表面温度反演方法 | |
CN103776532A (zh) | 一种基于遥感应用的高光谱成像仪指标优化方法 | |
CN104298844B (zh) | 获取点阵法测量光学遥感载荷在轨mtf测量精度的方法 | |
Zhu et al. | Potential and limits of non-local means InSAR filtering for TanDEM-X high-resolution DEM generation | |
CN105427267B (zh) | 获取刃边法测量光学遥感载荷在轨mtf精度的方法 | |
CN117409335A (zh) | 一种基于可见光图像的气象雷达降水率降尺度方法 | |
CN109726679B (zh) | 一种遥感分类误差空间分布制图方法 | |
CN106960443A (zh) | 基于全极化时序sar图像的非监督变化检测的方法及装置 | |
Refice et al. | On the use of anisotropic covariance models in estimating atmospheric DInSAR contributions |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |