无人机因其低廉的成本和良好的机动性与隐身性已经成为现代战场上广泛使用的一种侦察手段。配备有效载荷完成对任务区的侦察已经成为UAV最基本也是最主要的应用之一[1]。而多UAV协同侦察能够进行优势互补与行动协调, 实现单架UAV的任务能力扩展以及多UAV系统的整体作战效能提升, 具有单架UAV不可比拟的优点, 将会是未来UAV战场应用的重要趋势[2-4]。目前, 随着战场环境复杂性的增加, 多UAV异构载荷协同执行侦察任务已成为当前UAV领域的热点问题之一。目前很多学者提出了不少相应的解决方案, 如文献[3]提出应用矩阵遗传算法解决多任务分配问题, 在此问题中考虑了任务优先权、协调时间限制以及航迹约束等条件, 并利用蒙特卡罗模拟证明了该算法的连续性和快速收敛性; 文献[5]在通信受限场景下提出了一种基于市场机制的分布式任务分配算法, 该算法具有多项式时间复杂度, 比传统启发式算法具有明显优异的性能; 文献[6]针对通信受限时异构载荷下的多UAV任务分配问题, 提出了一种优于单任务拍卖和组合拍卖的动态组合竞拍算法; 文献[7]基于多UAV动态任务分配问题, 提出两种分布式算法:基于一致性竞拍算法(consensus-based auction algorithm, CBAA)和基于一致性束算法(consensus-based bundle algorithm, CBBA), 其中CBBA算法是CBAA算法在多任务环境下的扩展算法, 2种算法都是基于市场竞拍策略的分布式任务分配算法, 同时在局部通信环境下采用一致性策略作为冲突消解机制, 该算法能够保证快速收敛, 且能够产生免冲突的任务分配方案, 是目前性能最为优异的多任务分配处理算法。
通过对相关文献的分析可以看出, 目前对多UAV协同任务决策问题的研究已经取得了不少成果。但是由于UAV面临任务环境的特殊性, 针对特定任务环境下任务决策的研究仍然较少。本文主要针对多异构型UAV协同侦察任务决策问题展开研究, 首先根据UAV的侦察载荷特性及任务区的载荷资源需求, 建立相应的“资源-需求”矩阵, 然后基于多异构型UAV的总飞行航程与到达任务区的距离影响因素指标建立相应的协同侦察任务决策模型, 最后提出了采用分布式竞拍机制的扩展CBBA(extension consensus-based bundle algorithm, ECBBA)算法进行模型求解。仿真结果表明该算法在继承CBBA算法优良性能的基础上能够快速的给出问题的解决方案, 从而为多异构型UAV协同侦察任务决策问题提供了一种快速有效的求解策略。
1 问题描述设定各异构型UAV的性能及能够携带的侦察载荷类型、战场中待侦察任务区的位置、侦察载荷资源需求已知。异构型多UAV协同侦察任务决策主要研究如何在特定时间内, 为每架UAV分配相应的侦察任务区序列, 要求分配的侦察任务区序列满足相应约束条件且使某种特定目标函数达到最大化, 如图 1所示。
2 异构型多UAV多任务区侦察决策问题建模2.1 “资源-需求”矩阵2.1.1 异构型多UAV侦察载荷资源向量设定战场上有M种共Nu架UAV, 记为UAVi(i=1, 2, 3, …, Nu), 其中UAVi拥有侦察载荷资源m种, 则UAVi的侦察载荷资源向量可表示为
(1) |
式中, RipV为UAVi携带有第G0min=0种侦察载荷资源的数量, 如:RiV= < 5, 4, 0, 7, 0 > 表示一共有5种侦察载荷资源, 其中UAVi携带侦察载荷资源Ⅰ\Ⅱ\Ⅳ的数量分别为5\4\7, 而载荷资源Ⅲ\Ⅴ未携带。一般情况下, UAVi携带某种侦察资源的数量是有限的, 应满足如下约束条件
(2) |
设定战场上有Nt个任务区, 其中任务区j的侦察资源需求向量为
(3) |
式中:RjqT为任务区j对第q种侦察资源的需求数。
2.1.3 “资源-需求”矩阵根据UAV侦察载荷资源向量、任务区载荷资源需求向量以及相应的约束条件, 构建“资源-需求”矩阵ZNu×Nt如下
(4) |
为了便于综合考虑异构型多UAV对任务区的侦察收益, 提出综合考虑UAV到达任务区的距离、任务区优先权以及UAV任务飞行总时间3个因素的侦察收益目标函数, 定义如下
(5) |
式中:Xi为UAVi所分配任务区的二进制决策向量, 记为Xi={xi1, xi2, …, xij, …, xiNt}, 若UAVi被分配执行任务区j的侦察任务, 则xij=1, 否则xij=0。Pi为UAVi所分配任务区的任务侦察执行路径, 记为Pi={pi1, pi2, …, pik, …}, 若分配给UAVi的所有任务中, 其第k次需要执行的是任务区j的侦察任务, 则pik=j。Dji为UAVi到任务区j的距离; cij(Xi, Pi, Dji)为在Xi, Pi, Dji参数给定的情况下, UAVi执行侦察任务区j所获取的侦察收益。
多异构型UAV协同侦察任务决策过程就是如何给UAVi分配所需要执行的任务区集合Xi, 在UAVi所分配的任务区集合Xi确定后如何规划其任务执行路径Pi, 从而使UAV对所有任务区的综合侦察收益最大化。
2.3 约束条件1) 每架UAV所能执行侦察的任务区数量约束
(6) |
(6) 式表示每架UAVi能够执行的侦察任务区数不能超过其能够执行的最大任务区数Lit。
2) UAV可执行侦察的总任务区数量约束
(7) |
(7) 式表示全部UAV能够执行侦察任务区的总数量约束。
3) 单个任务区的侦察任务约束
(8) |
(8) 式限定每个待侦察任务区只能由1架UAV来执行侦察。
根据上述(5)~(8) 式可将异构型多UAV协同侦察任务决策问题模型描述如下
(9) |
针对多UAV协同任务分配问题的复杂性, Luc Brunet等人在文献[6, 8]中提出了一种能够实现快速计算, 同时具有免冲突、鲁棒性强等优点的分布式CBBA算法, 该算法的效率与通信网络规模、待分配任务数量相关, 算法能够在多项式时间内给出免冲突的各UAV分配的任务集及相应的最优任务执行路径。结合多任务区侦察任务决策的特点及标准CBBA算法对影响任务决策结果等相关因素考虑的不足, 本文提出一种针对异构型多任务区侦察决策的扩展CBBA算法来进行问题求解。
CBBA算法由2个任务阶段组成:任务选择阶段和冲突消解阶段。任务选择阶段, 每个UAV局部地构造其任务束bi及相应的任务执行路径Pi。各UAV将新任务添加至当前任务束来执行所获得的收益, 与当前获胜投标列表yi中该任务的价值相比, 若其收益值大于当前获胜投标价值, 则将该任务加入到当前UAV的任务束中, 然后继续添加下一个。冲突消解阶段, 各UAV与邻接UAV通信, 利用获胜投标列表yi、获胜UAV列表zi以及表示编队成员信息新旧的时间戳向量si来解决冲突, 以在任务分配上达成一致。冲突消解引起的获胜标价和获胜UAV信息更新最终将改变各UAV的任务包和任务执行路径[9]。算法在此两阶段循环迭代, 直至达到全部任务分配完成, 如图 2所示。
3.1 信息结构假定任务场景中有Nu个异构型的UAV, Nt个待侦察任务区, 分别记为V={1, 2, …, Nu}和T={1, 2, …, Nt}, 算法中的有关信息结构定义如下[10]:
1) 任务束集(bundle)
任务束定义为bi={bi1, bi2, …, bi(lb)}, 其中元素bin∈T, lb为UAVi当前竞拍得到的任务数, bin表示UAVi竞拍得到的第n个任务的序号, 如b15=3表示UAV1竞拍得到的第5个任务是序号为3的任务。UAVi当前被分配任务束的长度为lb, UAVi能够被分配的最多任务数为Lit, 则lb≤Lit。任务束bi表示了UAVi依次竞拍获得的任务序列, 若当前UAVi尚未竞拍到任何任务, 则bi为空集, bi=ϕ。
2) 任务路径集(path)
任务路径定义为Pi={pi1, pi2, …, pi(lb)}, 其中元素pin∈T, lb为UAVi当前竞拍得到的任务数, 任务路径Pi为UAVi所获得的任务束bi满足特定目标函数最优化的任务执行顺序, 即UAVi完成所分配任务执行的顺序依次为:pi1→pi2→…→pi(lb)。
3) 时间集(time)
时间集合定义为τi={τi1, τi2, …, τi(lb)}, 其中lb为UAVi当前竞拍得到的任务数, τin∈R+表示UAVi根据任务路径Pi到达任务区pin的时间。
4) 获胜者集(winning UAVs)
获胜者集定义为zi={zi1, zi2, …, ziNt}, 其中zin∈V∪{ϕ}用来表示根据UAVi与其相邻UAV通过信息交互所获得的关于任务n的获胜者信息。如z21=4表示UAV2通过信息交互获知UAV1在竞拍中获得任务了4。若UAVi已知任务n尚未被任何UAV竞拍到, 则zi=ϕ。
5) 获胜者出价集(winning bids)
获胜者出价集为yi={yi1, yi2, …, yiNt}, 其中yin∈R+, 用来表示当前时刻各UAV对任务n竞拍时的最大出价值, 若当前尚未有UAV竞拍到任务n, 则yin=0。
6) 时间戳集(time stamps)
时间戳为si={si1, si2, …, siNu}, 其中sin∈R+, 用来表示UAVi与其邻接UAV之间的最近一次信息交互时间。
3.2 “UAV-任务”时间指标CBBA算法作为解决多Agent任务分配问题的一种有效算法, 具有运算速度快、能够避免冲突、鲁棒性强等优点已成功应用于许多实际问题的求解中。但是在对收益函数Si(j, tij)的定义中, 仅考虑了Agent i到达任务j的时间对收益(竞拍出价)的影响。实质上, 影响Agent任务分配的因素还有很多, 在多UAV多任务区侦察决策问题中, 每架UAV都更趋向于去侦察离自身距离更近的任务区, 从而可以更好地为侦察任务载荷提供更多的任务工作时间, 因此我们提出了“UAV-任务”时间指标tji
(10) |
式中:Dij为UAVi到任务区j的距离; vi为UAVi的飞行速度。
把“UAV-任务”时间tji作为影响UAV收益的一个影响因素, 则可将收益函数Si(j, tij)更新为
(11) |
式中,λ1为竞拍收益随“UAVi沿当前任务路径到达任务区j”的时间tij而变化的影响因子; λ2为“UAV-任务”时间tji对UAV竞拍收益的影响因子。
(11) 式中的λ1+λ2=1, λ1与λ2的值可通过任务决策者的经验设定或专家评估系统确定, 以反映2种指标对收益函数的影响程度。
3.3 UAV的任务续航时间指标在多UAV执行任务的过程中, UAV的续航时间是否足够其执行完所分配的任务也是必须要考虑的指标。因此, 我们提出无人机最大续航时间指标Ti, 实际飞行时间为TR
(12) |
式中,tpinpi(n+1)为UAVi从任务区pin到任务区pi(n+1)飞行的时间; stpin为UAVi执行任务区pin的持续时间; t0为UAVi从起点到其执行的第一个任务区所飞行的时间; tend为UAVi执行完最后一个任务后返回基地所需的时间; 考虑UAV最大续航时间后的约束条件为
(13) |
每架UAV根据所携带探测载荷的工作能力将对自己来说能获取收益最大的任务依次加入到自己的任务束集bi中。若将任务j加入UAVi的束集后, zij=i, yij=cij, 循环增加任务直至UAVi的束集达到最大。UAVi执行任务区j的收益由下式计算
(14) |
式中,cij(pi)为UAVi按任务路径Pi执行任务j获得的收益; SiPi为UAVi按任务路径Pi执行分配的全部任务获得的总收益; 符号|Pi|为任务路径Pi的长度; Pi⊗n{j}为将任务j添加到任务路径集Pi的第n个位置, 若为⊗end{j}则表示将任务j插入到Pi的最后一个位置。
UAVi按任务路径Pi执行完任务束集中全部任务获得的收益SiPi计算如下
(15) |
式中:SiPin (·)为UAVi沿任务路径Pi执行第n个任务获得的收益; Δ(A, B)为从位置A到位置B的时间; aloci为UAVi的当前位置; tlocj为任务区j的位置; τcur为UAVi的当前时间。
文献[11-12]中将收益函数cij(pi)定义为Si(j, tij), 由两部分组成:1) 执行任务j的固定收益; 2) Agent i沿任务路径Pi执行任务j的边际收益函数, 定义为
(16) |
式中, tij为Agent i沿任务路径Pi达到任务j的时间; λ为收益系数, 一般为常数且λ < 1;Valjt为任务j的价值系数, 0 < Valjt < 1。
任务路径Pi为UAVi所获得的任务束bi满足收益最大化的任务执行顺序, 由任务束集bi产生任务路径Pi的迭代更新公式如下
(17) |
(18) |
(13) 式和(14) 式中
(19) |
(20) |
(21) |
UAVi进行任务选择的算法实现过程如下:
Step1 按(14) 式计算UAVi对任务区j的边际收益;
Step2 由(21) 式确定UAVi是否能够竞拍到任务j;
Step3 若竞拍成功, 则由(19) 式、(20) 式求得任务j在任务路径Pi中获得最大收益的插入位置;
Step4 由(17) 式、(18) 式进行UAVi的信息更新
Step5 更新共享信息向量yi和zi;
Step6 若lb=Lt, 则退出, 否则返回Step1继续执行。
3.5 冲突消解在所有的UAV都基于自身竞拍价构建完成任务束后, UAV之间需要相互通信来交换竞拍到的任务信息, 以解决任务分配中由于任务重复分配产生的冲突问题。
在冲突消解阶段, 各个相邻UAV之间通过相互通信共享以下消息:获胜者集合zi、相应的获胜者出价集yi以及时间戳si。当所有UAV从与其邻接的UAVi接收到获胜者集合zi、相应的获胜者出价集yi信息后, 将结合自身竞拍信息来确认自己是否能够中标, 若UAVi对任务j的出价(收益)最高, 表示此刻UAVi能够竞拍到任务j, 则需要对所有任务束集中含有任务j的其他UAV的任务束集进行修改, 将他们任务束集中的任务j及其以后所有任务都释放掉以便在下一次迭代过程中进行重新分配。UAV之间完成一次通信后将更新时间戳信息, UAVi时间戳sik更新公式如下
(22) |
式中,gik为若UAVi与UAVk之间能够直接通信, 则gik=1, 否则gik=0;m为通信拓扑关系中与UAVk能够直接通信的其他UAV;
(20) 式表示若UAVi与UAVk之间是通信连通的, 则时间戳sik为最后一次的信息接收时间τr, 否则sik为与UAVi相邻的某个UAVm从通信网络中收到的UAVk信息的最新时间。
冲突消解机制主要是基于UAV之间相互信息更新的基础上, 当信息发送者UAVk将信息传递给UAVi后, 作为接收者的UAVi将依据收到的现阶段任务竞拍信息zi、yi、si作出相应的反应并触发相应的信息更新机制。UAV i对当前添加的任务j有以下3种反应机制[13]:
1) 更新机制:yij=ykj, zij=zkj;
2) 复位机制:yij=0, zij=ϕ;
3) 离开机制:yij=yij, zij=zij;
根据相应的反应处理机制, 若UAVi的zi和yi随着接收到的其他UAV的信息发生了改变, 则UAVi需要检查自身的任务束集bi中是否包含有使zi和yi集合发生改变的任务, 如果有, 则意味着发生了任务冲突, 此时UAVi应当把该任务及其之后加入到任务束集bi中的所有任务都释放掉。因此若n是UAVi任务束集bi中第一个需要删除的任务所在的位置, 即n=min{n|zi(bin)≠i}, 则对于所有满足n≥n条件的任务bin来讲, 其后的任务更新机制应为
(23) |
(24) |
(25) |
可见, 对应的任务路径集Pi、时间集合τi的长度也将会相应减小。在完成任务冲突消解阶段后, 算法将再次回到任务选择阶段继续进行任务选择, 对由于冲突消解而删除掉的任务进行重新分配, 如此循环迭代, 当所有UAV的信息集合zi和yi不再发生改变时, 表明所有UAV都分配到了相互间无冲突的任务束集, 则ECBBA算法结束并输出任务分配结果。
4 仿真实验4.1 任务场景想定假设有多架异构型UAV需要在范围为300 km×300 km的任务场景中对多个待侦察任务区进行情报侦察, 每架UAV及待侦察任务区的特征信息已知, 具体想定如下:
1) UAV特性参数设置
设定任务场景中有2种类型的UAV, 其中Ⅰ型UAV共3架, 分别记为:UAV11、UAV12、UAV13, Ⅱ型UAV共3架, 分别记为:UAV21、UAV22、UAV23, 分别从不同基地起飞执行侦察任务, 2类UAV的参数如表 1、表 2所示:
UAV标号 | 基地坐标/km | UAV标号 | 基地坐标/km |
UAV11 | (224, 289) | UAV21 | (269, 161) |
UAV12 | (265, 23) | UAV22 | (207, 44) |
UAV13 | (71, 295) | UAV23 | (152, 179) |
2) 任务区特性参数设置
假设任务场景中存在有2种类型共25个任务区域需要被侦察, 其中标号T1~T12的任务区为第Ⅰ类型, T13~T25的任务区为第Ⅱ类型, 2类任务区的重要性/价值不同, 同种任务区的重要性相同, 同时25个任务区所在位置、待侦察面积、最小侦察收益(UAV侦察该任务区最少要获得的收益, 否则侦察该任务区将失去意义)均已知, 各任务区的价值系数由情报信息预先设定, 具体参数如表 3~表 5所示。
标号 | 坐标/km |
T1 | (143, 96) |
T2 | (191, 179) |
T3 | (298, 197) |
T4 | (207, 296) |
T5 | (249, 288) |
T6 | (245, 293) |
T7 | (30, 62) |
T8 | (171, 54) |
T9 | (201, 100) |
T10 | (153, 257) |
T11 | (62, 254) |
T12 | (184, 177) |
T13 | (149, 109) |
T14 | (95, 88) |
T15 | (45, 174) |
T16 | (17, 98) |
T17 | (136, 206) |
T18 | (252, 95) |
T19 | (199, 247) |
T20 | (94, 279) |
T21 | (260, 12) |
T22 | (19, 278) |
T23 | (44, 32) |
T24 | (83, 64) |
T25 | (226, 42) |
标号 | 待侦察面积/km2 | 最小侦察收益 |
T1 | 64 | 0.25 |
T2 | 20 | 0.25 |
T3 | 33 | 0.25 |
T4 | 43 | 0.25 |
T5 | 20 | 0.25 |
T6 | 15 | 0.25 |
T7 | 59 | 0.25 |
T8 | 46 | 0.25 |
T9 | 59 | 0.25 |
T10 | 62 | 0.25 |
T11 | 79 | 0.25 |
T12 | 35 | 0.25 |
T13 | 78 | 0.25 |
T14 | 40 | 0.25 |
T15 | 39 | 0.25 |
T16 | 46 | 0.25 |
T17 | 72 | 0.25 |
T18 | 72 | 0.25 |
T19 | 59 | 0.25 |
T20 | 48 | 0.25 |
T21 | 44 | 0.25 |
T22 | 59 | 0.25 |
T23 | 46 | 0.25 |
T24 | 64 | 0.25 |
T25 | 45 | 0.25 |
针对上述仿真任务场景, 采用Matlab 2010a软件进行了仿真计算, 仿真运行环境为Intel 2.53 GHz主频, 2G内存, Windows7操作系统的商用PC机。UAV之间的通信网络为全连通, 即任意2个UAV之间都能够直接通信。在不同影响因子下进行相应的任务决策分配计算, 仿真计算结果图 3所示。
由图 3可知, 当影响因子由λ1=0.9, λ2=0.1变为λ1=0.1, λ2=0.9时, 任务区分配结果发生了相应改变, 任务区2、10、12加入到了UAV11的任务集中, 任务区7则由UAV13执行, 任务区13加入到UAV23中, 任务区14、16则被UAV22执行。产生这种结果的原因主要是:在任务区分配过程中当指挥官对第1种指标“UAV到达时间”指标的重视程度增加时, 任务分配结果偏向于使多UAV的最终总飞行航路较短, 而当对第2种指标“UAV——任务区”时间指标的重视程度增加时, 分配结果更倾向于使各UAV执行自身附近的任务区, 即在任务分配中, 指挥官对2种指标的偏向程度不同将会影响任务区分配的结果。根据UAV及任务区的位置信息可得, 所分配结果是合理且符合预期的。同时, 扩展CBBA算法在上述仿真场景下的仿真计算时间为1.411 s, 同样场景运用混合遗传算法的仿真计算时间为6.381 s, 可见扩展CBBA算法对多任务区侦察决策问题的求解效率远高于传统的遗传算法。
5 结论本文针对异构型多UAV协同侦察任务分配问题进行了建模与分析, 提出了“资源-需求”矩阵以及“UAV到达时间”指标、“UAV-任务区”时间指标, 并对异构型多UAV多任务区分配决策问题建立了相应的数学模型。结合CBBA算法及多任务区分配的特点, 提出了一种扩展的CBBA算法, 对该算法进行了相应的仿真计算与分析, 仿真结果表明在多任务区侦察决策问题中, 决策者对2种指标的偏好程度不同将会影响任务分配的结果。同时通过与经典算法分析对比, 改进的CBBA算法在保证结果稳定可靠的同时, 运行时间要明显优于经典的启发式优化算法。
在真实的作战任务场景中常常会面临诸如地面突发威胁, 待侦察任务区数量的突然增加或减少、任务区的移动以及UAV因故障而失效等各种突发情况, 这些变化都是难以预料的。因此, 研究在线异构型多任务区侦察的最优化决策问题将是后续工作的研究重点。
[1] |
黄丁才. UAV侦察机航线与传感器规划方法研究[D]. 长沙: 国防科技大学, 2009 Huang Dingcai. Research on Path and Sensor Planning of Unmanned Reconnaissance Aerial Vehicle[D]. Changsha, National University of Defense Technology, 2009(in Chinese) |
[2] |
田菁. 多UAV协同侦察任务规划问题建模与优化技术研究[D]. 长沙: 国防科技大学, 2007 Tian Jing. Modeling and Optimization Methods for Multi-UAV Cooperative Reconnaissance Mission Planning Problem[D]. Changsha, National University of Defense Technology, 2007(in Chinese) |
[3] | Shima T, Rasmussen S J, Sparks A G, Passino K M. Multiple Task Assignments for Cooperating Uninhabited Aerial Vehicles Using Genetic Algorithms[J]. Elsevier Computers & Operations Research, 2006, 33(11): 3252-3269. |
[4] | Tim Bakker, Robert H Klenke. Dynamic Multi-Task Allocation for Collaborative Unmanned Aircraft Systems[J]. AIAA Sci Tech, 2014, 24(11): 30-54. |
[5] | Gyeongtaek Oh, Youdan Kim. Market-Based Task Assignment for Cooperative Timing Missions over Networks with Limited Connectivity[J]. AIAA Sci Tech, 2015, 1(15): 1-19. |
[6] | Brunet L, Choi H L, How J. Consensus-Based Auction Approaches for Decentralized Task Assignment[C]//AIAA Guidance, Navigation and Control Conference and Exhibit, 2008:1-24 |
[7] | Simon Hunt, Qinggang Meng, Chris Hinde. A Consensus-Based Grouping Algorithm for Multi-agent Cooperative Task Allocation with Complex Requirements[J]. Springer, 2014, 6(3): 338-350. |
[8] | Choi Hanlim, Luc Brunet, Jonathan P How. Consensus-Based Decentralized Auctions for Robust Task Allocation[J]. IEEE Trans on Robotics, 2010, 25(4): 1552-3098. |
[9] | Andrew Whitten. Decentralized Planning for Autonomous Agents Cooperating in Complex Missions[D]. Massachusetts, Massachusetts Institute of Technology, 2010:3-20 |
[10] | Tim Bakker, Robert H Klenke. Dynamic Multi-Task Allocation for Collaborative Unmanned Aircraft Systems[J]. AIAA Sci Tech, 2014, 1(13): 1-19. |
[11] | Matthew Argyle, David W.Casber, Randy Beard. A Multi-Team Extension of the Consensus-Based Bundle Algorithm[C]//American Control Conference, San Francisco, 2011:5376-5381 |
[12] | Travis Mercker. An Extension of Consensus-Based Auction Algorithm for Decentralized, Time-Constrainted Task Assignment[C]//American Control Conference, Baltimore, 2010:6324-6329 |
[13] | Matthew E. Argyle. Multi-Team Consensus Bundle Algorithm[J]. Springer Netherlands, 2015, 24(12): 1491-1507. |