2. 中国科学院太阳活动重点实验室, 北京 100101
2. Key Laboratory of Solar Activity, Chinese Academy of Sciences, Beijing 100101, China
太阳活动区是太阳强磁场区域，是太阳活动的主要能量来源。剧烈的太阳活动导致恶劣的空间天气，比如耀斑和日冕物质抛射对地球上电磁通信、电力系统、无线电传输等产生不良影响。因此，对太阳活动区进行检测与跟踪，分析其演化规律，对人类的空间探索和地球生活有重要的意义[1-2]。
目前，用来解决太阳活动区检测与跟踪的方法主要采用传统的图像处理技术[1-11]，检测主要采用强度阈值法或区域生长法等，这些方法一般需要设置参数，如强度阈值、开闭算子阈值和区域生长的边界阈值。这些参数对传统的方法来说至关重要，文[2, 6]对参数进行了详细讨论，并给出了设置参数的一些准则。跟踪太阳活动区主要根据纬向较差自转定律预测其位置，然后根据欧氏距离或特征的相关性进行目标关联[3-4]。
这些方法虽然较好地实现了太阳活动区的检测与跟踪，但是仍存在一些问题。比如，文[2-4]的方法容易导致一个太阳活动区被误检测为多个，或者多个相邻的太阳活动区被误检测为一个，从而导致误跟踪；文[11]的美国国家海洋和大气管理局(National Oceanic and Atmospheric Administration, NOAA)编号被广泛应用，但存在未能及时标注新浮现的太阳活动区和仍标注了已经在全日面像上消失的太阳活动区等现象。
近年来，随着深度学习[12]的普及和神经网络算法的发展，提出了一些基于深度学习的目标检测和跟踪算法。典型的目标检测算法有R-CNN[13], Faster R-CNN[14], YOLO[15], YOLOv3[16]等。其中，YOLOv3借鉴残差网络的思想采用了新的特征提取网络Darknet-53，并结合多尺度特征进行目标检测，取得了较好的检测性能。典型的多目标跟踪方法有EAMTT[17], STAM[18], DeepSort[19]等。其中，DeepSort具有较好的跟踪性能和较高的实时性能，越来越多的深度学习方法应用于天文领域[20-22]。
本文针对太阳活动区的演化特性，采用YOLOv3-spp作为检测算法，DeepSort作为跟踪算法并进行改进，提出了一种太阳活动区检测和跟踪方法。该算法可以有效检测并跟踪不同时间间隔序列图像中的太阳活动区。
1 数据和数据集本文采用太阳动力学天文台(Solar Dynamics Observatory, SDO)卫星的日震及磁场成像仪(Helioseismic and Magnetic Imager, HMI)与日球天文台(Solar and Heliospheric Observatory, SOHO)的迈克逊多普勒成像仪(Michelson Doppler Imager, MDI)两个设备的太阳全日面纵向磁图作为数据来源，制作了太阳活动区检测与跟踪的深度学习训练和验证数据集。
表 1列出了分别用于太阳活动区检测与跟踪的数据集。D1共选取3 017张全日面图像用于制作太阳活动区检测网络的训练和验证集，共选取16 361个太阳活动区样本，其中80%作为训练样本，用于训练深度学习网络模型，20%作为验证样本，在训练过程中用于防止过拟合的验证。D2共选取33 515个太阳活动区样本用于跟踪过程中重识别网络的训练和验证，将太阳活动区局部区域缩放为128 × 128像素尺寸，同样80%用于训练，20%用于验证。本文根据SolarMonitor网站上标注的太阳活动区编号，结合人工辅助判断的方式，采用labelImg软件标注数据集。
| Purpose | Dataset | Data source | Stage | Time/Interval | Sample number |
| Training/ Validation set |
D1 | HMI | Detection | 2011.01.01~2017.12.31 | 16 361 |
| D2 | HMI | Tracking | 2011.11.01~2011.11.30/12 min 2017.03.20~2017.04.20/12 min |
33 515 |
活动区检测和跟踪算法主要分为两部分：检测与跟踪。在检测部分采用YOLOv3-spp模型，该模型的网络结构如图 1，由Darknet-53(图 1(a))与yolo层(图 1(b))两部分组成，分别用于提取图像特征与多尺度预测。主要方法步骤如下：(1)将全日面图像压缩为608 × 608的输入网络；(2)输入图像经过一个卷积单元(Convolutional)得到32个608 × 608的特征图，卷积单元包含卷积、批量归一化和leaky relu激活函数3部分；(3)特征图经过步长为2、64个卷积核的卷积单元实现降采样，得到64个304 × 304的特征图，然后送入残差单元(Residual)，输出64个304 × 304的特征图(1x残差块)，1x表示此残差块由1个残差单元组成，残差单元是将输入与两个卷积单元进行残差操作；(4)上一步得到的特征图经过降采样和2x残差块得到128个152 × 152的特征图，其中，2x表示此残差块由2个残差单元组成；(5)依次执行降采样、8x残差块、降采样、8x残差块、降采样，4x残差块操作，3个残差块依次得到256个76 × 76、512个38 × 38和1 024个19 × 19的特征图；(6)将8x，8x和4x 3个残差块的结果送入yolo层(图 1(b))；(7)19 × 19的特征图经过金字塔卷积集(Convolutional-spp Set)等操作得到第1次的检测结果(Output1)，同时，将金字塔卷积集的结果上采样至38 × 38与第4个残差块38 × 38的特征图融合，经过卷积集(Convolutional Set)等操作得到第2次的检测结果(Output2)，依次类推，得到第3次的检测结果(Output3)。其中，卷积集操作包含5个卷积单元，金字塔卷积集是在卷积集基础上添加了空间金字塔池化[23]；(8)将以上3次检测结果使用非极大抑制算法去除交并比较大的预测框，保留置信得分较高的预测框作为太阳活动区的检测结果。该网络的损失值由3部分组成：目标定位偏移量损失、目标置信度损失和目标类别损失，其中，目标置信度损失和类别损失采用二元交叉熵损失，目标定位损失则采用均方误差进行计算。
|
| 图 1 Yolov3-spp网络结构示意图 Fig. 1 YOLOv3-spp network structure diagram |
在跟踪部分，基于DeepSort算法对其进行改进，主要步骤如下：
(1) 根据太阳活动区的特性改进DeepSort的特征提取网络，将其称之为太阳活动区重识别网络，如图 2。首先将检测到的太阳活动区缩放为128 × 128，经过2次3 × 3的卷积、1次最大池化以及6个残差块，特征图减小到16 × 16，再经过一个全连接层得到128维的全局特征向量，最后采用L2归一化得到特征向量。在训练过程中，采用Cosine Softmax函数得到分类结果，采用交叉熵计算损失值。需要说明的是，检测部分的网络和太阳活动区重识别网络虽然都提取太阳活动区的特征，但是作用不同。检测部分的网络主要用于检测太阳活动区，而太阳活动区重识别网络主要为了判定太阳活动区在演化过程中的特征相似性。
|
| 图 2 太阳活动区重识别网络结构示意图 Fig. 2 Active regions re-identification network structure diagram |
(2) 用太阳活动区重识别网络提取已跟踪到的太阳活动区(称为轨迹，假设数量为L)和当前帧检测到的太阳活动区(假设数量为J)的特征向量，计算之间的余弦距离，得到关联矩阵CL×J。然后，根据前一帧所有太阳活动区的位置，采用纬向较差自转[24-25]预估在当前帧的位置。接着，计算所有预测位置和所有检测位置的欧氏距离，若过大则将关联矩阵中对应的值设置为无穷大。
(3) 对关联矩阵CL×J采用匈牙利匹配算法[26]，得到3类匹配结果：匹配到检测结果的轨迹、未匹配到轨迹的检测结果和未匹配到检测结果的轨迹。对于匹配到检测结果的轨迹则继续跟踪；未匹配到轨迹的检测结果则为之分配一个新的轨迹，当新轨迹连续3帧匹配到检测结果时，则认为该太阳活动区产生，否则，删除该轨迹；对于未匹配到检测结果的轨迹则提取使用纬向较差自转预测结果的特征向量，与前一帧对应的特征向量计算余弦距离，若其值在一定范围内(本文实验结果为0.7)则采纳该预测结果作为轨迹的一部分，否则该活动区结束，但如果连续3帧均为预测结果，也判定该活动区结束。
3 实验结果与分析本文采用日震及磁场成像仪和迈克逊多普勒成像仪的4组全日面序列图像进行测试。这4组数据分别代表不同的观测设备和活动区数量。表 2列出了数据的基本信息，为了便于描述，分别用D3~D6表示。其中，D3和D4选取迈克逊多普勒成像仪纵向磁图，D5和D6选取日震及磁场成像仪纵向磁图。D4和D6选取太阳活动区数量较多的序列图像，而D3和D5选取太阳活动区数量较少的序列图像。由于时间间隔较长的数据太阳活动区的特征变化大，跟踪难度较高，因此4组数据均采用时间间隔较长的序列图像(为了便于与SolarMonitor结果进行比较，数据与SolarMonitor提供的数据时间基本一致)。
| Dataset | Data source | Time/Interval | NOAA | ARDTM | ARDTM∩NOAA | NOAA-(ARDTM∩NOAA) | ARDTM-(ARDTM∩NOAA) | (ARDTM∩NOAA)/NOAA |
| D3 | MDI | 2005.03.01~2005.03.31/24 h | 68 | 72 | 62 | 6 | 10 | 91.18% |
| D4 | MDI | 2000.07.12~2000.08.11/24 h | 410 | 415 | 398 | 12 | 17 | 97.07% |
| D5 | HMI | 2010.11.05~2010.12.04/24 h | 103 | 108 | 91 | 12 | 17 | 88.35% |
| D6 | HMI | 2014.10.01~2014.10.31/24 h | 225 | 277 | 204 | 21 | 73 | 90.67% |
| Total | 806 | 872 | 755 | 51 | 117 | 91.82% |
表 2列出了4组数据集的跟踪结果，并与NOAA进行了对比。包括SolarMonitor上标注的NOAA活动区数量、ARDTM方法标注的太阳活动区数量、ARDTM结果与NOAA结果相同的活动区的数量(ARDTM∩NOAA)、NOAA活动区数量与交集的差值(NOAA-(ARDTM∩NOAA))、ARDTM活动区数量与交集的差值(ARDTM-(ARDTM∩NOAA))、交集与NOAA活动区数量的百分比((ARDTM∩NOAA)/NOAA)。NOAA一共标注了806个活动区，ARDTM标注了872个，两者的交集为755个。ARDTM标注的活动区比NOAA的多一些，主要由以下4种情况产生：(1)ARDTM及时标注了新出现的活动区，而NOAA存在未能及时标注的现象；(2)ARDTM较完整地标注了活动区的演化末期，而NOAA已经提前结束标注；(3)ARDTM标注了一些没有被NOAA标注的活动区，这些活动区的特点均为磁场强度较弱、磁结构特征不太明显，但日面已有黑子产生且持续了3帧以上；(4)一些已具有磁场特征但未形成黑子的活动区被ARDTM误跟踪。这4种情况所占比例分别为52%、10%、32%和6%。这意味着一半以上的情况是ARDTM及时标注了新出现的活动区。反过来，NOAA标注的活动区比交集多了51个，这意味着ARDTM没有标注51个被NOAA标注的活动区。造成此情况的主要原因有：(1)太阳活动区的特征已经在纵向磁图上消失但仍被NOAA标注；(2)ARDTM漏标注一些磁场强度较弱且特征不明显的活动区。两种情况所占比例分别为88%和12%。在这4组数据中，NOAA在D4数据集上及时地跟踪了新活动区和终止跟踪消失的活动区，所以本文方法和NOAA的结果相似度最高。
为了说明以上几种情况，表 3详细列出了D6测试集中前10天的跟踪结果，第2~6列的数据含义与表 2相同，第7列详细列出了ARDTM方法与NOAA标注结果不同的具体情况。比较典型的是AR12185于10月5日出现，但NOAA于7日才开始标注；AR12183左侧有个活动区，于10月5日形成了黑子，一直到14日才消失，但该活动区一直未被NOAA标注；AR12173已于10月3日从日面消失，但NOAA仍做了标注。
| Date | NOAA | ARDTM | ARDTM∩NOAA | NOAA-(ARDTM∩NOAA) | ARDTM-(ARDTM∩NOAA) | Remark |
| 2014-10-01 | 10 | 10 | 9 | 1 | 1 | AR12171 has disappeared in the solar full-disk image, however, it is labeled by NOAA; AR12182 is not labeled by NOAA |
| 2014-10-02 | 9 | 9 | 8 | 1 | 1 | AR12175 has disappeared in the solar full-disk image, however, it is labeled by NOAA; AR12182 is not labeled by NOAA |
| 2014-10-03 | 8 | 10 | 7 | 1 | 3 | AR12173 has disappeared in the solar full-disk image, however, it is labeled by NOAA; AR12182, AR12183 and AR12184 are not labeled by NOAA |
| 2014-10-04 | 10 | 9 | 9 | 1 | 0 | AR12172 has disappeared in the solar full-disk images, however, it is labeled by NOAA |
| 2014-10-05 | 9 | 10 | 8 | 1 | 2 | AR12176 has disappeared in the solar full-disk image, however, it is labeled by NOAA; AR12185 is not labeled by NOAA; An AR is labeled by ARDTM but not by NOAA |
| 2014-10-06 | 9 | 9 | 7 | 2 | 2 | AR12176 and AR12180 have disappeared in the solar full-disk image, however, they are labeled by NOAA; AR12185 is not labeled by NOAA; An AR is labeled by ARDTM but not by NOAA |
| 2014-10-07 | 8 | 10 | 8 | 0 | 2 | AR12186 is not labeled by NOAA; An AR is labeled by ARDTM but not by NOAA |
| 2014-10-08 | 9 | 10 | 9 | 0 | 1 | An AR is labeled by ARDTM but not by NOAA |
| 2014-10-09 | 8 | 8 | 7 | 1 | 1 | AR12181 has disappeared in the solar full-disk image, however, it is labeled by NOAA; An AR is labeled by ARDTM but not by NOAA |
| 2014-10-10 | 7 | 6 | 5 | 2 | 1 | AR12178 and AR12179 have disappeared in the solar full-disk image, however, they are labeled by NOAA; An AR is labeled by ARDTM but not by NOAA |
| Total | 87 | 91 | 77 | 10 | 14 |
ARDTM较好地解决了文[4]提到的一个太阳活动区被误检测为多个，或者多个太阳活动区被误检测为一个的问题。图 3(a)采用了文[4]图11(a)的迈克逊多普勒成像仪同时间数据(2002年5月4日12:48UT)，红色实线框为ARDTM的检测结果，蓝色虚线框为文[4]的结果。其中，1和2两个区域分别为AR09934和AR09936，而文[4]误把他们识别成了一个活动区。(b)为日震及磁场成像仪2010年11月12日的数据。同样，3和4分别为AR11123和AR11121，而文[4]误识别成了一个活动区。(c)为日震及磁场成像仪2010年11月7日的AR11120的图像，ARDTM顺利检测为一个活动区，但文[4]将其检测成了两个活动区。
|
| 图 3 ARDTM的检测结果图 Fig. 3 The detection results of ARDTM |
图 4展示了D6数据集中连续3帧的跟踪结果(a)，(b)和(c)，时间分别为2014年10月20日22:48UT、21日23:00UT和22日22:48UT。图中粗框代表检测到的太阳活动区匹配到了轨迹；虚框表示该区域未被检测或未被跟踪匹配，而采用预测机制填补的结果。跟踪框的左上角数字表示ARDTM对太阳活动区的跟踪编号。通过与对应时间的NOAA对比，ARDTM具有良好的跟踪结果。比如，(1)ARDTM及时捕获到了新的活动区。ARDTM在(a)及时捕获了AR12194，而NOAA在(c)才被标注，ARDTM比NOAA提前了2天。NOAA中AR12183，AR12186，AR12189，AR12190，AR12191，AR12192和AR12195这几个活动区，ARDTM均比NOAA更早检测到并及时开始跟踪。(2)ARDTM及时终止跟踪消失的活动区。ARDTM在22日(c)及时终止跟踪AR12189，而NOAA直到26日才终止跟踪该活动区。同样，(a)图的AR12186也是类似的情况。(3)预测机制较好地弥补了检测算法漏检带来的漏跟踪或误跟踪问题。ARDTM在(a)和(b)都通过预测机制连续跟踪了AR12189。(4)ARDTM检测并跟踪到了NOAA漏检的活动区。(a)图中编号为23和29的活动区和(b)中编号为33的活动区，在演化过程中都形成了黑子，生命期分别为8、6和6帧，但都未被NOAA标注。
|
| 图 4 ARDTM在D6数据集中连续3帧的跟踪结果。(a) 2014年10月20日22：48UT；(b) 2014年10月21日23：00UT；(c) 2014年10月22日22：48UT Fig. 4 ARDTM tracking results of three consecutive frames in the D6 dataset. (a) Oct. 20, 2014 22:48UT; (b) Oct.21, 2014 23:00UT; (c) Oct.22, 2014 22:48UT |
图 5展示了图 4中AR12194以及D4数据集中AR09093两个活动区的完整跟踪结果。(a)到(l)子图是AR12194从出现到消失的跟踪过程，时间从2014年10月20日22:48:00UT到31日19:48:00UT。该活动区穿越整个日面，磁场结构和面积变化都非常大。(n)到(t)子图是AR09093从出现到消失的跟踪结果，时间从2000年7月16日12:48:00UT到22日12:48:00UT。(m)和(u)是AR09093产生的前一帧和消失的后一帧对应位置图。该活动区演化周期时间较短，特征变化也很大。ARDTM及时捕获且及时停止跟踪该活动区。其中，(t)时刻通过预测机制填补了漏检，(u)时刻及时停止了跟踪。图 6显示了图 5中AR09093和AR12194两个活动区的磁通量和面积的变化曲线(均已进行了投影改正)，可以看出其磁通量和面积都具有较强的渐变性。
|
| 图 5 (a) 到(l)是D6数据集中AR12194完整的跟踪结果。(n)到(t)是D4数据集中AR09093完整的跟踪结果 Fig. 5 Panels from (a) to (l) are the complete tracking results of the number AR12194 in D6 dataset. Panels from (n) to (t) are the complete tracking results of the number AR09093 in D4 dataset |
本文基于深度学习提出了一种太阳活动区检测与跟踪的方法ARDTM。该方法首先使用YOLOv3-spp网络检测太阳活动区，然后结合较差自转改进DeepSort算法跟踪太阳活动区。经过网络训练后，分别用4组数据进行了测试，并与NOAA的结果进行了比较，ARDTM跟踪标注的结果与NOAA的交集约为NOAA的92%。
本文的训练集只采用了日震及磁场成像仪数据，但无论是测试日震及磁场成像仪的数据，还是迈克逊多普勒成像仪的数据，都能很好地检测和跟踪太阳活动区，表现了较好的泛化性。为了与SolarMonitor的结果进行逐一比较，本文主要采用与SolarMonitor网站公布的数据时间基本吻合的数据进行测试，时间间隔约为24小时。本文还对不同时间间隔的日震及磁场成像仪序列图像进行了测试，时间间隔选取为12分钟、1小时、4小时……，一直到48小时。总体来说，间隔为12分钟的数据由于活动区特征变化不大，跟踪结果最好；一直到24小时的时间间隔跟踪效果都还不错；超过24小时的时间间隔就会出现较多的误跟踪现象。总的来说，该方法较好地解决了传统方法中一个太阳活动区被误检测为多个，或者多个相邻的太阳活动区被误检测为一个的问题。而且，该方法能及时捕获新出现的活动区，也能及时终止跟踪消失的活动区，并且通过预测机制弥补检测算法漏检带来的误跟踪或漏跟踪问题。当然，该方法还存在一定的问题，比如由于单一地根据磁场结构的特征检测和跟踪活动区，导致一些具有磁场结构但未形成黑子的特征被误跟踪情况。在未来的工作中，考虑结合其他波段观测的图像进一步提高算法的准确率。
| [1] | QAHWAJI R, COLAK T. Automatic detection and verification of solar features[J]. International Journal of Imaging Systems & Technology, 2005, 15(4): 199–210. |
| [2] | ZHANG J, WANG Y, LIU Y. Statistical properties of solar active regions obtained from an automatic detection system and the computational biases[J]. The Astrophysical Journal, 2010, 723(2): 1006–1018. |
| [3] | CABALLERO C, ARANDA M C. Automatic tracking of active regions and detection of solar flares in solar EUV images[J]. Solar Physics, 2014, 289(5): 1643–1661. |
| [4] | HIGGINS P A, GALLAGHER P T, MCATEER R T J, et al. Solar magnetic feature detection and tracking for space weather monitoring[J]. Advances in Space Research, 2011, 47(12): 2105–2117. |
| [5] | TURMON M, PAP J M, MUKHTAR S. Statistical pattern recognition for labeling solar active regions:application to SOHO/MDI imagery[J]. The Astrophysical Journal, 2002, 568(1): 396–407. |
| [6] | ZHARKOVA V V, ABOUDARHAM J, ZHARKOV S, et al. Solar feature catalogues in EGSO[J]. Solar Physics, 2005, 228: 361–375. |
| [7] | MCATEER R T J, GALLAGHER P T, IRELAND J, et al. Automated boundary-extraction and region-growing techniques applied to solar magnetograms[J]. Solar Physics, 2005, 228(1/2): 55–66. |
| [8] | COLAK T, QAHWAJI R. Automated mcIntosh-based classification of sunspot groups using MDI images[J]. Solar Physics, 2008, 248(2): 277–296. |
| [9] | COLAK T, QAHWAJI R. Automated solar activity prediction:a hybrid computer platform using machine learning and solar imaging for automated prediction of solar flares[J]. Space Weather, 2009, 7: S06001-1–S06001-12. |
| [10] | LABONTE B J, GEORGOULIS M K, RUST D M. Survey of magnetic helicityinjection in regions producing X-class flares[J]. The Astrophysical Journal, 2008, 671(1): 955–963. |
| [11] | GALLAGHER P T, MOON Y J, WANG H. Active-region monitoring and flare forecasting-I. data processing and first results[J]. Solar Physics, 2002, 209(1): 171–183. |
| [12] | LECUN Y, BENGIO Y, HINTON G. Deep learning[J]. Nature, 2015, 521: 436–444. |
| [13] | GIRSHICK R, DONAHUE J, DARRELLAND T, et al. Rich feature hierarchies for object detection and semantic segmentation[C]//2014 IEEE Conference on Computer Vision and Pattern Recognition. 2014. |
| [14] | REN S, HE K, GIRSHICK R, et al. Faster R-CNN:towards real-time object detection with region proposal networks[J]. IEEE Transactions on Pattern Analysis & Machine Intelligence, 2015, 39(6): 1137–1149. |
| [15] | REDMON J, DIVVALA S, GIRSHICK R, et al. You only look once: unified, real-time object detection[C]//IEEE Conference on Computer Vision and Pattern Recognition. 2015: 779-788. |
| [16] | REDMON J, FARHADI A. YOLOv3: an incremental improvement[C]//IEEE Conference on Computer Vision and Pattern Recognition. 2018: 89-95. |
| [17] | SANCHEZ-MATILLA R, POIESI F, CAVALLARO A. Online multi-target tracking with strong and weak detections[M]//Computer Vision-ECCV 2016 Workshops. 2016: 84-99. |
| [18] | CHU Q, OUYANG W, LI H S, et al. Online multi-object tracking using CNN-based single object tracker with spatial-temporal attention mechanism[C]//IEEE International Conference on Computer Vision. 2017: 4836-4845. |
| [19] | WOJKE N, BEWLEY A, PAULUS D. Simple online and realtime tracking with a deep association metric[C]//IEEE International Conference on Image Processing. 2017: 3645-3649. |
| [20] | 付小娜, 廖成武, 白先勇, 等. 基于LeNet-5卷积神经网络的太阳黑子检测方法[J]. 天文研究与技术, 2018, 15(3): 87–93 |
| [21] | 刘辉, 季凯帆, 金振宇. 机器学习在太阳物理中的应用[J]. 中国科学:物理学力学天文学, 2019, 49(10): 105–117 |
| [22] | 崔顺, 许允飞, 苏丽颖, 等. 基于卷积神经网络的全天空地基云图分类研究[J]. 天文研究与技术, 2019, 16(2): 98–108 |
| [23] | HE K, ZHANG X, REN S, et al. Spatial pyramid pooling in deep convolutional networks for visual recognition[J]. IEEE Transactions on Pattern Analysis & Machine Intelligence, 2014, 37(9): 1904–1919. |
| [24] | HOWARD R. The rotation of the sun[J]. Reviews of Geophysics, 1978, 16(4): 721–732. |
| [25] | HOWARD R F, HARVEY J W, FORGACH S. Solar surface velocity fields determined from small magnetic features[J]. Solar Physics, 1990, 130(1/2): 295–311. |
| [26] | KUHN H W. The hungarian method for the assignment problem[J]. Naval Research Logistics, 1955, 2(1/2): 83–97. |