基于曲线拟合的大气参量查找表生成方法与流程

本发明属于物理技术领域，更进一步涉及物理技术领域中的一种基于曲线拟合的大气参量查找表生成方法。本发明可生成光电三维场景中的大气参量查找表，用于大规模光电三维场景仿真的实时渲染。

背景技术：

随着计算机技术和仿真技术的飞速发展，大规模光电三维场景的实时仿真一直是研究的热点。光电三维场景的实时仿真包括对大气传输效应的实时仿真，大气传输效应的实时仿真需要根据仿真帧频和分辨率计算仿真图像中的每个像素所对应的辐射传输几何路径下的包括大气透过率和大气路径辐射的大气参量，计算量庞大。针对上述计算量庞大的问题，通过对大气透过率和大气路径辐射进行预计算，生成大气参量查找表是一种有效的技术方法。

中国航天科工集团第二研究院二ο七所在其申请的专利文献“红外特性测量中大气传输快速修正方法”(专利申请号201210362216.3，申请号公布号cn103674904a)中公开了一种红外特性测量中大气传输快速修正方法。该方法的主要步骤为：第一，设置计算参数；第二，选择气象模式；第三，根据所选的气象模式调用modtran软件中相应的计算模型来进行大气传输参数计算；第四，计算距离测量位置不同观测高度和水平投影距离位置点的大气透过率和背景辐射亮度；第五，调用计算出的大气透过率和背景辐射亮度网格数据表格，对数据进行查分，得出不同观测仰角和观测距离点的大气透过率和背景辐射亮度变化的曲线。该方法存在的不足之处是：该方法辐射传输几何参数计算间隔太小，计算的数据量太大，计算耗时。

孙明明在其发表论文“红外场景的大气效应高效计算与仿真研究”(西安电子科技大学硕士论文2014年)中提出了一种计算大气参量的方法。该方法的主要步骤为：第一，设置气象参数；第二，modtran预计算所选辐射传输几何参数下的大气透过率和大气路径辐射；第三，选取拟合函数形式；第四，采用麦夸特算法拟合大气透过率和大气路径辐射随路径长度变化的曲线；第五，对路径长度进行密集取值，生成仿真场景的大气参量查找表。该方法存在的不足之处是：大气参量的计算精度低，计算耗时，计算量大，计算场景中所有辐射传输几何路径的大气参量时，需遍历所有探测器海拔高度值。

技术实现要素：

本发明的目的在于克服上述现有技术存在的不足，提出一种基于曲线拟合的大气参量查找表生成方法，本发明可以显著减小生成大规模光电三维场景的大气参量查找表的计算量，提高生成大气参量查找表的速度，计算精度高，方法简单实用，能满足大规模光电三维场景仿真的实时性要求。

为实现上述目的，本发明的具体步骤如下：

(1)设置曲线拟合仿真参数：

(1a)设置曲线拟合的气象条件参数；

(1b)设置曲线拟合仿真波段；

(1c)将与大气参量关联的探测器海拔高度、视线天顶角、路径长度设置为modtran软件中的辐射传输几何路径参数；

(1d)将待仿真场景的几何尺度作为探测器海拔高度的最大值、视线天顶角的最大值和路径长度的最大值；

(2)将探测器视线天顶角的初始值赋值为0；

(3)划分探测器海拔高度子区间：

在探测器海拔高度取值范围[0.2，15.4]的区间内，将探测器海拔高度的划分为[0.2，4]、(4，7.8]、(7.8，11.6]、(11.6，15.4]四个子区间；

(4)按照下式，分别计算四个高度分段子区间中的每个探测器海拔高度值：

其中，hjm表示在第m个探测器海拔高度子区间内的第j个探测器海拔高度值，j在1≤j≤12中取整数值，hminm表示探测器海拔高度中第m个子区间的左端点；

(5)计算路径长度最大有效值：

(5a)按照下式，计算视线天顶角大于90度时，辐射在穿过大气的传输路径上不与地球表面相交的路径长度的最大有效值：

其中，l'max表示辐射在穿过大气的传输路径上不与地球表面相交的路径长度的最大有效值，h表示探测器海拔高度，-表示负操作，θ表示视线天顶角，cos表示余弦操作；

(5b)视线天顶角不大于90度时，令辐射在穿过大气的传输路径上不与地球表面相交的路径长度最大有效值等于路径长度最大值；

(6)划分路径长度子区间：

(6a)判断视线天顶角是否大于90度；若是，则执行步骤(6b)；否则，执行步骤(7)；

(6b)按照下式，计算路径长度子区间分段点处的路径长度值：

其中，lpart表示路径长度子区间分段点处的路径长度值；e表示自然常数，其值为2.71828；

(6c)在路径长度的取值范围[0,l'max]的区间内，将路径长度划分为

[0,lpart]和(lpart,l'max]两个路径长度子区间；

(7)计算路径长度值：

(7a)按照下式，计算区间[0,lpart]的路径长度的起始长度控制因子：

其中，pstart表示区间[0,lpart]的路径长度子的起始路径长度控制因子；

(7b)按照下式，计算在区间[0,lpart]中选取的路径长度控制因子的个数：

其中，q表示在区间[0,lpart]中选取的路径长度因子个数，[]表示向下取整；

(7c)按照下式，计算视线天顶角大于90度，区间[0,lpart]中选取的每个路径长度控制因子中的第q个路径长度控制因子：

其中，pq表示视线天顶角大于90度，在路径长度子区间[0,lpart]中选取的q个路径长度控制因子中的第q个路径长度控制因子，pq中的q在1≤q≤q中取整数值，pq-1表示在区间[0,lpart]中总共选取的q个路径长度控制因子中的第q-1个路径长度控制因子，pq-1中的q在1＜q≤q中取整数值；

(7d)按照下式，计算视线天顶角大于90度时，区间(lpart,l'max]的终止路径长度控制因子：

其中，pend表示路径长度子区间(lpart,l'max]的终止路径长度控制因子；

(7e)按照下式，计算视线天顶角大于90度时，区间[0,lpart]中的每个路径长度值：

其中，lq表示在区间[0,lpart]中所选取的每个路径长度值中的第q个路径长度值，q在1≤q≤q取整数值；

(7f)按照下式，计算视线天顶角大于90度，在区间(lpart,l'max]中选取的每个路径长度控制因子中的第r个路径长度控制因子：

其中，pr表示视线天顶角大于90度，在路径长度子区间(lpart,l'max]中选取的每个路径长度控制因子中的第r个路径长度控制因子，pr中的r在1≤r≤r中取整数值，pr-1表示视线天顶角大于90度，在区间(lpart,l'max]中选取的每个路径长度控制因子中的第r-1个路径长度控制因子，pr-1中的r在1＜r≤r中取整数值，r表示在区间(lpart,l'max]中选取的路径长度控制因子的个数，当h＜2.6km时，r等于7，当h≥2.6km时，r等于9；

(7g)按照下式，分别计算视线天顶角大于90度时区间(lpart,l'max]中的每个路径长度值：

其中，lr表示在区间(lpart,l'max]中所选取的每个路径长度值中的第r个路径长度值，r在1≤r≤r取整数值；

(7h)按照下式，计算视线天顶角不大于90度时路径长度的取值范围[0,l'max]中的每个路径的长度：

其中，ln表示路径长度的取值范围[0,l'max]中所选取的每个路径长度值中的第n个路径的长度值；

(8)计算大气参量：

(8a)将辐射传输几何路径参数(θ,hjm,lq)和(θ,hjm,lr)输入modtran软件中，分别计算视线天顶角大于90度时大气参量e(θ,hjm,lq)和大气参量e(θ,hjm,lr)；

(8b)将辐射传输几何路径参数(θ,hjm,ln)输入modtran软件中，计算视线天顶角不大于90度时的大气参量e(θ,hjm,ln)；

(9)计算拟合曲线参数；

(9a)根据最小二乘方法，将视线天顶角不大于90度的ln和e(θ,hjm,ln)作为原始数据，计算使得拟合数据与原始数据的误差平方和最小时的θ和hjm取确定值时的大气参量随路径长度变化的曲线的八个曲线参数，其中，表示的第s+1项参数，l表示路径长度；

(9b)根据最小二乘方法，将视线天顶角大于90度的路径长度子区间[0,lpart]中的lq和e(θ,hjm,lq)作为原始数据，计算使得拟合数据与原始数据的误差平方和最小时的θ和hjm取确定值时的大气参量随路径长度变化的曲线的七个曲线参数，其中，表示的第s+1项参数；

(9c)根据最小二乘方法，将视线天顶角大于90度的路径长度子区间(lpart,l'max]中的lr和e(θ,hjm,lr)作为原始数据，计算使得拟合数据与原始数据的误差平方和最小时的θ和hjm取确定值时的大气参量随路径长度变化的曲线的七个参数，其中，表示的第s+1项参数；

(9d)根据最小二乘方法，将视线天顶角不大于90度的hjm和作为原始数据，计算s+1条使得拟合数据与原始数据的误差平方和最小时的θ取确定值时的大气参量随路径长度变化的曲线系数随探测器海拔高度变化的曲线a(θ)s(h)的七个曲线参数，d(θ)st表示a(θ)s(h)的第t+1项参数，h表示探测器海拔高度；

(9e)按照下式，计算视线天顶角大于90度，θ和hjm取确定值，路径子区间[0,lpart]的任意高度的大气参量随路径长度变化的曲线的各项系数：

其中，b(θ)s(h)表示视线天顶角大于90度，θ和hjm取确定值，路径子区间[0,lpart]任意探测器海拔高度的大气参量随路径长度变化的曲线的第s+1项参数，表示视线天顶角为θ，探测器海拔高度值为h(j-1)m时大气参量随路径长度变化的曲线的第s+1项系数，其中，h(j-1)m表示第m个探测器海拔高度子区间中第j-1个探测器海拔高度值；

(9f)按照下式，计算视线天顶角大于90度，θ和hjm取确定值时路径子区间(lpart,l'max]的任意高度的大气参量随路径长度变化的曲线的各项系数：

其中，c(θ)s(h)表示视线天顶角大于90度，θ和hjm取确定值，探测器海拔高度等于h时路径子区间(lpart,l'max]的大气参量随路径长度变化的曲线的第s+1项参数，表示θ和h(j-1)m取确定值时大气参量随路径长度变化的曲线的第s+1项系数；

(10)对视线天顶角进行累加5操作，得到累加后视线天顶角；

(11)判断累加后视线天顶角是否大于180度，若是，则执行步骤(12)；否则，执行步骤(3)；

(12)生成大气参量查找表：

(12a)按步骤(2)和步骤(10)中的方法对视线天顶角取值，在零至探测器海拔高度最大值范围内对探测器海拔高度随机取值，在零至路径长度最大有效值范围内对路径长度随机取值；

(12b)将选取的视线天顶角值、探测器海拔高度值和路径长度值代入步骤(9)确定的曲线中，计算大气参量，生成以视线天顶角、探测器海拔高度和路径长度为索引的大气参量查找表。

本发明与现有技术相比具有以下优点：

第一，由于本发明划分了路径长度子区间，按照不同的计算方法计算了不同路径长度子区间的路径长度值，克服了现有技术中辐射传输几何参数计算间隔太小，计算的数据量太大，计算耗时的缺陷。使得本发明具有能用更少的计算量，生成网格间隔更密的大气参量查找表，计算效率高的优点。

第二，由于本发明计算了路径长度最大有效值，对路径长度在零至路径长度最大有效值范围内取值，并且计算了大气参量随路径长度变化的曲线参数和大气参量随路径长度变化的曲线参数随探测器海拔高度变化的曲线参数，克服了现有技术中大气参量的计算精度低，计算耗时，计算大，计算场景中所有辐射传输几何路径的大气参量时，需遍历所有探测器海拔高度值的缺陷，使得本发明具有计算精度高，计算量少的优点。

附图说明

图1为本发明的流程图；

图2为本发明生成的大气透过率查找表中的大气透过率值与modtran计算的标准大气透过率值的对比图；

图3为本发明生成的大气路径辐射查找表中的大气路径辐射值与modtran计算的标准大气路径辐射值的对比图。

具体实施方式

下面结合附图对本发明做进一步的详细描述。

参照附图1，本发明的具体步骤如下。

步骤1，设置曲线拟合仿真参数。

第1步，设置仿真的气象条件参数。

本发明的实施例中，气象条件参数设置如下：

大气模式为1976年美国标准大气，温度和气压的高度廓线使用模式默认值，水汽的高度廓线使用模式默认值，臭氧的高度廓线使用模式默认值，甲烷的高度廓线使用模式默认值，一氧化二氮的高度廓线使用模式默认值，一氧化碳的高度廓线使用模式默认值，二氧化碳混合率330，气溶胶模式为乡村型能见度23km，气溶胶的季节性修正为由模式决定的修正，高层大气气溶胶为同温层背景气溶胶。

第2步，设置曲线拟合的仿真波段。

本发明的实施例中，仿真波段为3μm～5μm。

第3步，将与大气传输效应物理量关联的探测器海拔高度、视线天顶角、路径长度设置为modtran软件中的辐射传输几何路径参数。

第4步，将待仿真场景的几何尺度作为探测器海拔高度的最大值、视线天顶角的最大值和路径长度的最大值。

本发明的实施例中，仿真场景的几何尺度为70km*70km，探测器海拔高度的最大值为15.4km，视线天顶角的最大值为180度，路径长度最大值为89.32km。

步骤2，将探测器视线天顶角的初始值赋值为0。

步骤3，划分探测器海拔高度子区间。

在探测器海拔高度取值范围[0.2，15.4]的区间内，将探测器海拔高度的划分为[0.2，4]、(4，7.8]、(7.8，11.6]、(11.6，15.4]四个子区间；

步骤4，确定探测器海拔高度值。

按照下式，分别对四个高度分段子区间中的每个探测器海拔高度取值：

其中，hjm表示在第m个探测器海拔高度子区间内的第j个探测器海拔高度值；hminm表示探测器海拔高度中第m个子区间的左端点；j在1≤j≤12中取整数值。

步骤5，获得路径长度最大有效值。

按照下式，计算视线天顶角大于90度时辐射在传输路径上不与地球表面相交的路径长度最大有效值：

其中，l'max表示辐射在穿过大气的传输路径上不与地球表面相交的路径长度的最大有效值，θ表示视线天顶角，cos表示余弦操作；h表示探测器海拔高度。

视线天顶角不大于90度时，令辐射在传输路径上不与地球表面相交的路径长度最大有效值等于路径长度最大值。

步骤6，划分路径长度子区间。

第1步，判断视线天顶角是否大于90度，若是，则执行步骤6中的第2步；否则，执行步骤7。

第2步，按照下式，确定路径长度子区间分段点处的路径长度值，得到划分路径长度后的路径长度子区间：

其中，lpart表示路径长度子区间分段点处的路径长度值；e表示自然常数，其值为2.71828。

第3步，在路径长度的取值范围[0,l'max]的区间内，将路径长度划分为[0,lpart]和(lpart,l'max]两个路径长度子区间。

步骤7，计算路径长度值。

按照下式，计算区间[0,lpart]的路径长度的起始长度控制因子：

其中，pstart表示区间[0,lpart]的路径长度子的起始路径长度控制因子；

按照下式，计算在区间[0,lpart]中选取的路径长度控制因子的个数：

其中，q表示在区间[0,lpart]中选取的路径长度因子个数，[]表示向下取整；

按照下式，计算视线天顶角大于90度，区间[0,lpart]中选取的每个路径长度控制因子中的第q个路径长度控制因子：

其中，pq表示视线天顶角大于90度，在路径长度子区间[0,lpart]中选取的q个路径长度控制因子中的第q个路径长度控制因子，pq中的q在1≤q≤q中取整数值，pq-1表示在区间[0,lpart]中总共选取的q个路径长度控制因子中的第q-1个路径长度控制因子，pq-1中的q在1＜q≤q中取整数值；

按照下式，计算视线天顶角大于90度时，区间(lpart,l'max]的终止路径长度控制因子：

其中，pend表示路径长度子区间(lpart,l'max]的终止路径长度控制因子

按照下式，计算视线天顶角大于90度时，区间[0,lpart]中的每个路径长度值：

其中，lq表示在区间[0,lpart]中所选取的每个路径长度值中的第q个路径长度值，q在1≤q≤q取整数值；

按照下式，计算视线天顶角大于90度，在区间(lpart,l'max]中选取的每个路径长度控制因子中的第r个路径长度控制因子：

其中，pr表示视线天顶角大于90度，在路径长度子区间(lpart,l'max]中选取的每个路径长度控制因子中的第r个路径长度控制因子，pr中的r在1≤r≤r中取整数值，pr-1表示视线天顶角大于90度，在区间(lpart,l'max]中选取的每个路径长度控制因子中的第r-1个路径长度控制因子，pr-1中的r在1＜r≤r中取整数值，r表示在区间(lpart,l'max]中选取的路径长度控制因子的个数，当h＜2.6km时，r等于7，当h≥2.6km时，r等于9；

按照下式，分别计算视线天顶角大于90度时区间(lpart,l'max]中的每个路径长度值：

其中，lr表示在区间(lpart,l'max]中所选取的每个路径长度值中的第r个路径长度值，r在1≤r≤r取整数值；

按照下式，计算视线天顶角不大于90度时路径长度的取值范围[0,l'max]中的每个路径的长度：

其中，ln表示路径长度的取值范围[0,l'max]中所选取的每个路径长度值中的第n个路径的长度值；

步骤8，计算大气参量。

将辐射传输几何路径参数(θ,hjm,lq)和(θ,hjm,lr)输入modtran软件中，分别计算视线天顶角大于90度时大气参量e(θ,hjm,lq)和大气参量e(θ,hjm,lr)；

将辐射传输几何路径参数(θ,hjm,ln)输入modtran软件中，计算视线天顶角不大于90度时的大气传输效应物理量e(θ,hjm,ln)；

步骤9，计算拟合曲线参数。

根据最小二乘方法，将视线天顶角不大于90度的ln和e(θ,hjm,ln)作为原始数据，计算使得拟合数据与原始数据的误差平方和最小时的θ和hjm取确定值时的大气参量随路径长度变化的曲线的八个曲线参数，其中，表示的第s+1项参数，l表示路径长度。

根据最小二乘方法，将视线天顶角大于90度的路径长度子区间[0,lpart]中的lq和e(θ,hjm,lq)作为原始数据，计算使得拟合数据与原始数据的误差平方和最小时的θ和hjm取确定值时的大气参量随路径长度变化的曲线的七个曲线参数，其中，表示的第s+1项参数。

根据最小二乘方法，将视线天顶角大于90度的路径长度子区间(lpart,l'max]中的lr和e(θ,hjm,lr)作为原始数据，计算使得拟合数据与原始数据的误差平方和最小时的θ和hjm取确定值时的大气参量随路径长度变化的曲线的七个参数，其中，表示的第s+1项参数。

根据最小二乘方法，将视线天顶角不大于90度的hjm和作为原始数据，计算s+1条使得拟合数据与原始数据的误差平方和最小时的θ取确定值时的大气参量随路径长度变化的曲线系数随探测器海拔高度变化的曲线a(θ)s(h)的七个曲线参数，d(θ)st表示a(θ)s(h)的第t+1项参数，h表示探测器海拔高度。

按照下式，计算视线天顶角大于90度，θ和hjm取确定值，路径子区间

[0,lpart]的任意高度的大气参量随路径长度变化的曲线的各项系数：

其中，b(θ)s(h)表示视线天顶角大于90度，θ和hjm取确定值，路径子区间[0,lpart]任意探测器海拔高度的大气参量随路径长度变化的曲线的第s+1项参数，表示视线天顶角为θ，探测器海拔高度值为h(j-1)m时大气参量随路径长度变化的曲线的第s+1项系数，其中，h(j-1)m表示第m个探测器海拔高度子区间中第j-1个探测器海拔高度值。

按照下式，计算视线天顶角大于90度，θ和hjm取确定值时路径子区间

(lpart,l'max]的任意高度的大气参量随路径长度变化的曲线的各项系数：

其中，c(θ)s(h)表示视线天顶角大于90度，θ和hjm取确定值，探测器海拔高度等于h时路径子区间(lpart,l'max]的大气参量随路径长度变化的曲线的第s+1项参数，表示θ和h(j-1)m取确定值时大气参量随路径长度变化的曲线的第s+1项系数。

步骤10，对视线天顶角进行累加5操作，得到累加后视线天顶角。

步骤11，判断累加5后视线天顶角θ是否大于180度，若是，则执行步骤12；否则，执行步骤3。

步骤12，生成大气参量查找表。

按步骤2和步骤10中的方法对视线天顶角取值，在零至探测器海拔高度最大值范围内对探测器海拔高度随机取值，在零至路径长度最大有效值范围内对路径长度随机取值；

将选取的视线天顶角值、探测器海拔高度值和路径长度值代入步骤9确定的曲线中，计算大气参量，生成以视线天顶角、探测器海拔高度和路径长度为索引的大气参量查找表。

本发明的实施例中，对视线天顶角在[0，180]范围内，每间隔5度进行取值；对探测器海拔高度在[0.2，15.4km]范围内，每间隔0.2km进行取值；对路径长度在[0.01218，l'max]范围内，每间隔进行取值，25≤k≤114。

下面结合仿真实验对本发明的效果做进一步的描述。

1.仿真内容

为了验证基于曲线拟合的大气参量查找表的有效性，生成了以视线天顶角、探测器海拔高度和路径长度为索引的大气透过率查找表和大气路径辐射查找表。

气象条件参数设置如下：大气模式为1975年美国标准大气；天气状况为晴天；近地面大气温度为0；边界层气溶胶模式中的消光系数和气象视距类型为乡村消光系数，缺省气象视距为23km；地面能见度为23km；光谱波段为3um-5um。大气透过率和大气路径辐射查找表中视线天顶角的范围为[0，180]，间隔为5度，进行取值探测器海拔高度范围为[0.2，15.4km]，间隔为0.2km，路径长度最大值为[0.01218，l'max]，间隔为25≤k≤114。

2.仿真效果分析

图2是本发明生成的大气透过率查找表中的大气透过率值与modtran计算的标准大气透过率值的对比图。其中图2中的横坐标表示路径长度，纵坐标表示大气透过率。图2中的所有的圆形点表示视线天顶角为30度，探测器海拔高度为13.6km时modtran软件计算的标准大气透过率值。所有的十字点表示视线天顶角为30度，探测器海拔高度为13.6km时本发明计算的大气透过率值。所有的正方形点表示视线天顶角为30度，探测器海拔高度为6.4km时modtran软件计算的标准大气透过率值。所有的星形点表示视线天顶角为30度，探测器海拔高度为6.4km时本发明计算的大气透过率值。

比较图2中的数据，可以发现所有的圆形点与所有十字点几乎重合，所有的正方形点与所有的星形点几乎重合，可见本发明生成的大气参量查找表中大气透过率值与modtran软件计算的标准大气透过率值具有很高的吻合度，计算精度高。

图3是本发明生成的大气路径辐射查找表中的大气路径辐射值和modtran软件计算的标准大气路径辐射值的对比图，其中图2的横坐标表示路径长度，纵坐标表示大气路径辐射。图3中的所有的圆形点表示视线天顶角为30度，探测器海拔高度为13.6km时modtran软件计算的标准大气路径辐射值。所有的十字点表示视线天顶角为30度，探测器海拔高度为13.6km时本发明计算的大气路径辐射值。所有的正方形点表示视线天顶角为30度，探测器海拔高度为6.4km时modtran软件计算的标准大气路径辐射值。所有的星形点表示视线天顶角为30度，探测器海拔高度为6.4km时本发明计算的大气路径辐射值。

比较图3中的数据，可以发现所有的圆形点与所有十字点几乎重合，所有的正方形点与所有的星形点几乎重合，可见本发明生成的大气参量查找表中大气路径辐射值与modtran软件计算的标准大气路径辐射值具有很高的吻合度，计算精度高。