Abstract
Aiming at the simultaneous measurement of the size, shape, and fall velocity of precipitation particles in the natural environment, we present here a new ground-based precipitation microphysical characteristics sensor (PMCS) based on the particle imaging velocimetry technology. The PMCS can capture autocorrelated images of precipitation particles by double-exposure in one frame, by which the size, axis ratio, and fall velocity of precipitation particles can be calculated. The PMCS is calibrated by a series of glass balls with certain diameters under varying light conditions, and a self-adaptive threshold method is proposed. The shape, axis ratio, and fall velocity of raindrops were calculated and discussed based on the field measurement results of PMCS. The typical shape of large raindrop is an oblate ellipsoid, the axis ratio of raindrops decreases linearly with the diameter, the fall velocity of raindrops approaches its asymptote, and the above observed results are in good agreement with the empirical models; the synchronous observation of a PMCS, an OTT PARSIVEL disdrometer, and a rain gauge shows that the PMCS is able to measure the rain intensity, accumulated rainfall, and drop size distribution with high accuracy. These results have validated the performance of PMCS.
1. Introduction
The size, shape, and fall velocity of precipitation particles and their spatial distributions are important microphysical parameters playing key roles in many branches of research, such as quantitative precipitation estimation (QPE) by weather radar [1], assessment of electromagnetic wave propagation in troposphere [2], study on kinetic energy of raindrops and soil erosion [3], validation of simulated microphysics in numerical weather models [4], and parameterization of washout efficiency of particle-bound radionuclides and atmospheric pollutants in general [5, 6]. Due to the force balance on a freely falling drop surface between the gravity, aerodynamic pressure, surface tension, and internal hydrostatic pressure [7], the size, shape, axis ratio, terminal velocity, phase, oscillation, and orientation of hydrometeors show complicated distribution, which varies with the synoptic systems, cloud types, weather conditions, and wind fields of different regions and seasons. These microphysical properties of hydrometeor still remain unclear.
Since the past century, there have been a number of different methods and instruments for measuring hydrometeors [8]. Impact disdrometer uses either acoustic or displacement methods to measure the raindrop size distribution; however, it cannot obtain the accurate measurement of velocity, kinetic energy, and drop shape. The accuracy of acoustic disdrometer is limited by the nonuniform acoustic response over the entire diaphragm; the displacement disdrometer underestimates the number of small drops in heavy rain due to the “dead time”; therefore difficulties remained especially during high rainfall intensity measurement. Optical disdrometer uses nonintrusive scattering or imaging methods to measure the drop size distribution of hydrometeors, and the scattering disdrometer estimates the fall velocity of particles by the empirical assumption of raindrops’ shapes [9], but the errors remain due to the difference between assumption and real rainfall; the imaging disdrometer can obtain the images and shapes of drops, by which the type of hydrometeors can be identified, and the fall velocity of hydrometeors can be obtained by two line-scan cameras [10], or autocorrelation images from a double-pulse light source [11]. However, the air movement and wind around sampling area might influence drops trajectory, which influence the accuracy of particles spatial distribution and precipitation intensity. Although several correction algorithms have been proposed, the errors remain.
Aiming at the simultaneous measurement of the size, shape, and fall velocity of precipitation particles, we present a precipitation microphysical characteristics sensor (PMCS) based on particle image velocimetry techniques in this paper; the critical parts of PMCS are a double-pulse light source and a planar array CMOS camera. The prototype of PMCS was called Video Precipitation Sensor (VPS) [11]; its performance has been validated by field experiments; however, the free fall of precipitation particles is easily affected by the instrument structure; especially in the strong wind environment, the turbulence of air flow, the splash of raindrops on the surface, and the complex variation of ambient light can induce considerable error that cannot be negligible. To address these problems, a new instrument structure which can diminish the turbulence of air flow and the splash of raindrops is designed, and an adaptive calibration method is proposed to reduce the errors of imaging. The following sections contain descriptions of the instrument PMCS, of the image processing and calibration method, and of the field experiments results. The last part summarizes the main features and future work of PMCS.
2. Precipitation Microphysical Characteristics Sensor
2.1. Main Components and Measurement Principle
The PMCS consists of four units: optical unit, imaging unit, acquisition and control unit, and data processing unit, as shown in Figure 1(a). The optical unit (OU) contains a light-emitting diode (LED), a multimode fiber cluster, expanded beam lens, and concentration lens, which can provide a parallel cylindrical light beam for imaging unit; the imaging unit (IU) contains one complementary metal oxide semiconductor (CMOS) image sensor and driving circuit, which record digital images with 640 pixels × 480 pixels and 50 frames per second (fps), and the sampling volume is 300 mm × 40 mm × 30 mm; hence the pixel size of an image is 0.0625 mm2; the acquisition and control unit contains a digital signal processor (DSP) and field programmable gate array (FPGA), which generates and outputs timing signals of double exposures in one frame, controls the optical unit and imaging unit, records and preprocesses the raw images, and encodes and transmits the raw data; the data processing unit (DPU) is a PC terminal that communicates with the acquisition and control unit using a coaxial network cable via the Internet TCP/IP protocol. Software running on the terminal receives the data obtained by the sensor, processes the hydrometeors’ images, and calculates the size, velocity, and shape of hydrometeors. A photography of the PMCS is as shown in Figure 1(b), and the optical unit and imaging unit are integrated into one tunnel housing, on which there are a couple of metal splash elimination grids, to prevent raindrops splashing into the sampling area and to minimize the wind disturbances about the instruments.
(a)
(b)
The double-exposure in one frame (DEOF) of the PMCS plays a key rule on the simultaneous measurement of the size, shape, and fall velocity of precipitation particles, as shown in Figure 2. The CMOS camera runs in 50 frames per second, a pulse synchronous signal generator is used to generate two exposures from pulse light source in each frame, and the single exposure time is 20 us, so blurry that the particle motion is insignificant; besides, the interval of two exposures is 2 ms. During the period of exposures in a single frame, the CMOS camera captures the double-exposure image of the particles crossing the space, by which the particle shape information can be obtained, and the size, axis ratio, and canting angle of raindrops can be calculated; considering the displacement and time internal, its velocity is calculated; further, the velocity and size distribution of particles can be obtained by the time integration method; furthermore, the precipitation intensity and accumulated rainfall can be calculated.
2.2. Image Processing and Feature Extraction
In an attempt to acquire the effective particle images and improve the processing speed of the PMCS, we have proposed an image detection method based on dimension-reduced processing: firstly, every region with 5 pixel × 5 pixel in raw image was treated as 1 pixel, and the image resolution decrease as 128 pixel × 96 pixel; secondly, the particle images in single frame were traversed and labeled by the digital signal processor (DSP), and the distribution of exposure energy is obtained according to its gray information; then the dimension-reduced images are mapped to the raw images; the single frame image without exposure particles is ignored by PMCS; for the single frame image with exposure particles, the location, area, and gray value of each particle are labeled, and its centroid coordinates are calculated. Finally, the particles are paired based on the area and gray value of each particle image in single frame. A diagram of the dimension-reduced processing is shown in Figure 3.
The core of the PMCS is a planar CMOS image sensor, but the CMOS has the characteristics of pixel nonuniformity, voltage instability, light refraction, reflection, scattering, and absorption, which would cause noise and defocusing blur in the image (Figure 4(a)). As shown in Figure 4, the image processing can be divided into 5 steps as follows: in step 1, it extracts valid particle according to the pixels 4-connected principle; in step 2, a median filter algorithm is used to remove the noise; in steps 3 and 4, the point spread function and erosion are applied to restore the filtered images, respectively; finally, an appropriate threshold is selected to obtain the binary images.
Tremendous information of the particles can be obtained from the captured digital images. As shown in Figure 5, the minimal bounding rectangle (MBR) is applied to the feature detection of the particles, and the detailed algorithm is summarized:(1): the short axis of MBR of binary image(2): the long axis of MBR of binary image(3)Perimeter: the summation of pixel values on the contour of binary image(4)Area: the number of pixels, where the value equals 1(5)Equivalent diameter: the diameter of spherical particle with the same volume of nonspherical particle(6)Axis ratio: the ratio of height and width(7)Canting angle: representing the angle between the vertical direction and symmetry axis
The vertical and horizontal velocity can be obtained by the displacement and exposure interval, as shown in Figure 6
Given that only the particles passing through the sampling space during the period of exposure time can be captured, the complete capturing of precipitation particles has a certain probability, and the probability that COMS sensor can capture the fully double-exposure images of precipitation particles iswhere is the height of sampling space, is the vertical size of particles, is the fall velocity of particles, which is related to the equivalent diameter, represents the exposure interval, and is the time of single frame. The velocity of particles varies with its type and size, so a suitable exposure interval should be considered to ensure that the particles passing through the sampling space might be fully photographed twice in a single frame; in this paper, is 2 mms, which is suitable for precipitation of different types (rain, snow, and hailstone). As shown in Figure 7, when the diameter of raindrop is larger than 7.8 mm, the capture probability is below 0%.
Considering the influence of the capture probability on the sampling results, the drop size distribution (denoted by DSD) collected directly by PMCS should be corrected:where represents the DSD collected directly by PMCS, represents the effective sampling space, and represents the capture probability, which are calculated according to the velocity of particle.
Finally, the rain intensity can be calculated by the following integration formula:
The raindrop has an irregular shape with a slightly serrated border caused by the finite resolution of PMCS; therefore an appropriate contour smoothing procedure is used to reconstruct the real shape of each individual drop. As shown in Figure 8, the inner and outer corners were detected, and a transform from Cartesian coordinates to Polar coordinate is conducted in case of distortion of the shape during smoothing; then the contour smoothing procedure is conducted by the Akima Interpolation method [13], and the average contour between inner and outer smoothed contour is calculated; finally, the real shape is obtained by a transform from polar coordinate to Cartesian coordinates.
3. Calibration
The influence of divergence and ambient light is higher in the field observation than in the laboratory. To address this issue, we need to calibrate the instrument again to optimize the image processing. As shown in Figure 9, small glass balls with diameters from 0.5 to 5 mm were dropped through the different location of sampling space in sunny, cloudy, and night conditions. For each size (0.5 mm, 0.8 mm, 1.0 mm, 1.2 mm, 1.5 mm, 2.0 mm, 2.38 mm, 2.5 mm, 3.0 mm, 3.5 mm, 4.0 mm, 4.76 mm, and 5.0 mm), calibration spheres were released from the same height above the center location of sampling space in sunny, cloudy, and night conditions, and at least 5 complete images were collected for each size. For each location (0~5 cm, 5~10 cm, 10~15 cm, 15~20 cm, and 20~25 cm from the light source), calibration spheres (1 mm, 2 mm, 3 mm, and 4 mm) were released in the cloudy condition.
Now, we assume the radius of median filtering is , the radius of point spread function (PSF) is , and the radius of erosion is , and the binary threshold is defined as follows:where represents the brightness of image and is the max value of . The axis ratio and equivalent diameter are used to quantify the shape of drops. The axis ratio and equivalent diameter of the calibration sphere are 1 and , respectively. Hence the object of image processing can be defined as follows: where is the object value, is the weight factor, is the long axis, is the short axis, and .
The image descriptors represent the characteristics of gray image, and the gray level histogram which represents intensity distribution of gray image is given by the following formula:
represents the image function of two variables and , and , where and are image resolution. This function takes the values , where represents the total number of intensity levels of an image. In this research gray-scale images are concerned; therefore equals 256; then a set of image descriptors can be calculated based on the following formula:
The mean brightness value can quantify the shadowing level caused by particles, and the pixel area represents the size of the shadowing area by particles, and both of them are important descriptors of the particle image, which could be applied in the calibration scheme based on the relationship between image features and image processing parameters.
From the above analysis, we can design the calibration algorithm, as shown in Figure 10. Firstly, the calibration algorithm extracts the valid images from the original images. Secondly, the image features are calculated from valid images and the optimized parameters of image processing are obtained from optimized processing with enumeration method where all the ranges , , and are [] pixels and the range of the is from 0.0001 to 1. Then, the relationship between optimized parameters and image features can be determined. Finally, the calibration scheme could be checked by the other calibration images.
A total number of 2116 particle images were collected in the condition of different ambient light and 932 particle images in the condition of different sampling location. The observations are divided into two parts: half of them are used as an analysis set to determine the image processing parameters, while the remaining observations are kept for validation. The results are shown in Table 1. For the radius of median filtering and erosion, 3 pixels have highest probability of 49.46% and 89.51%, respectively. As for the radium of PSF modeling, [] have been chosen with similar probability and 7 pixels mode has highest probability of 27.84%. In fact, none of the three image processing models can completely diminish the blur in images. Therefore, the appropriate threshold value plays a key role in the image processing.
So we discuss the correlation among the image descriptors and optimize binary threshold, as shown in Table 2; significantly negative correlation was found among the descriptors and the binary threshold, while it could be influenced by ambient light and sampling location; an additional analysis should be discussed.
As shown in Figure 11, the calibration experiment in different ambient light condition shows that the average and range of and increase with the diameter of particles; besides, the mean measured in night condition is higher than those measured in cloudy condition, and the mean measured in sunny condition is the lowest. This means that the mean brightness value increases with the light intensity decreasing. However, there are no significant changes in for ; one possibility is that though ambient light could affect the exposure brightness, it has little influence on the shielding area in the COMS sensor caused by particles.
(a)
(b)
While the calibration experiment in different sampling location indicates that the mean decreases with the distance between the sampling location and light source decreasing, as shown in Figure 12, one explanation could be that more ambient light integral in the sensitive surface of CMOS image sensor causes a lower imaging luminance value, and it is especially significant in the location of and . The mean , by contrast, increases with the distance between the sampling location and light source decreasing; this might occur due to the existence of divergence angle of the light source, and the closer they are to the light source, the larger the image area will be recorded for particles with the same diameter.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Since the effect from the ambient is minor, is a more suitable descriptor for the establishment of the self-adaptive threshold than , and to assuage the influence from the sampling location, a fitter relationship between the threshold and pixel area of particles captured in different ambient condition and sampling location is established with the least square method; as shown in Figure 13, the scattered points of and thresh uniformly distributed around the fitted curve, and its root-mean-square-error (RMSE) is 0.076; the correlation coefficient (R2) is 0.8988, since a valuable self-adaptive threshold method based on the image pixel area is established.
The evaluation of the accuracy of the new calibration method is conducted, as shown in Figure 14; the diameters and axis measured by the new calibration method are closer to the real diameter of the small ball; besides, the standard deviation and mean absolute deviation for the diameter and axis obtained by new method are lower compared with the value obtained before calibration, particularly for small-scale balls, since the self-adaptive threshold based on the image pixel area is conductive to the accuracy and applicability of the image algorithm.
(a)
(b)
4. Filed Observations
A PMCS was deployed in the meteorological observation station of Nanjing, China, during 2015, and several rainfall events have been recorded. In this paper, we collect the precipitation data from seven rainfall events in June 2015, and 1643551 images were obtained.
Some anomalous particles should be rejected from the images, and these particles appeared as discrete pixel, incomplete capture, and complex edge profile caused by the spot in the surface of sensor, captured in the edge of sampling space and some insects. We define some descriptors which are useful to characterize these anomalous particles.(a)Dispersion: , the dispersion is a descriptor that reflects dispersion of the pixels. is the total brightness of the binary image (after image processing), and is the number of the raw image’s brightness (greater than 5) in the region where the brightness is equal to 1 in the binary image.(b)Edge pixels: , the edge pixels are the number of brightness in the edge of binary image and show the integrity of the particles captured by PMCS.(c)Shape parameter: , the shape parameter describes the circularity of particle. represent the pixels of the edge profile of the particles, and is the total number of pixels of the binary image.
The , , and can be calculated according to the valid rain images selected manually by their visualizing particle shapes with a diameter interval of 0.2 mm. For each interval, 100 images were obtained. The ranges of are 0.9 to 1, and for all of the valid images, are less than 1; hence 0.9 and 1 are the fixed thresholds of and . Because of the effect of the noise to the image of particles with different diameter, a fixed threshold is inapplicable to the . So we chose images data from one rainfall event; for each interval, its mean and standard deviation after being rejected by and were calculated. Then a fitter relationship between and diameter is established with the least square method, and any particles beyond from were removed, by which a valid range of the is obtained. As shown in Figure 15, the mean of raw images is 1.699, and the mean standard deviation is 0.076; almost all of the calculated from valid images (red dots) lie between the upper boundary (green line) and lower boundary (blue line) of the fitted curves obtained from raw images (black dots) rejected by and ; besides, most of the ambiguous dots are distributed in the top of upper boundary. Finally, the criterion of valid image can be obtained.
The microcharacteristics of raindrops recorded by PMCS were listed in Table 3, we obtained 1203393 images after elimination by the criterion of valid image, and the ratio between the valid images and raw images is 73.2%. These valid images were reprocessed with the image processing algorithm, then the binary images were obtained, and the shape, axis ratio, and fall velocity of raindrops were analyzed and discussed.
The raindrops after undergoing contour smoothing were divided into various diameters internal with a width of 0.2 mm, and the probability contours [14] of these raindrops were discussed. Figure 16 shows the probability contours (in log scale) for equivalent drop diameters in the range (a) 0.8~1.0 mm, (b) 1.4~1.6 mm, (c) 2.4~2.6 mm, and (d) 3.4~3.6 mm. The black lines superimposing on the images are the results from TB model [15] (solid line) and oblate model [16] (dashed line). The probability contours agree closely with the empirical models. The 0.8–1.0 mm raindrop contours are ideally spherical; its probability contours agree closely with the empirical models; the 1.4~1.6 mm raindrop contours are basically elliptical, which is slightly flatter than the oblate model; the 2.4~2.6 mm raindrop contours are more flattened than the 1.4~1.6 mm raindrop contours, with a smaller base than the Thurai model; the 3.4~3.6 mm raindrop contours agree closely with the Thurai model, with a flatter base and a more raised top. The width of the probability contours indicates the drop oscillations; such phenomenon was quantified according to the parameters of axis distributions in Thurai and Bringi [15].
(a) 0.8~1.0 mm
(b) 1.4~1.6 mm
(c) 2.4~2.6 mm
(d) 3.4~3.6 mm
The axis ratio is a suitable parameter to qualify the shape of raindrops, and its distribution reflects the varied processes of raindrop shape. For each raindrop, its axis ratio is calculated, then the axis ratio distributions considering all the internals were obtained, and in an attempt to evaluate the accuracy of the results, the present result was compared with the results from BC model [17], Goddard model [18], Brandes model [19], and TB model [17]. Figure 17(a) shows the variation of axis ratio with the equivalent diameter, the measured mean axis ratio (black circles) and standard deviation ( vertical bars) decrease with equivalent diameter, and the black curve represents a fitted relationship between the measured axis ratio and equivalent diameter by using 4th degree polynomial fitting method for mm:The RMSE is 0.0046, and the coefficient is 0.9983. This formula agrees closely with the empirical models. Figure 17(b) shows the variances among the measured mean axis ratio with the fitted curve and the empirical curves were insignificant, with a finite range of ±0.04, the measured mean axis ratios of the 1.0 mm–2.2 mm drops were slightly larger than the results from BC model and smaller than the results from others models, the axis ratios of the 2.2 mm–4.8 mm drops were slightly larger than the results from empirical models, which were closer to the results from TB model, and the peak values were measured in the internal of 3.4 mm–4.0 mm, which indicates that the measured shape of raindrops in this internal is rounder than the results from empirical models; one reason is that the larger drops in the condition of convective precipitation were affected by environment turbulence, low-level wind shear, and collision processes, by which its oscillation frequencies are likely to increase, and much rounder raindrops would be measured.
(a) Axis ratio distribution
(b) Variances of the axis ratio
For each raindrop, its fall velocity is calculated, and a certain number of outliers were measured, especially for smaller raindrops, in which relatively higher velocities were measured, since we used the following equation to remove these outliers: where is the literature model for fall velocity of raindrops by Atlas et al. [12], and any raindrops crossing the velocity range were rejected. Figure 18 clearly shows that the fall velocity approaches its asymptote (of 9.65 m/s), and the velocity distribution of raindrops obtained by PMCS agrees closely with the Atlas velocity curve (black line), validating the accuracy of PMCS on the fall velocity measurement.
To evaluate the measured performance of PMCS on the precipitation, we compared the rain rate and accumulated rainfall calculated from the PMCS to colocated OTT PARSIVEL disdrometer (denoted by OTT) and rain gauge (denoted by Gauge) measurements. As shown in Figure 19, the general distribution of rain intensity from PMCS agrees well with the results from the OTT, and some noises were observed in the results of PMCS and OTT, which could be caused by the instruments, and the distribution of rain rate calculated from the rain gauge is flatter than the results from PMCS and OTT. The general distribution of accumulated rainfall among the results of PMCS, rain gauge, and OTT is similar, while the accumulated rainfall growth by PMCS is much closer to the result of rain gauge; the accumulated rainfall recorded by PMCS, OTT, and rain gauge is 426.8 mm, 480.8 mm, and 412.7 mm.
(a)
(b)
The variances of rain intensity among these instruments can be quantified by the following formula:where represents the rain intensity observed by two different instruments, respectively, represents the relative deviation, and represents the absolute deviation. Given that the measurement performance of instruments on the precipitation would be influenced by the rain intensity, the samples were divided into four categories based on the rain intensity recorded by the rain gauge. As shown in Table 4, the between the PMCS and Gauge results is 3.3%, which is smaller than the bias between the OTT and Gauge results, while between the PMCS and Gauge results is moderately larger than between the OTT and Gauge results. In general, observed by different instruments is gradually decreased with the increase of rain intensity, and the correlation coefficients between the OTT and PMCS results are larger than 0.9 except when 5 mmh−1 ≤ < 10 mmh−1.
The variances of accumulated rainfall among these instruments can be quantified by the following formula:where and represent the accumulated rainfalls obtained by two different instruments, respectively. As shown in Table 5, the accumulated rainfall variance between the PMCS and Gauge results is 3.3% (−3.4%), while the variance between the OTT and Gauge results is 14.2% (−16.5%), validating the measurement performance of PMCS.
The raindrop spectrum data obtained by PMCS and OTT are mainly divided into four categories based on the rain intensity level recorded by rain gauge. As shown in Figure 20, the drop size distribution obtained by PMCS is well consistent with that obtained by OTT, especially when < 3 mm; in addition, much smaller particles were observed by PMCS with the increasing of the rain intensity, while much larger particles were observed by OTT, especially when > 3 mm. One possible reason is that the laser attenuation of OTT caused by small particles is insignificant; on the contrary, the small particles can be easily imaged by PMCS; besides, the large raindrops are not spherical, so the diameter recorded by OTT is the long axis of drops, and the capture probability of PMCS is below 0.3 when > 3 mm.
(a)
(b)
(c)
(d)
5. Conclusions
In an attempt at the simultaneous measurement of the precise information of precipitation particles, a ground-based optical instrument, the precipitation microphysical characteristics sensor (PMCS) was proposed and developed, which provides a promising alternative to measuring precipitation. The PMCS can record the double-exposure image of each particle in a single frame, by which the size, equivalent diameter, and axis ratio of raindrops can be calculated and their fall velocity can be calculated according to their displacement and time interval. A set of image processing techniques were applied to the PMCS raw images, by which some noises caused by sampling environment and by instrument were removed and the quality of the images was improved. A calibration experiment was conducted by a series of glass balls with certain diameters in sunny, cloudy, and night conditions and in different sampling locations, based on which a self-adaptive threshold method considering the relationship between the image pixel area and the binary threshold of particle image is discussed; it can ensure the measurement accuracy of raindrop sizes.
The field measurement of PMCS was conducted at Nanjing, China, during 2015, based on more than 1203393 valid images of raindrops, the shape, axis ratio, and fall velocity of raindrops are analyzed and discussed, and the accumulated rainfall by PMCS was compared with the results by rain gauge and OTT PARSIVEL disdrometer. The shapes of raindrops are averaged by contour smoothing procedure, the typical shape of large raindrop is an oblate ellipsoid with a flatter base and a more raised top, and the observed shapes agree well with the Thurai model and oblate model for raindrop shape; the axis ratio of raindrops decreases linearly with the diameter, and the fitted line is close to the empirical relationship from Beard, Goddard, Brandes, and Thurai; the fall velocity of raindrops approaches its asymptote (of 9.65 m/s), and the observed results are in good agreement with the empirical line from Atlas; the general distribution of rain intensity from PMCS agrees well with the results from the OTT, and the accumulated rainfall by PMCS is much closer to the result by a rain gauge than the result by OTT; the drop size distribution obtained by PMCS is well consistent with that obtained by OTT. The above results validate the performance of PMCS, which can be widely applied for interpreting weather radar data, investigating regional precipitation features, and studying the physical process of precipitation.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work is supported by the National Natural Science Foundation of China (Grant nos. 41327003, 41505135, BK20150708, and 41475020). The PMCS was developed with support from Ying EnTe Environment Technique Co. Ltd., Nanjing, China.