﻿ 过载约束下的大机动目标协同拦截
 文章快速检索 高级检索

Cooperative interception against highly maneuvering target with acceleration constraints
XIAO Wei, YU Jianglong, DONG Xiwang, LI Qingdong, REN Zhang
Science and Technology on Aircraft Control Laboratory, Beihang University, Beijing 100083, China
Abstract: In this paper, the cooperative interception of multiple inferior missiles with acceleration constraints to intercept a highly maneuvering target is studied using the nonlinear interception geometry. Firstly, based on the establishment of the reachable region of the missile, the feasible region of the missile, and the escape region of the target, a cooperative interception strategy based on the coverage of escape region using nonlinear interception geometry and a design method based on the standard trajectory are proposed. Secondly, after giving the form of co-operative guidance law, the relationship between the initial position of the terminal guidance phase, the parameters of the guidance law, and the coverage area of the maneuvering range of the target is studied. A numerical algorithm is designed to realize the allocation of the coverage area, the design of the cooperative guidance law, and the configuration of the initial interception position. Finally, the theoretical results are simulated. The results show that multiple missiles with weak maneuverability can intercept highly maneuvering target through a reasonable configuration of the initial interception position and a reasonable design of parameters in the cooperative guidance law.
Keywords: guidance    cooperative guidance    cooperative interception    acceleration constraints    highly maneuvering target

1 比例导引局限性分析

 图 1 平面拦截几何 Fig. 1 Geometry of planar engagement
 $\left\{ {\begin{array}{*{20}{l}} {\dot r = {V_T}{\rm{cos}}{\eta _T} - {V_M}{\rm{cos}}{\eta _M}}\\ {\dot r = {V_M}{\rm{sin}}{\eta _M} - {V_T}{\rm{sin}}{\eta _T}}\\ {{\theta _M} = {\eta _M} + {\sigma _M}}\\ {{\theta _T} = {\eta _T} + {\sigma _T}}\\ {{{\dot \sigma }_M} = {a_M}/{V_M}}\\ {{{\dot \sigma }_T} = {a_T}/{V_T}} \end{array}} \right.$ （1）

 ${a_M} = N|\dot r|\dot \theta$ （2）

 $\dot r\ddot \theta = - (N|\dot r|{\rm{cos}}{\eta _M} + 2\dot r)(\dot \theta - {\dot \theta ^*})$ （3）

 $|{a_{T,f}}| \le \frac{{N{\rm{cos}}{\eta _{M,f}} - 2}}{{N|{\rm{cos}}{\eta _{T,f}}|}}{a_{M, {\rm{Max}} }}$ （4）

2 基于覆盖的协同拦截策略与方法 2.1 基于覆盖的协同拦截策略

 图 2 基于逃逸域覆盖的协同拦截策略 Fig. 2 Cooperative interception strategy based on coverage of escape region

1) ∀aT, i, sUT:=[-aT, Max, aT, Max]，确定初始位置Pi, s=[ηM, i, 0, s, ηT, i, 0, s, ri, 0, s]T和制导律gM, i，以满足标准弹道的要求:拦截过程中过载aM, i=gM, i的评价函数J=J(aM, i)最小。

2) 确定导弹i的覆盖区域UM, i(aT, i, s)。UM, i(aT, i, s)是$\hat U$M, i(aT, i, s)的闭包。其中$\hat U$M, i(aT, i, s)由3个条件定义: ① $\hat U$M, i(aT, i, s)是aT, i, s的邻域；②∀${\tilde a}$T, i, s$\hat U$M, i(aT, i, s)机动的目标，拦截过程中导弹的过载不饱和，即|aM, i(t)| < aM, Max(∀0≤t < tf); ③ ∀ε>0，$\exists \delta$∈(0, ε)，当目标以sup$\hat U$M, i+δ或inf$\hat U$M, i-δ机动时，$\exists {t_1}$∈[0, tf)，使得|aM, i(t1)|≥aM, Max

3) 分配多弹的覆盖区域UM, i(aT, i, s)以均匀地覆盖目标机动范围UT，即UT⊂∪UM, i

2.2 协同拦截制导律

 ${a_{M,i}} = N(|{\dot r_i}|{\dot \theta _i} - {B_i})$ （5）

aT, i为目标的实际逃逸加速度，将式(5)代入式(1)，整理可得

 ${r_i}{\ddot \theta _i} = - (N|{\dot r_i}|{\rm{cos}}{\eta _{M,i}} + 2{\dot r_i})({\dot \theta _i} - \dot \theta _i^*)$ （6）

 $\dot \theta _i^* = \frac{1}{{N|{{\dot r}_i}|{\rm{cos}}{\eta _{M,i}} + 2{{\dot r}_i}}}({a_{T,i}}{\rm{cos}}{\eta _{T,i}} + N{B_i}{\rm{cos}}{\eta _{M,i}})$ （7）

 ${\dot \theta _{i,f}}: = {\dot \theta _i}{|_{t = {t_f}}} = \dot \theta _i^*{|_{t = {t_f}}}$ （8）

 ${a_{M,i,0}} = {a_{M,i,f}} = 0$ （9）

 ${B_i} = - 0.5{a_{T,i,s}}{\rm{cos}}{\eta _{T,i,f,s}}$ （10）

 $\begin{array}{*{20}{c}} {({V_T}{\rm{cos}}{\eta _{T,i,0,s}} - {V_M}{\rm{cos}}{\eta _{M,i,0,s}}) \cdot ({V_M}{\rm{sin}}{\eta _{M,i,0,s}} - }\\ {{V_T}{\rm{sin}}{\eta _{T,i,0,s}}) + {B_i}{r_{i,0,s}} = 0} \end{array}$ （11）

 ${V_M}{\rm{sin}}{\eta _{M,i,f}} - {V_T}{\rm{sin}}{\eta _{T,i,f}} = 0$ （12）

 ${t_{i,f}} = \frac{{{r_{i,0}}}}{{|{{\dot r}_{i,0}}|}} = \frac{{{r_{i,0}}}}{{{V_M}{\rm{cos}}{\eta _{M,i,0}} - {V_T}{\rm{cos}}{\eta _{T,i,0}}}}$ （13）

 ${\dot \theta _i} = {\dot \eta _{T,i}} + \frac{{{a_{T,i}}}}{{{V_T}}}$ （14）

 ${\dot \theta _i} = {\dot \eta _{M,i}} + {a_{M,i}}/{V_M}$ （15）

 ${\dot \eta _{M,i}} + \left( {\frac{{N|{{\dot r}_i}|}}{{{V_M}}} - 1} \right){\dot \theta _i} - \frac{{N{B_i}}}{{{V_M}}} = 0$ （16）

 $\begin{array}{l} {\eta _{M,i}} + (N|{{\dot r}_i}|/{V_M} - 1)({{\dot \eta }_{T,i}} + {a_T}/{V_T}) - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} N{B_i}/{V_M} = 0 \end{array}$ （17）

 $\begin{array}{l} {\eta _{M,i,f}} - {\eta _{M,i,0}} - N{B_i}{t_{i,f}}/{V_M} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {K_{i,0}}({\eta _{T,i,f}} - {\eta _{T,i,0}} + {a_{T,i}}{t_{i,f}}/{V_T}) = 0 \end{array}$ （18）

 $\left\{ {\begin{array}{*{20}{l}} {{V_M}{\rm{sin}}{\eta _{M,i,f}} - {V_T}{\rm{sin}}{\eta _{T,i,f}} = 0}\\ {{\eta _{M,i,f}} + {K_{i,0}}{\eta _{T,i,f}} + z = 0} \end{array}} \right.$ （19）

 $\begin{array}{*{20}{l}} {z = - {\eta _{M,i,0}} + {K_{i,0}}( - {\eta _{T,i,0}} + {a_{T,i}}{t_{i,f}}/{V_T}) - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} N{B_i}{t_{i,f}}/{V_M}} \end{array}$

 $\left\{ \begin{array}{l} ({V_T}{\rm{cos}}{\eta _{T,i,0,s}} - {V_M}{\rm{cos}}{\eta _{M,i,0,s}}) \cdot ({V_M}{\rm{sin}}{\eta _{M,i,0,s}} - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {V_T}{\rm{sin}}{\eta _{T,i,0,s}}) - 0.5{a_{T,i,s}}{r_{i,0,s}}{\rm{cos}}{\eta _{T,i,f,s}} = 0\\ {V_M}{\rm{sin}}{\eta _{M,i,f,s}} - {V_T}{\rm{sin}}{\eta _{T,i,f,s}} = 0\\ {\eta _{M,i,f,s}} + {K_{i,0,s}}{\eta _{T,i,f,s}} + 0.5{a_{T,i,s}}N{t_{i,f,s}}{\rm{cos}}{\eta _{T,i,f,s}}/\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {V_M} + \bar z = 0 \end{array} \right.$ （20）

 $\left| {\frac{{N({{\tilde a}_{T,i}}{\rm{cos}}{{\tilde \eta }_{T,i,f}} - {a_{T,i,s}}{\rm{cos}}{\eta _{T,i,f,s}})}}{{N{\rm{cos}}{{\tilde \eta }_{M,i,f}} - 2}}} \right| \le {a_{M, {\rm{Max}} }}$ （21）

UM, i(aT, i, s)是aT, i, s的邻域，且满足: ① ∀${\tilde a}$T, iUM, i(aT, i, s)，${\tilde a}$T, i满足式(21)；② ∀ε>0，$\exists$δ∈(0, ε)，机动值${\tilde a}$T, i取为UmaxM, i(aT, i, s)+δ或minUM, i(aT, i, s)-δ时，${\tilde a}$T, i皆不满足式(21), 则UM, i(aT, i, s)是命中点处导弹的过载约束对覆盖区域UM, i(aT, i, s)的一个限定。因此有UM, i(aT, i, s)⊂UM, i(aT, i, s)。在线性度较好的情况下，由于初始时刻过载为零，故可近似认为导弹的过载绝对值在拦截过程中递增至命中点的过载值，在拦截过程中没有饱和，则该命中点过载约束下的邻域UM, i(aT, i, s)即覆盖区域UM, i(aT, i, s)本身。

2.3 多弹覆盖区域的分配算法

 图 3 2种理想覆盖模式 Fig. 3 Two ideal coverage modes

Table
 算法1  理想覆盖模式1的实现算法 Algorithm 1  Algorithm for implement ideal coverage Mode 1 1 输入:导弹、目标参数, 初始阵位约束 2 求解[B1, P1, sT, uT, 1, s, UM, 1]T使得u1, low=-1 3 令i=2 4 while ui-1, up < 1 5 求解[Bi, Pi, sT, uT, i, s, UM, i]T使得ui, low=ui-1, up 6 if ui, up>1 7 求解(Bi, Pi,sT, uT, i, s, UM, i)T使得ui, up=1 8 i=i+1 9 令Nm=i-1 10 令i=2, Δ=(uNm-1, up-uNm, low)/(Nm-1) 11 while i≤Nm-1 12 uT, i, s=uT, i, s-(i-1)Δ 13 以新的uT, i, s为期望值, 求解[Bi, Pi, sT, UM, i]T 14 i=i+1 15 输出:{Bi, Pi, s, uT, i, s, UM, i|i=1, 2, …, Nm}和Nm
2.4 制导律参数动态调整策略

 ${\rm{ZE}}{{\rm{M}}_i} = {r_i}\sqrt {\frac{{{{({r_i}{{\dot \theta }_i})}^2}}}{{{{({{\dot r}_i})}^2} + {{({r_i}{{\dot \theta }_i})}^2}}}}$ （22）

 ${\rm{ZE}}{{\rm{M}}_{i,s}} = \left\{ {\begin{array}{*{20}{c}} {{r_{i,s}}\sqrt {\frac{{{{({r_{i,s}}{{\dot \theta }_{i,s}})}^2}}}{{{{({{\dot r}_{i,s}})}^2} + {{({r_{i,s}}{{\dot \theta }_{i,s}})}^2}}}} }&{{\rm{ if}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {r_{i,s}} > \kappa }\\ 0&{{\rm{ otherwise }}} \end{array}} \right.$ （23）

 $\left\{ \begin{array}{l} {{\dot r}_{i,s}} = {V_T}{\rm{cos}}(({\theta _{i,s}} - {\theta _{i,s,0}}) - (t \cdot {a_{T,i,s}}/{V_T} - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\eta _{T,i,0,s}})) - {V_M}{\rm{cos}}({\theta _{i,s}} - {\sigma _{M,i,s}})\\ {r_{i,s}}{{\dot \theta }_{i,s}} = {V_M}{\rm{sin}}({\theta _{i,s}} - {\sigma _{M,i,s}}) - {V_T}{\rm{sin}}(({\theta _{i,s}} - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\theta _{i,s,0}}) - (t{a_{T,i,s}}/{V_T} - {\eta _{T,i,0,s}}))\\ {{\dot \sigma }_{M,i,s}} = N(|{{\dot r}_{i,s}}|{{\dot \theta }_{i,s}} - {B_i})/{V_M} \end{array} \right.$ （24）

ZEMi, s可在拦截过程中实时计算得到，也可以在末制导开始前计算完成并存储起来而在拦截过程中通过插值获得。

 ${B_i}({t_{k + 1}}) = \left\{ {\begin{array}{*{20}{l}} {{c_i}{B_{i,0}} + (1 - {c_i}){B_{{i_1}}}({t_k})}&{{\rm{ if}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \varLambda > \rho }\\ {{B_i}({t_k})}&{{\rm{ otherwise }}} \end{array}} \right.$ （25）

3 仿真验证

 $\left\{ {\begin{array}{*{20}{l}} {{\eta _{M,i,0,s}} = 0}\\ {{r_{i,0,s}} = {r_0}} \end{array}} \right.$ （26）

 $\left\{ {\begin{array}{*{20}{l}} {{\eta _{T,i,0,s}} = - \pi }\\ {{r_{i,0,s}} = {r_0}} \end{array}} \right.$ （27）

 导弹 Bi ηT, i, 0, s/(°) UM, i uT, i, s M1 -11.44 -183.6 [-1.00, -0.17] -0.59 M2 0 -180.0 [-0.41, 0.41] 0 M3 11.44 -176.4 [0.17, 1] 0.59

 图 4 u=±1, 0时的拦截弹道 Fig. 4 Interception trajectory when u=±1, 0
 图 5 u=±1, 0时的过载 Fig. 5 Acceleration when u=±1, 0

 图 6 脱靶量分布情况 Fig. 6 Distribution of miss distance

 导弹 Bi ηM, i, 0, s/(°) UM, i uT, i, s M1 -11.36 -3.10 [-1.00, -0.17] -0.58 M2 0 0 [-0.41, 0.41] 0 M3 11.36 3.10 [0.17, 1] 0.58

 图 7 u=±1, 0时的拦截弹道 Fig. 7 Interception trajectory when u=±1, 0
 图 8 u=±1, 0时的过载 Fig. 8 Acceleration when u=±1, 0

 图 9 多弹齐射式作战拦截bang-bang机动目标的弹道 Fig. 9 Trajectory of multiple missiles against bang-bang maneuvering target
 图 10 子母弹分离式作战拦截bang-bang机动目标的弹道 Fig. 10 Trajectory of separated missiles against bang-bang maneuvering target
 图 11 多弹齐射式作战拦截bang-bang机动目标的ZEM Fig. 11 ZEM of multiple missiles against bang-bang maneuvering target
 图 12 子母弹分离式作战拦截bang-bang机动目标的ZEM Fig. 12 ZEM of separated missiles against bang-bang maneuvering target
4 结论

1) 研究了多枚过载受限的弱机动能力导弹拦截强机动能力目标的协同拦截问题，针对平面非线性拦截模型，在建立可达域、可行域和逃逸域这3个概念的基础上提出了基于覆盖的协同拦截策略，并提出了基于标准弹道的设计方法。

2) 研究了协同制导律的形式，并给出了拦截导弹的末制导初始阵位、制导律参数以及导弹对目标机动的覆盖区域的关系。

3) 给出了多弹对目标的2种理想覆盖模式，并分别设计了实现这2种覆盖模式的数值求解算法。

4) 针对多弹齐射作战和子母弹分离作战两种作战模式，对所提的方法进行了仿真。仿真结果验证了本文研究的基于覆盖的协同拦截策略、基于标准弹道的设计方法以及协同制导律参数和中末制导交班阵位的数值求解算法的有效性。

 [1] JEON I S, LEE J I, TAHK M J. Impact-time-control guidance law for anti-ship missiles[J]. IEEE Transactions on Control Systems Technology, 2006, 14(2): 206-266. Click to display the text [2] HE S M, LIN D F. Three-dimensional optimal impact time guidance for antiship missiles[J]. Journal of Guidance, Control, and Dynamics, 2019, 42(4): 941-948. Click to display the text [3] ZHAO Q L, DONG X W, LIANG Z X, et al. Distributed cooperative guidance for multiple missiles with fixed and switching communication topologies[J]. Chinese Journal of Aeronautics, 2017, 30(4): 1570-1581. Click to display the text [4] KUMAR S R, GHOSE D. Impact time guidance for large heading errors using sliding mode control[J]. IEEE Transactions on Aerospace and Electronic Systems, 2015, 51(4): 3123-3138. Click to display the text [5] LEE J K, JEON I S, TAHK M J. Guidance law to control impact time and angle[J]. IEEE Transactions on Aerospace and Electronic Systems, 2007, 43(1): 301-310. Click to display the text [6] CHEN Y D, WANG J N, WANG C Y, et al. Impact time and angle control based on constrained optimal solutions[J]. Journal of Guidance, Control, and Dynamics, 2016, 39(10): 2445-2451. Click to display the text [7] KANG S, KIM H J. Differential game missile guidance with impact angle and time constraints[C]//Proceedings of the International Federation of Automatic Control World Congress, 2011: 3920-3925. [8] JUNG B, KIM Y. Guidance laws for anti-ship missiles using impact angle and impact time[C]//AIAA Guidance, Navigation, and Control Conference and Exhibit. Reston: AIAA, 2006. [9] HARL N, BALAKRISHNAN S N. Impact time and angle guidance with sliding mode control[J]. IEEE Transactions on Control Systems Technology, 2012, 20(6): 1436-1449. Click to display the text [10] SHAFERMAN V, SHIMA T. Cooperative differential games guidance laws for imposing a relative intercept angle[J]. Journal of Guidance, Control, and Dynamics, 2017, 40(10): 2465-2480. Click to display the text [11] 赵世钰, 周锐. 基于协调变量的多导弹协同制导[J]. 航空学报, 2008, 29(6): 1605-1611. ZHAO S Y, ZHOU R. Multi-missile cooperative guidance using coordination variables[J]. Acta Aeronautica et Astronautica Sinica, 2008, 29(6): 1605-1611. (in Chinese) Cited By in Cnki (110) | Click to display the text [12] JEON I S, LEE J I, TAHK M J. Homing guidance law for cooperative attack of multiple missiles[J]. Journal of Guidance, Control, and Dynamics, 2010, 33(1): 275-280. Click to display the text [13] SHAFERMAN V, OSHMAN Y. Cooperative interception in a multi-missile engagement[C]//AIAA Guidance, Navigation, and Control Conference and Exhibit. Reston: AIAA, 2009. [14] YUKSEK B, URE N K, INALHAN G. Cooperative interception of a highly manoeuvrable aerial target[C]//AIAA Guidance, Navigation, and Control Conference. Reston: AIAA, 2018. [15] SU W S, SHIN H S, CHEN L, et al. Cooperative interception strategy for multiple inferior missiles against one highly maneuvering target[J]. Aerospace Science and Technology, 2018, 80: 91-100. Click to display the text [16] SU W S, LI K B, CHEN L. Coverage-based cooperative guidance strategy against highly maneuvering target[J]. Aerospace Science and Technology, 2017(71): 147-155. Click to display the text [17] SU W S, LI K B, CHEN L. Coverage-based three-dimensional cooperative guidance strategy against highly maneuvering target[J]. Aerospace Science and Technology, 2019, 85: 556-566. Click to display the text [18] 钱杏芳, 林瑞雄, 赵亚男. 导弹飞行力学[M]. 北京: 北京理工大学出版社, 2008: 90-135. QIAN X F, LIN R X, ZHAO Y N. Missile flight me-chanics[M]. Beijing: Beijing Institute of Technology Press, 2008: 90-135. (in Chinese) [19] NESLINE F W, ZARCHAN P. A new look at classical vs modern homing missile guidance[J]. Journal of Guidance, Control, and Dynamics, 1981, 4(1): 78-85. Click to display the text [20] ZARCHAN P. Tactical and strategic missile guidance[M]. 6th ed. Reston: AIAA, 2012. [21] YUAN P J, CHERN J S. Analytic study of biased proportional navigation[J]. Journal of Guidance, Control, and Dynamics, 1992, 15(1): 185-190. Click to display the text
http://dx.doi.org/10.7527/S1000-6893.2019.23777

0

#### 文章信息

XIAO Wei, YU Jianglong, DONG Xiwang, LI Qingdong, REN Zhang

Cooperative interception against highly maneuvering target with acceleration constraints

Acta Aeronautica et Astronautica Sinica, 2020, 41(S1): 723777.
http://dx.doi.org/10.7527/S1000-6893.2019.23777