Developing an Automated Detection, Tracking and Analysis Method for Solar Filaments Observed by CHASE via Machine Learning

Our filament detection method is a supervised machine learning method, which means we need to provide the data sets with filament labeled. Here we utilize the CHASE full-disk H$\alpha$ spectra to identify and label the filament. It is difficult to identify filaments by spectra manually since there are 118 wavelength points in CHASE H$\alpha$ profile. We employed K-Means algorithm (MacQueen, 1967), an unsupervised clustering algorithm which has been widely applied in spectral classification of solar physics (Viticchié & Sánchez Almeida, 2011; Panos et al., 2018; Asensio Ramos et al., 2023). It requires the initial setting of the number of clusters, denoted as $k$. Subsequently, it will assign $k$ centroids $(\mu_{1},\mu_{2},\cdots,\mu_{k})$ and divide the data into $k$ clusters $(S_{1},S_{2},\cdots,S_{k})$ based on the distance to the centroids. Then, it continually updates the centroids and clusters to minimize the intracluster distance by solving following equation:

where $x$ is the contrast of the spectral intensities of each pixel with the average spectral intensity $I/I_{avg}$.

The K-means method is very sensitive to the uneven radiation intensity distribution across the solar disk, which is adversely affecting spectral classification for labeling filaments. We use a second-order cosine function (Pierce & Slaughter, 1977) to remove the limb darkening:

where $I_{\lambda}^{*}(R)$ is the mean observed radiation intensity at the radius $R$ in wavelength $\lambda$, $R_{S}$ is the radius of the Sun, $cos\theta=\sqrt{1-(R/R_{s})^{2}}$ and $a_{\lambda}$, $b_{\lambda}$, $c_{\lambda}$ are the parameters to be fitted. We used the radiation intensity along the solar equator as the fitting data and applied a least-squares fit to minimize the deviation rate. The radiation intensity along the solar equator is selected since there are few activities and it is longer enough to fit the whole disk. For each wavelength $\lambda$ of the H$\alpha$ profile, we need find the best $a_{\lambda}$, $b_{\lambda}$, $c_{\lambda}$ to minimize the following equation:

where $I_{\lambda}^{*}(p_{i})$ is the observed radiation intensity at point $p_{i}$ of the equator and $R(p_{i})$ is the radius at $p_{i}$. We limited the fitting to 880 pixels (about $0.95R_{s}$) since the pixels alone the edge may vary with the observation. Figure 2(a) and (b) give an example to show the H$\alpha$ line center image before and after limb-darkening removal, respectively.

After a series of trial and error we set the $k=30$, i.e., the K-Means algorithm automated categorizes the spectra into 30 classes, in which they belong to sunspot, plage, filament and so on. Then we manually select the class that most closely matches the filament as the output label candidate. As shown in Figure 2(c) and (d), the K-Means spectral classification could effectively segment filaments from the solar disk. However, small chromospheric fibers can not be removed since they are also cold material and have similar H$\alpha$ spectral profiles to the filaments. Here, we set a area threshold of 64 pixels ( about 69 arcsec²) to sieve the chromospheric fibers. We can also find that the filament classified by K-Means are more mottled with holes as shown in Figure 2(d). We adopt the morphological close operation with round structure element (radius is 5 pixel) to fill these holes. Figure 2(e) and (f) show the results after adjustment. Then the labeled data are used for model training. In addition, we found that the data with morphological close operation can improve the ability of U-Net model to distinguish the sunspot and filament. It is because the input of U-Net model is the normalized H$\alpha$ line center images, which can not distinguish the sunspot and filament only by intensity since they have similar dark appearance.

We collected 120 sets of H$\alpha$ spectral observation from CHASE and constructed our labeled dataset by applying the approaches mentioned above. The data span from December 2022 to July 2023, with a time interval of about 2 days. In a supervised machine learning method, the data usually are divided into training, validation, and testing set, respectively. The validation set was utilized during the training phase to choose models, which can prevent overfitting on the training set. The testing set was employed to assess the performance of trained models. Therefore, these 120 sets are divided into 20 groups, the third and the sixth sets in each group are chosen as the validation and testing sets, respectively, while the others are chosen as the training set, with the ratio being about $1:1:4$. In this way, we can ensured a evenly distribution in our data sets so that the trained model could be robust and increase its generalization performance.

The data augmentation technique is usually applied to enlarge the dataset by flipping and rotating transformations since the convolution kernel does not have rotational symmetry and axis-symmetry, which can reduce the risk of overfitting (Mikołajczyk & Grochowski, 2018). Here, we also adopt data augmentation during the training process to enhance the robustness of our model. In each training iteration, the data are randomly applied transformations, including vertical flipping (with a $50\%$ probability) and rotation with a angle within $[-45^{\circ},45^{\circ}]$ (with an $80\%$ probability).

The U-Net model derives its name from its U-shaped architecture, as shown in Figure 3. Upon receiving image inputs, it initiates an encoding process to progressively extract high-level semantic information while reducing the image dimensions, as shown in left part of Figure 3. Following this, the decoding process reconstructs the feature maps to the original image size (right part of the U-shaped architecture), delivering detailed pixel-to-pixel level prediction results.

In our model, the encoding and decoding parts each has four blocks, which are repeating patterns of layers. During the encoding process, each block has two $3\times 3$ convolution layers, followed by a $2\times 2$ max-pooling layer, as indicated by the blue, red, and black arrows in the left side of Figure 3. A convolution layer works as a feature detector to get the feature maps. As the layers of the neural network go deeper, the number of the output channels gets larger, which can regard as the feature numbers increasing. Max-pooling layers output a feature map containing the most prominent features of the previous feature map, which play a critical role in compressing the size of feature maps while preserving important features and relationships in the input images. For example, the input is an image with 1 channel and dimension size of $2048\times 2048$, the output after the first block becomes a feature map with 64 channels and dimension size of $1024\times 1024$, as shown in Figure 3. The feature maps from the convolution layer passes on to activation functions to make the network adapt to nonlinear features from the previous layers. We use LeakyReLU and Sigmoid functions as our activation functions, which are defined as:

where $\alpha$ is a hyperparameter and set to $0.01$ in our model. The Sigmoid function is only applied in the output layer of the model to transform the output into the range between 0 and 1 for binary classification, enabling the distinction between targets, i.e., filaments and the background.

In the decoding process, the feature and spatial information through a series of transpose convolutions return to the original image dimension. Transpose convolution, often referred to as deconvolution or up-sampling, involves the introduction of zero values through interpolation in the image, followed by a convolution operation, allowing for effective image magnification. After this block, the output feature map size are 4 times larger, from $128\times 128$ to $256\times 256$, $256\times 256$ to $512\times 512$ and so on, as indicated by Figure 3. Copy and crop, or skip connections, indicated by the gray arrows in Figure 3, entails the concatenation of semantic information derived from the encoding process with the corresponding feature maps during the decoding process. This operation facilitates the transmission of augmented semantic information, thereby bolstering the segmentation performance of the model.

Loss functions play an important role in machine learning. They define an objective by which the performance of the model is evaluated. The parameters learned by the model are determined by minimizing a chosen loss function. In other words, loss functions are a measurement of how good our model is at segmenting the filaments. We utilize the Focal Loss function (Lin et al., 2017), a well-established choice in binary classification scenarios:

where $y$ represents the values in the label, and $\hat{y}$ is the probability values output by the model. $\gamma$ is a hyperparameter for reducing the relative loss for well-classified examples, putting more focus on hard, misclassified examples, i.e., filament regions. $w$ stands for the weight of filament regions. Focal Loss proves beneficial in addressing the challenges posed by imbalanced distribution of positive and negative samples, corresponding to the scenario where the filament regions are much smaller than the non-filament regions. The training of the U-Net model involves the iterative adjustment of the weights of various convolution kernels within the convolution layers to minimize the Focal Loss by the optimizer, which we choose ADAM optimizer (Kingma & Ba, 2014) in our model. A dropout layer is a regularization technique employed during model training to stochastically zero out a subset of weights, which can mitigate overfitting in the model. We also adapt dropout layers in our model. The detailed parameter settings are list in Table 1.

The precision $P$ and recall ratio $R$ are common semantic segmentation evaluation strategy which we also adopt in our model evaluation. These two metrics are defined as:

where TP, FP, and FN denote the true positive, false positive, and false negative measurements, respectively. In addition, the intersection over union (IoU) is often used to evaluate the performance of a model (Rezatofighi et al., 2019), which is defined as:

In practical, the IoU over $60\%$ is a good enough score. The average precision and recall ratios of the training, validation, and testing are listed in Table 2. The precision, recall, and IoU ratios of each data file are also plotted in Figure 5. These results indicated our model is a viable strategy for solar filament recognition.

The dynamic evolution of filaments has consistently been a key issue in the study of filament (Chen et al., 2020). CHASE provides the full-disk H$\alpha$ spectral data, which bring us a chance to extract the dynamic information of filaments. We applied the cloud model (Beckers, 1964) to fit the H$\alpha$ spectra of filaments in order to derive its LOS velocities. The contrast function of cloud model refers to:

where $\tau=\tau_{0}e^{-(\lambda-\lambda_{0}-v\lambda_{0}/c)^{2}/w^{2}}$, $I$ represents the spectral profile of the filament region, $I_{0}$ corresponds to the quiet region, $\lambda_{0}$ is the line center wavelength, and $c$ is the light speed. The parameters to be fitted include the LOS velocity of the cloud $v$, optical depth $\tau_{0}$, source function $S$, and Doppler width $w$.

The cloud model inversion is applied to derive the spectral profiles of each individual filament. As shown in Figure 6(a), the blue contour represents the filament detected by our U-Net model. The red contour represents the extension of the blue one by 10 pixels. The average spectrum between the red and green contours is extracted as $I_{0}$. we determine the line center wavelength $\lambda_{0}$ by using the moment method (Yu et al., 2020). Then, the contrast profile of the spectrum in blue contour is applied by a least-squares fitting, providing the inversion result of the spectral profile. Figure 6(b) presents the inversion result of LOS velocity distribution of the entire filament region.

Filament spine defines the skeleton of a filament along its full length, with several extend barbs. In order to derive the filament spine, many authors employed iterations of the morphological thinning and spur removal operations (Fuller et al., 2005; Qu et al., 2005; Hao et al., 2013). However, after iterated spur removal operations the spines often become shorter than the original ones. Bernasconi et al. (2005) used the Euclidean distance method to overcome the shortcoming by morphological spur operation. Yuan et al. (2011) and Hao et al. (2015, 2016) used the algorithm based on the graph theory by calculating the paths from end to end points of filament skeleton in pixels, where the longest path is kept as the spine. This method is also employed in our work.

Following the acquisition of the LOS velocity field distribution of filaments, we want to systematically analyze the evaluations of various filaments. However, different filaments have different length, we need compare them in a standardized manner. Thus, we straighten and normalize the filaments along their spines. This approach can help to uncover unified patterns in the dynamical evolution of various kind of filaments. After extracting the spine of the filament, we could obtain the distribution of the coordinates $\{(x_{i},y_{i})\}_{i=1}^{n}$ of points along the spine with respect to the distance set $\{l_{i}\}_{i=1}^{n}$, where $n$ is the number of pixels along the spine. Subsequently, we perform spline interpolation then obtain the parametric equation of the main spine:

Then, we differentiate the parametric equations and apply Gaussian filter (with $\sigma=5$) for smoothing. With the tangent equation $(x^{\prime},y^{\prime})=(\mathrm{d}X/\mathrm{d}l,\mathrm{d}Y/\mathrm{d}l)$, we can obtain the corresponding normal equation $(y^{\prime},-x^{\prime})$. Subsequently, we shift each point on the major axis in the direction perpendicular to the normal by the size of $2S_{f}/L_{f}$, where $S_{f}$ is the filament area and $L_{f}$ is the spine length, respectively. This process can straighten filaments along their spines according to their irregular shapes. Figure 6 gives an example of our process. The yellow and orange curves shows the derived filament skeleton and spine in Figure 6(c) and (d), respectively. Figure 6(e) shows the straightened filament, which also effectively preserves the morphological characteristics of filaments, such as brab structures.

Moreover, the straightening approach enables a quantitative study of the distribution of physical information along and perpendicular to the spine. As shown in Figure 6(f) and (g), the LOS velocity distribution along and perpendicular to its spine. Note that the majority of filaments exhibit a rhombic or elliptical LOS velocity distribution, suggesting that the locations with the maximum LOS velocity in filaments are typically near the spine.

The CSRT tracker is C++ implementation of the CSR-DCF (Channel and Spatial Reliability of Discriminative Correlation Filter) tracking algorithm (Lukežic et al., 2017) in OpenCV library (Bradski, 2000). The tracked object is localized by summing the responses of the learned correlation filters and weighted by the estimated channel reliability scores. In other words, the CSRT tracker distinguish the target and background based on adjusting the weights of different channels, e.g., 10 HoGs (histogram of oriented gradients) channels and intensity channel for grayscale-images. After that, the target is localized by the probability map and its region is updated. Furthermore, Farhodov et al. (2019) integrated the region-based CNN pre-trained object detection model and the CSRT tracker, and got better tracking results since the detection model has already separated the traget and the background. Therefore, we integrate the U-Net model and the CSRT algorithm, i.e., the filaments detected by the U-Net model are used as inputs of the CSRT tracker.

The third part of Figure 1 shows our tracking processing scheme, which consists of two part, i.e., the initialization step and the update step. During the initialization step, we need to set a series of boxes which contain the detected filaments in each frame since the input of the CSRT tracker is an image with targets with bounding boxes. A simple way is using the minimum boundary boxes of the detected filaments directly as the input. However, filament sometimes may split as several fragments, or several fragments may merge into one filament. We adopt a distance criterion of 50 pixels (about 52^′′) to combine filament fragments and a area criterion of 200 pixels (about 216 arcsec²) to filter relative small filament fragments. These operation can enhance the stability and accuracy of the CSRT tracker. Figure 7 shows two frames of H$\alpha$ line center images observed at 04:16 UT and 05:51 UT on 2023 September 16 as the previous and subsequent frames, respectively. The orange boxes in Figure 7(a) and (b) are the results after the initialization step, where each box indicates the merged single filament. The previous frame has been marked with a unique tracking ID, as shown by the numbers above the boxes in Figure 7(a) and (b). Then the CSRT tracker tracks the subsequent frame and output the tracked filaments based on the previous tracking IDs. Figure 7(c) shows the tracked filament indicated by the blue boxes, which have the same tracking IDs as that in previous frame in Figure 7(b).

However, in addition to the fragmentation of filaments, it is possible that the filament maybe disappear after eruption as well as the formation of a new filament. Therefore, in the update step, we adopt the IoU ratio to compare the tracking result of CSRT tracker with the results of the subsequent frame after initialization step. If the filament (merged boundary boxes) of the subsequent frame after initialization step has the largest IoU with a certain tracked filament by the CSRT tracker, it will be set the same tracking ID; if the IoU is empty, i.e., there is no corresponding filament, it will be set a new ID. In this way, the update step is finished and the result will be input to CSRT tracker again for tracking the next subsequent frame until all frames are tracked. The red boxes in Figure 7(d) are the results after the initialization step, which has no IDs yet. After the update step, the red boxes in Figure 7(d) and blue boxes Figure 7(c) are compared by their IoU ratios and finally output the tracking results indicated by the yellow boxes in Figure 7(e) and (f). The update step can effectively handle with the situation of the eruption of a filament. We extracted a filament with ID No.12 within field-of-view of the blue box in Figure 8(a) as an example. We can see its dynamic evolution and eruption from Figure 8(b) to (p) at different time, and finally disappeared after eruption in Figure 8(q).

We selected 10 groups of CHASE observation for testing the tracking performance. Each group is the first day of a month and has about 15 frames of H$\alpha$ line center images with time interval about one hour. If the filament are tracked by our tracking method in more than 3 frames, it will be counted for the tracking accuracy testing. If the tracking ID changed but the ground-truth is the same filament, it would be regarded as false tracking. The testing results are plotted in Figure 9 and summarized in Table 3. The average tracking accuracy is $81.7\%$, which confirm the good performance of our tracking method.