一种基于蜂群算法的遥感影像有理函数模型结构优化方法

摘要

一种基于蜂群算法的遥感影像有理函数模型优化方法，涉及遥感成像几何模型优化领域。针对基于有理函数模型的遥感成像几何模型因参数项过冗余造成的模型精度低问题，用蜂群算法对该模型进行优化，步骤如下：一、构造待求解的二值向量x和收益度函数f(x)；二、初始化蜂群算法参数；三、蜂群初始化；四、开始迭代，E蜂优化；五、计算E蜂招募概率；六、O蜂优化；七、更新全局最优解；八、判断迭代终止条件，若满足，优化完成；反之进行步骤九；九、判断各解是否满足S蜂启动条件，若不满足，返回步骤四；若满足，进行步骤十；十、启动S蜂，产生一个新解替换原解，返回步骤四。该方法能简化遥感成像几何模型结构，使其更准确地逼近成像几何关系。

说明书

技术领域

本发明涉及遥感影像有理函数模型结构优化方法，属于遥感成像几何模型优化技术领 域。

背景技术

遥感影像成像几何模型通常分为严格几何模型和通用几何模型两种。为了降低对用户 专业水平的需求，同时保护敏感的卫星设计参数和技术机密，包括精确的传感器几何参数 和卫星星历数据等。因而通常情况下，用户仅能获得一种通用几何模型——有理函数模型 (Rational Function Model，RFM)，以此替代以共线条件为基础的严格几何模型。然而，有 理函数模型高度依赖于模型结构，又很大程度上取决于地面控制点(Ground Control Points, GCP)数量和精度。通常RFM模型各分式的多项式次数可选择为1,2,3次，次数越高，RFM 所能描述的几何关系的精细变化程度越精确，同时需要的GCP数量也越高。一般情况下， 要各分式中各多项式的全部项均参与组成模型。经近年来研究发现，这种选择全部项的模 型构架方式，不仅增加了必要GCP数量，还可能导致模型精度下降，这是参数项过冗余 的结果。其原因是，很多参数项所描述的坐标关系并不适合当前影像所拍摄地物的真实情 况。因而，需要针对参与构成RFM的各多项式的参数项进行优化选择，实现RFM模型的 最佳结构，即目标是为了建立在保证模型精度的前提下尽可能减少参数项。然而，由于 RFM的各参数项并没有确切的物理解释，则需要实现的目标函数与决策变量之间难以架 构明确、直接的函数关系，难以用传统方法进行优化。针对单幅影像的2次RFM模型就 有超过30个参数项有待选择，则会形成成千上万种不同组合方式的解，因而一般的尝试 法进行优化是不实际的，这些问题使获得最佳结构的RFM模型成为技术难题。

多年来针对RFM参数项优化选择的方法以如遗传算法、粒子群算法等群智能算法的 处理结果相对理想，但仍存在较大的局部极值收敛和求解精度低等问题，因而存在一定的 改进空间。

发明内容

本发明为了解决现有的遥感RFM模型在描述成像几何关系时因参数项过冗余造成的 模型精度低问题，提出了一种基于蜂群算法的遥感影像有理函数模型结构优化方法。

本发明为解决上述技术问题采取的技术方案是：

步骤一：构造一个待求解的二值向量x={x_i},(i=1,2,…,D,x_i={0,1})，其长度D为待优化 的遥感影像的有理函数模型RFM(Rational Function Model，RFM)参数项系数的个数；根据 已知地面控制点(Ground Control Points,GCP)的实际坐标，将根据RFM模型解算出对应影 像点的影像坐标，与实际GCP对应影像坐标之间的误差，作为收益度函数f(x)，则显然f(x) 越小，则优化的收益度越好；进而通过以下步骤找到使f(x)达到最小值的解x；

步骤二：初始化蜂群算法参数，所述参数包括蜜蜂总数PN，将蜜蜂分为三类，S蜂 PN/2个，E蜂0个，O蜂PN/2个，设定S蜂启动阈值Limit，解空间维数D；设定迭代终 止条件：最大迭代次数(maximum iteration number,MCN)和收益度阈值(Fitness Threshold， FT)；

步骤三：蜂群初始化，令每一个蜜蜂对应一个解x，首先由PN/2个S蜂找到PN/2 个初始解x_i(i=1,2,…,PN/2)，为每个解设置一个S蜂启动计数器Failure(i)，并将其置零，计 算各解的收益度f(x_i)，存储当前全局最优解x_best，，随后S蜂全部变为E蜂，即此时PN/2 个E蜂对应PN/2个初始化解，随后进入E蜂优化阶段；

步骤四：迭代开始：首先是E蜂优化阶段：每个E蜂修正其对应解x_i(i=1,2,…,PN/2)， 计算修正后的解所对应的收益度f(x'_i)，比较f(xi)和f(x'_i)，若f(xi)≥f(x'_i)，则保留当前 E蜂所对应的解xi不变，同时Failure(i)=Failure(i)+1；若f(xi)<f(x'_i)则用x'_i替换当前E蜂 所对应的解x_i，替换完成后统一仍用x_i表示；

步骤五：计算各个E蜂招募概率$p\left(i\right)=1-\frac{f\left(x_i\right)}{\underset{i=1}{\overset{\mathrm{PN}/2}{\Sigma }}f\left(x_i\right)}(i=\mathrm{1,2},···,\mathrm{PN}/2);$

步骤六：O蜂优化阶段：对PN/2个E蜂优化后所获得的PN/2个解x_i(i=1,2,…,PN/2)， 每一个O蜂根据轮盘赌准则，首先生成一个[0,1]之间的随机数p_o，逐次与各E蜂所对应 的招募概率p(i)进行比较，若p_o<p(i)，则当前O蜂对该E蜂进行修正，每一个O蜂只修正 一个E蜂。以此原则，能保证招募概率大的E蜂，即更好的解获得更大概率的修正次数； 计算每个被O蜂修正后的解x'_i所对应的收益度f(x'_i)，与步骤四中类似，比较f(x'_i)与当 前解的收益度f(x_i)，若f(xi)≥f(x'_i)，则保留解x_i不变，同时Failure(i)=Failure(i)+1；若 f(x_i)<f(x'_i)则用x'_i替换当前解x_i，替换完成后统一仍用x_i表示；

步骤七：计算并比较经过O蜂优化后的PN/2个解的收益度f(x_i)(i=1,2,…,PN/2)，选出 对应最小f(x)值的解x_iter，进而比较f(x_iter)与当前全局最优解所对应的收益度f(x_best)，若f(x_iter) ≥f(x_best)，则保留原最优解x_best不变；若f(x_iter)<f(x_best)，则用x_iter替换x_best成为新的全局最 优解，替换完成后统一仍用x_best表示；

步骤八：判断迭代终止条件情况，若迭代次数已到达最大迭代次数MCN，或者当前 全局最优解对应收益度f(x_best)已经达到或小于收益度阈值FT，则优化完成，将所得的解对 应原RFM模型参数项，保留数值为1的项，舍弃数值为0的项，即为优化后的RFM模型； 若不满足，则进行步骤九；

步骤九：判断各解对应的S蜂启动计数器Failure(i)值是否超过阈值Limit，若未超过， 则返回步骤四；若超过，则进行步骤十；

步骤十：启动S蜂，此时设定一个额外的S蜂，产生一个新的解x，取代对应Failure(i) 超出Limit的解x_i，成为新的一个E蜂，并返回步骤四，进行下一轮迭代，直至步骤八满 足迭代终止条件。

本发明的有益效果是：

本发明针对基于有理函数模型的遥感成像模型由于参数项过冗余 (over-parameterization)而造成的模型精度低问题，利用蜂群算法针对RFM模型进行优化。 利用该方法可以有效地解决传统方法在RFM模型参数项过冗余问题，可获得使得所获有 理函数模型更准确地逼近遥感成像几何关系。

本方法用于RFM模型结构优化，它将RFM模型的各参数项对应于一个二值向量以解 决目标函数构架问题，每一解向量即对应一种参数项组合方案，进而充分利用了蜂群算法 的群智能特性，E蜂广泛搜索广域解，O对较好的结构优化方案精益求精的改进，S蜂帮 助跳出局部极值。利用该方法可以有效选择出遥感影像的RFM模型结构中关键的参数项， 去除冗余项，从而简化RFM模型结构的同时又提高了模型精度。

附图说明

图1是本发明方法的流程图；图2将全部参数项对应于一个二值向量的示意图；图3 为从地面分辨率为1m/pixel的IKONOS卫星遥感影像中选择了30GCP；图4为从地面分辨 率为0.5m/pixel的WorldView2卫星遥感影像中选择了30个GCP。

具体实施方式

具体实施方式一：结合图1说明本实施方式，本实施方式的步骤如下：

步骤一：构造一个待求解的二值向量x={x_i},(i=1,2,…,D,x_i={0,1})，其长度D为待优化 的遥感影像的有理函数模型RFM(Rational Function Model，RFM)参数项系数的个数；根据 已知地面控制点(Ground Control Points,GCP)的实际坐标，将根据RFM模型解算出对应影 像点的影像坐标，与实际GCP对应影像坐标之间误差，作为收益度函数f(x)，则显然f(x) 越小，则优化的收益度越好；进而通过以下步骤找到使f(x)达到最小值的解x；

步骤二：初始化蜂群算法参数，所述参数包括蜜蜂总数PN，将蜜蜂分为三类，S蜂 PN/2个，E蜂0个，O蜂PN/2个，设定S蜂启动阈值Limit，解空间维数D；设定迭代终 止条件：最大迭代次数(maximum iteration number,MCN)和收益度阈值(Fitness Threshold， FT)；

步骤三：蜂群初始化，令每一个蜜蜂对应一个解x，首先由PN/2个S蜂找到PN/2 个初始解x_i(i=1,2,…,PN/2)，为每个解设置一个S蜂启动计数器Failure(i)，并将其置零，计 算各解的收益度f(x_i)，存储当前全局最优解x_best，，随后S蜂全部变为E蜂，即此时PN/2 个E蜂对应PN/2个初始化解，随后进入E蜂优化阶段；

步骤四：迭代开始：首先是E蜂优化阶段：每个E蜂修正其对应解x_i(i=1,2,…,PN/2)， 计算修正后的解所对应的收益度f(x'_i)，比较f(xi)和f(x'_i)，若f(xi)≥f(x'_i)，则保留当前 E蜂所对应的解x_i不变，同时Failure(i)=Failure(i)+1；若f(x_i)<f(x'_i)则用x'_i替换当前E蜂 所对应的解x_i，替换完成后统一仍用x_i表示；

步骤五：计算各个E蜂招募概率$p\left(i\right)=1-\frac{f\left(x_i\right)}{\underset{i=1}{\overset{\mathrm{PN}/2}{\Sigma }}f\left(x_i\right)}(i=\mathrm{1,2},···,\mathrm{PN}/2);$

步骤六：O蜂优化阶段：对PN/2个E蜂优化后所获得的PN/2个解x_i(i=1,2,…,PN/2)， 每一个O蜂根据轮盘赌准则，首先生成一个[0,1]之间的随机数p_o，逐次与各E蜂所对应 的招募概率p(i)进行比较，若p_o<p(i)，则当前O蜂对该E蜂进行修正，每一个O蜂只修正 一个E蜂。以此原则，能保证招募概率大的E蜂，即更好的解获得更大概率的修正次数； 计算每个被O蜂修正后的解x'_i所对应的收益度f(x'_i)，与步骤四中类似，比较f(x'_i)与当 前解的收益度f(x_i)，若f(x_i)≥f(x'_i)，则保留解x_i不变，同时Failure(i)=Failure(i)+1；若 f(x_i)<f(x'_i)则用x'_i替换当前解x_i，替换完成后统一仍用x_i表示；

步骤七：计算并比较经过O蜂优化后的PN/2个解的收益度f(x_i)(i=1,2,…,PN/2)，选出 对应最小f(x)值的解x_iter，进而比较f(x_iter)与当前全局最优解所对应的收益度f(x_best)，若f(x_iter) ≥f(x_best)，则保留原最优解x_best不变；若f(x_iter)<f(x_best)，则用x_iter替换x_best成为新的全局最 优解，替换完成后统一仍用x_best表示；

步骤八：判断迭代终止条件情况，若迭代次数已到达最大迭代次数MCN，或者当前 全局最优解对应收益度f(x_best)已经达到或小于收益度阈值FT，则优化完成，将所得的解对 应原RFM模型参数项，保留数值为1的项，舍弃数值为0的项，即为优化后的RFM模型； 若不满足，则进行步骤九；

步骤九：判断各解对应的S蜂启动计数器Failure(i)值是否超过阈值Limit，若未超过， 则返回步骤四；若超过，则进行步骤十；

步骤十：启动S蜂，此时设定一个额外的S蜂，产生一个新的解x，取代对应Failure(i) 超出Limit的解x_i，成为新的一个E蜂，并返回步骤四，进行下一轮迭代，直至步骤八满 足迭代终止条件。

具体实施方式二：本实施方式与具体实施方式一不同点在于步骤一中的RFM模型 如公式(1):

$r=\frac{P_1(X,Y,Z)}{P_3(X,Y,Z)},c=\frac{P_2(X,Y,Z)}{P_4(X,Y,Z)}---\left(1\right)$

其中r,c为遥感影像坐标，X,Y,Z为真实地理坐标。P_i(X,Y,Z)组成形式如公式(2):

$\left(\begin{array}{cc}{P}_i(X,Y,Z)=a_{i0}+a_{i1}X+a_{i2}Y+a_{i3}Z+a_{i4}\mathrm{XY}& \\ +a_{i5}\mathrm{XZ}+a_{i6}\mathrm{YZ}+a_{i7}{X}^2+a_{i8}{Y}^2& i=\mathrm{1,2,3,4}\\ +a_{i9}{Z}^2+a_{i10}\mathrm{XYZ}+...& \end{array}\right)---\left(2\right)$

这里各分式的多项式最高次数取2，并令P₃(X,Y,Z)=P₄(X,Y,Z)，如图2，将RFM模型中 所有参数项的系数ai_j组合起来，每个系数对应一个二值数，并组成一个长度为D的向量 x^j(j=1,2…D)，若x_j取0，表示舍弃对应项；若x_j取1，表示保留对应项。

其它步骤与具体实施方式一相同。

具体实施方式三：本实施方式与具体实施方式一不同点在于步骤一中的适宜度函 数f(x)，计算公式如式(3)：

$\left(\begin{array}{c}f\left(x\right)=\sqrt{\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}({(r_n\prime -r_n)}^2+{(c_n\prime -c_n)}^2)}\\ =\sqrt{\frac{1}{N}\underset{n=1}{\overset{N}{\Sigma }}({(\frac{P_1\prime (X_n,Y_n,Z_n)}{P_3\prime (X_n,Y_n,Z_n)}-r_n)}^2+{(\frac{P_2\prime (X_n,Y_n,Z_n)}{P_4\prime (X_n,Y_n,Z_n)}-c_n)}^2)}\end{array}\right)---\left(3\right)$

将N个GCP代入优化后的RFM模型，解得的相应影像行列坐标r_n'和c_n'；式(3)中为 P₁'(X_n,Y_n,Z_n),P₂'(X_n,Y_n,Z_n),P₃'(X_n,Y_n,Z_n),P₄'(X_n,Y_n,Z_n)为优化后RFM模型的分式项，如 式(1)中形式。X_n,Y_n,Z_n(n=1,2,…,N)为N个已知GCP的实际三维坐标，r_n，c_n分别为(X_n,Y_n,Z_n) 点所对应的遥感影像行、列坐标；该收益度函数值越小，则所对应的解越好。

其它步骤与具体实施方式一或二相同。

具体实施方式四：本实施方式与具体实施方式一不同点在于步骤三中形成初始蜂 群的方法可采用随机法或经验法：

随机法即可根据公式(4)实现随机生成PN/2个初始解向量x_i^(j)(i=1,2…PN/2，j=1,2…D)， 式中为[0,1]间满足均匀分布的随机数；

经验法即可根据需求或专家经验对PN/2个初始解中的部分或全部，设定一组固定的

其它步骤与具体实施方式一、二或三相同。

具体实施方式五：实施方式与具体实施方式一不同点在于步骤四中，对于第i个E 蜂的解向量x_i的修正根据公式(5)进行：

$\left(\begin{array}{ccc}{x}_i^{\prime \left(j\right)}=[x_i^{\left(j\right)}+{\lambda }_i^{\left(j\right)}]& {\mathrm{ifx}}_i^{\prime \left(j\right)}<0& x_i^{\prime \left(j\right)}=0\\ & {\mathrm{ifx}}_i^{\prime \left(j\right)}& x_i^{\prime \left(j\right)}=1\end{array}\right)---\left(5\right)$

j为随机选择的一个欲修正决策变量；修正后的解向量为[-1,1]间满足均匀 分布的随机数。

其它步骤与具体实施方式一、二、三或四相同。

具体实施方式六：实施方式与具体实施方式一不同点在于步骤十中，S启动时，其 产生新解的方法是利用公式(4)，其中i取值为满足Failure(i)超过S蜂启动阈值Limit的对 应值。

其它步骤与具体实施方式一、二、三、四或五相同。

本发明内容不仅限于上述各实施方式的内容，如公式(4)或公式(5)，亦可采用向上取 整或向下取整的策略，针对不同的遥感影像有机会获得同样或更好的RFM优化结构。

应用效果：

将本发明方法应用于实际遥感影像的RFM模型优化，针对图3，图4中分别从IKONOS 和WorldView2卫星遥感影像中选择的30个GCP进行测试，根据公式(3)的像点平均误差 进行评价。在利用本发明所提方法优化有理函数模型前，IKONOS影像中各GCP对应像点 平均误差为1.71，优化后评价误差为0.60；WorldView2影像中各GCP对应像点平均误差 为0.87，优化后评价误差为0.62。