Earthquake Reaearch in China  2018, Vol. 32 Issue (3): 388-399
Research on Regional Crustal Deformation Characteristics Using Displacement Time Series Data of GNSS Reference Stations in Xinjiang
Li Guirong, Wang Xiaoqiang, Li Jie, Liu Daiqin, Ailixiati Yushan, Chen Li, Li Rui
Earthquake Agency of Xinjiang Uygur Autonomous Region, Urumqi 830011, China
Abstract: Using GNSS data from the Crustal Movement Observation Network of China (CMONOC), and PODAP software which was developed by the Satellite Navigation Institute of Wuhan University, the authors calculated data from 31 GNSS stations from July 1, 2011 to December 31, 2014, sampling at 30 seconds, and studied regional crustal deformation characteristics. Analysis results showed that in southwestern Xinjiang, the NS movement rate was influenced by Indian plate pushing. Under the blocking effect of the Tarim Basin, the EW movement rate was slightly smaller. In the north Tianshan area, the vertical dimension movement rate was quite different, which shows as a high gradient zone in the combination area between the basin and the mountain. With regards to regional overall characteristics, the authors considered that the intersection region from south Tianshan and Kunlun Mountains was prone to strong earthquakes, especially moderate-strong earthquakes, even more than M ≥ 7.0 earthquakes. Middle of North Tianshan was the turning point of the vertical movement rate around Bayanbulak, and was also the high gradient zone of vertical movement. The area is also prone to strong earthquakes in the future.
Key words: GNSS     Displacement     Time series variation     Characteristics

The Global Navigation Satellite System (GNSS for short) incorporates satellite navigation systems such as the US GPS, Russian GLONASS, European GALILEO and the Chinese BeiDou Navigation Satellite System which are being developed. It has made notable progress in the study of crustal movement, tectonic deformation and earthquake prediction using GNSS for earth observations, and many important research results have been obtained. The time series of crustal movements observed at continuous observation stations can clearly reflect the nonlinear process of crustal movement, which can not only get global observation results, but also many regional observation results (Gu Guohua, 2007; Li Jie et al., 2010, 2015a, 2015b; Jiang Zaisen et al., 2003; Wu Zemin et al., 2015; Liu Guangming et al., 2012; Wang Xiaoqiang et al., 2012; Gu Guohua et al., 2002). The analysis of various observed time series is the basis for earthquake analysis and prediction, and displacement time series are the most basic time series in GNSS deformation measuring data (Li Jie et al., 2015a). The study of time series obtained by GNSS observations has become a developing trend of research on crustal movement. The rapid development of GNSS observation systems, a swift increase in the number of observation stations and large increases in sampling rates of receivers have dramatically expanded the quantity of observations, which can help improve the accuracy of crustal movement observations and other measurements. The rate of vertical crustal movement calculated by using years of continuous GNSS observation data can visually reflect ongoing dynamic changes of crustal movement. According to technical analysis indicators in the Initial Design of Crustal Movement Observation Network of China, the accuracy of each component of geocentric coordinates of reference stations is better than 2cm in a single day, and the accuracy of each component of horizontal movement velocity of observation stations is better than 2mm/a. High-precision observation results provide the most convictive method and evidence for the study of change of crustal movement using GNSS observations (Gu Guohua, 2007; Li Jie et al., 2010, 2015a, 2015b; Jiang Zaisen et al., 2003; Wu Zemin et al., 2015; Liu Guangming et al., 2012; Wang Xiaoqiang et al., 2012; Gu Guohua et al., 2002).

1 A BRIEF INTRODUCTION TO DATA COMPUTATION METHODS

The Precise Satellite Orbit Determination and Positioning (PODAP) developed by the Satellite Navigation Institute of Wuhan University is adopted as a data-computing program. PODAP software is based on inertial system, adopts the non-differential processing mode, and uses the single-station data preprocessing method which is consistent with Blewitt, to detect and repair cycle slips and eliminate outliers as much as possible. New fuzzy parameters are introduced for unrepaired cycle slips, and undetected cycle slips and gross errors are processed in quality control of the estimation module, for which the computation method can be found in related studies of Shi Chuang et al.(2009)and Li Guirong et al. (2013).

Displacement time series are processed using formula (1), the simplest linear fitting method, to get linear velocities of three components (NS, EW and vertical) of displacement of each reference station.

 $y = ax + b$ (1)

where, y is the displacement value of each component, a is the linear velocity, x the time and b the initial displacement of each component.

2 AN OVERVIEW OF GNSS DATA

The Crustal Movement Observation Network of China (or CMONOC for short) is a key scientific and engineering project for China (Wang Xiaoqiang et al., 2012). It comprises 6 ministries and commissions coordinated by the China Earthquake Administration during the Eleventh "Five-year Plan" period. A total of 260 reference stations were built and upgraded nationwide, including 31 stations in Xinjiang, among which 28 were newly built and 3 were upgraded. The civil works and equipment installation of the 31 reference stations (Fig. 1) were completed in 2010, and were put into trial operation after completion and acceptance of the project in 2011. With the exception of the Urumqi reference station, whose name was a duplicate of the existing IGS station name and was thusly renamed URU2, where the data recorded during January 4, 2011 - April 19, 2013 doesn't match up with data from April 19, 2013 - December 31, 2014 for the time series analysis, other reference stations have been standardized to record and store GNSS data since January 4, 2011, and have accumulated continuous observation data for nearly 4 years.

 Fig. 1 The Xinjiang reference stations distribution map of the Crustal Movement Observation Network of China

The 31 reference stations of the Xinjiang CMONOC are distributed evenly in seismically active zones to monitor constant changes of crustal movement. There are 4 reference stations in the Altay region, closely spaced in the range of 97km-145km. The middle section of the Tianshan Mountains and the western section of the northern Tianshan Mountains is the main economic belt of Xinjiang, where seismic activities are intensive. To better track and monitor the changes of crustal deformation in this region, a total of 13 reference stations are laid out in the region bounded by Urumqi, Wenquan, Zhaosu and Korla, and the distance from the nearest reference station varies in the range of 98km-207km, among which, spacing between Dushanzi and Shihezi station is the shortest, 98km, and spacing between Karamay and Wenquan station is the largest, 207km, with an average spacing of 143km. The spacing between stations in the western section of the southern Tianshan Mountains is slightly larger, and the average distance from the nearest reference station is 223km. The spacing between stations in the Junggar Basin, eastern Tianshan Mountains, Altun Mountains and Tarim Basin is larger, with an average of more than 250km. The spacing between Jijitaizi and Shanshan reference station is 333km, which is the maximum.

Affected by the surrounding environment, the effectiveness rate of data recorded at Yining reference station is only about 75%. Influenced by the rise of groundwater, data reliability of the Karamay reference station is reduced. Except these two reference stations which have slightly lower data quality, the other 29 reference stations have good data quality. Besides, 21 of the 31 reference stations are based on bedrock, 9 on soil layer, and are all buried under 1.5m, therefore the precision of data is relatively high.

According to the design of the CMONOC, each station uses dual frequency geodesic receivers of Trimble R8 type or above. By setting the sampling rate of the receiver, GPS and GLONASS satellite data sampling at 30s, 1Hz, 10Hz, 20Hz and 50Hz can be obtained. In general, we use data sampling at 30s for crustal deformation monitoring, and use high frequency data sampling at 1Hz, 10Hz or even 50Hz for the study of coseismic deformation. In this study, we mainly calculate GPS observation data sampling at 30s. The GLONASS data is not processed, and computed results are single-day solutions.

For unified coordinate frame systems, we use brdc broadcast ephemeris, igs final precise ephemeris, igs final precise clock bias sampling at 5s and igs final precise clock bias sampling at 30s released by the CDDIS (The Crustal Dynamics Data Information System) all the time. Dynamical ephemeris files such as tables of polar motion and leap second provided by Satellite Navigation Institute of Wuhan University are used. Ephemeris that needs to be updated can be found in references (Li Guirong et al., 2013).

PODAP is a high-precision single point positioning software, and data processing results obtained using this software can meet the demands in crustal deformation monitoring. In this study, data recorded by 31 CMONOC GNSS reference stations sampling at 30s during the period from January 4, 2011 to December 31, 2014 are selected for calculation. By calculating repeated accuracy of a single-day solution, repeated accuracy of data calculations can be obtained. Results in Table 1 can be obtained by statistical analysis. Leaving out gross errors, the maximum RMS values are 6.8mm in EW direction, 4.7mm in NS direction and 8.4mm on vertical dimension, and the minimum RMS values are 0.8mm in EW direction, 0.4mm in NS direction and 4.0mm on vertical dimensions. After eliminating gross errors, the average repeated accuracy is calculated to be 2.4mm in the EW direction, 1.3mm in the NS direction and 4.0mm on the vertical dimension. The repeated accuracy is high, all in millimeters. Values in the EW direction with repeated accuracy is better than the average account for 55.9%, 51.1% in the NS direction, and 55.1% on the vertical dimension. Ratios of excellent and good results are all above 50%. With the exception of the Karamay and Yining reference stations, three-component data results of point displacement of 29 reference stations are analyzed, and thus the current change features of crustal deformation in Xinjiang are acquired, providing reliable data for the study of the relationship between crustal deformation and seismic activity.

Table 1 Statistics of repeatable accuracy of data calculation
3 DATA PROCESSING AND RESULTS

In the Windows system, PODAP software has a friendly visual interface to facilitate data calculators to perform post-treatments. Results of three components of displacement from the time series analysis results file are split up with a self-programming batch program to obtain NS components, EW components and UP component data and such data can be applied directly to general data analysis software for earthquake system, Mapsis. Fig. 2 is a map of NS components from 29 reference stations of Xinjiang CMONOC obtained by the PODAP computing software, drawn using Mapsis to make it easy to compare and comprehend the trend of point movements and to find special time-sequence curves that is different from other reference stations. It also allows us to explore whether the characteristic time of background change tendency corresponds to the origin time of a moderate-strong earthquake, and the same for EW components (Fig. 3) and UP components (Fig. 4).We can see from the maps that the displacement time series of each component of the reference stations present a more regular change in trend, the change of the NS component of each station varies slightly, characterized by gradual decreases from south to north and west to east.The changes of EW components are basically the same, and the changes of vertical components are more stable, showing a certain law of annual variation.

 Fig. 2 NS component displacement time series of GNSS reference stations of Xinjiang "CMONOC"

 Fig. 3 EW component displacement time series of GNSS reference stations in the Xinjiang "CMONOC"

 Fig. 4 UP component displacement time series of GNSS reference stations in the Xinjiang "CMONOC"

As we can see from the displacement-time changing curves of the NS components in Fig. 2, pushed by the Indian plate over time, the Pamir arc structure and its adjacent western section of the southern Tianshan Mountains are strongly influenced by NS stress, thus the displacements of the NS components of reference stations in this area are much larger than that of other regions. Specifically, since 2011 when observation was first carried out, the time series curve of the NS component of the Taxkorgan reference station shows an increment of NS coordinates of about 85mm based on the ITRF frame, while the increment of NS coordinates of the Jijitaizi and Qinghe reference stations are almost 0 based on the ITRF frame. That is, most reference stations in western Xinjiang, especially west of 80°E, have larger cumulative NS displacements, such as Bulungkol, Taxkorgan and Yecheng reference stations along the west Kunlun earthquake belt. However, reference stations in the east of 85°E east longitude, especially in Altay and the eastern section of the northern Tianshan Mountains, have the maximum change of displacement of only 20mm and cumulative variation of about 0mm.Can the above results indicate that the west Kunlun earthquake belt is most affected in the NS direction by the pushing of the Indian plate? A large amount of displacement from the south is gathered in the arc region, and constant displacement accumulation also brings about considerable strain energy, resulting in major deformation changes, and its effect continues to work northward. It has been six years since the Wuqia MS6.8 earthquake in 2008, and in the aftermath of the Yutian MS7.3 earthquake on February 12, 2014, the possibility and urgency of a major earthquake in this area is evident.

It also can be seen in Fig. 2 that the NS components showed accelerated northward movement around December 2011, and for nearly six months after February 2012, the NS movement stalled or even went backwards to the south. In this period of time, the Lop M6.0 earthquake on March 9, 2012, the Xinyuan-Hejing MS6.6 earthquake on June 30, 2012 and the Yutian MS6.2 earthquake on August 12, 2013 took place. In December 2012, accelerated northward movement appeared after fast southward movement, which corresponds to the M6.1 earthquake occurring 43.7km from Kazakhstan on January 29, 2013. The anomalous changes from the end of 2013 to May 2014 corresponds to the Yutian MS7.3 earthquake on February 12, 2014. The turns and changes in trends of curves basically correspond to the occurrence of earthquakes with magnitude over 6.0 in Xinjiang and its adjacent areas, which is with a certain implication.

Fig. 3 shows that sequential variations of EW components are relatively stable. Except for permanent deformation caused by the MS7.3 earthquake on February 12, 2014 at the Yutian reference station, sequential variations of EW components of reference stations in Xinjiang are not much different. We calculate cumulative changes of all reference stations, and get that most of the displacements are about 120mm, that is, based on ITRF frame, the EW components of reference stations have been moving eastward since 2011 and have shifted about 120mm to the east so far. Of course, there are stations with relatively small displacements. Displacements of the two reference stations, Taxkorgan and Bulungkol, are 98mm and 103mm respectively. But at the same time, these two reference stations have the largest displacement of NS components. Therefore, it is believed that the region where Taxkorgan and Bulungkol reference stations are located is affected by strong NS stress and relatively weaker EW stress, and possible faults of stress accumulation and release should be EW-striking thrust faults or NS-striking strike-slip faults. Correspondingly, for the Jijitaizi reference station (located in east section of northern Tianshan Mountains) with larger EW displacement and smaller NS displacement, the stress triggering earthquakes may be EW.

It also can be seen in Fig. 3 that the change rate of EW displacement accelerated since October 2011, and the EW movement began to stagnate by February 2012, and then the overall trend of eastward movement of EW component was restored around November 2012. Three earthquakes with magnitude over 6.0 occurred during this period of time. From August 2013 to May 2014, the variation trend of time series curves of reference stations showed a group turns again, which was manifested as accelerated eastward movement and then accelerated backward movement to the west. This change of state corresponds to the Yutian MS7.3 earthquake on February 12, 2014.

Fig. 4 shows time series curves of UP components. On the whole, it turns to uplift at the beginning of each year and subside in July and August of each year, and there is no significant difference in the amount of uplift and subsidence of most reference stations. Except the Yining and Karamay reference stations, affected by the environment, which have low data quality and show abnormal changes in trend, the variation patterns of UP component time series curves are generally consistent.

The NS and EW components of the three-component time series change curve can show a certain trend change before and after an earthquake, which has an impact on earthquake prediction. However, the accuracy of calculation results of vertical deformation is slightly lower, and it is still impossible to find possible changes related to the occurrence of earthquakes. It may be possible to extract some abnormal information by filtering, fitting or segmentation processing, therefore, further research will be carried out.

Formula (1) is used to calculate linear motion velocity of each component of stations for three components of displacement of 29 reference stations (Table 2). It is found that the velocity of the NS component is the highest at Taxkorgan reference station, 23.08mm/a, and the minimum is at Jijitaizi reference station, -0.06mm/a.A negative value indicates a slight southward movement of the station. The velocity of the EW component is the highest at Wulasitai reference station, with a speed of 32.70mm/a, and the minimum at Taxkorgan reference station, 24.77mm/a. It can be seen that the difference between the maximum and the minimum movement velocity of the EW component is not much at 7.93mm/a. The UP component has the highest rate of rising at Bayanbulak reference station, 1.75mm/a, Qiemo reference station shows the highest rate of subsiding, -2.17mm/a, and the UP component has the minimum vertical change rate at Hotan reference station, which is -0.05mm/a.

Table 2 The linear motion rate of 3 components of GNSS reference stations in the Xinjiang "CMONOC"

By distributing data in the table over Figs. 5, 6 and 7, we can clearly see that the distribution law of velocity field of the GNSS reference stations in Xinjiang, which is also the overall trend of regional movement change in Xinjiang. In Fig. 5, the isoline map of movement velocities of NS components, shows that the southwest of Xinjiang is pushed by the Indian plate, thus the change rate of its SN movement is larger, while the linear motion rate of NS components in the Altay earthquake belt and the east section of northern Tianshan Mountains is small, that is, from the southwest of Xinjiang to the northeast, the change of NS movement decreases gradually, showing a smooth transition.

 Fig. 5 Average velocity field distribution of NS component in reference stations

 Fig. 6 Average velocity field distribution of EW component in reference stations

 Fig. 7 Average velocity field distribution of UP component in reference stations

The isoline map of movement velocities of EW components (Fig. 6) also presents a gradual change from southwest to northeast, but unlike the NS components, it shows a change trend of progressive increase from southwest to northeast, but in the middle section of Tianshan Mountains, different from the stair-step smooth transition of NS components, the average movement velocity distribution of EW components is affected by tectonics here, and a distinct closed region is formed. The region from the northern Tianshan Mountains to Altay also follows the law of gradual change, however, it is not consistent with the change trend of isolines of movement velocity of EW components, a change trend of gradual decrease is observed in this region. Therefore, it is believed that stress effects on the middle section of the Tianshan Mountains may be quite complex, which leads to the formation of a symmetrical pattern of northern and southern foothills of the Tianshan Mountains, centering on a large transform fault, the Bolokenu-Aqikekuduke fault. Since 2011, moderate-strong earthquakes have mainly occurred in regions from Dushanzi to Xinyuan, from Shihezi to Wulasitai and from Urumqi to Wulasitai where movement velocity of EW components varies greatly. Based on the same principle, we believe this may also be one of the reasons for strong seismic activities in the region from Taxkorgan to Wuqia.

The results of vertical motion change rate (Fig. 7) shows that there are rapid upward movement changes in Bayanbulak and Zhaosu, while at Kuqa reference station, which is close to Byanbulak, its vertical displacement change rate shows fast descending movement. There are also obvious high gradient zones of velocity changes between Urumqi and Shihezi, Urumqi and Wulasitai. It can be seen from the map that the interior of the Tianshan Mountains is rising up, and the Junggar Basin and Tarim Basin are descending. Relative motion between adjacent blocks makes the mountain uplift more obvious and the orogenic movement of Tianshan Mountains results in large-gradient change of vertical movement in the basin-mountain junction zone. Wang Xiaoqiang et al. (2009) believed that high gradient zones of vertical velocities are areas with intense vertical shearing activities, regions with large vertical shear strain energy and also regions with strong crustal tectonic activities. High gradient zones of vertical movement velocities are also favorable places for earthquakes. Most earthquakes occur at the intersection of different high gradient zones or at the turning part of high gradient zones. In the area covered by leveling networks, from 1950 to 2008, 34 MSn earthquakes occurred in Xinjiang. Except 3 earthquakes in the eastern Kunlun Mountains and 1 earthquake in the Altay Mountains, which occurred just beyond the edge of high gradient zones, the other 17 earthquakes occurred within or near high gradient zones. The vertical movement velocity gradient zones have a good control of contemporary seismology (Wang Xiaoqiang et al., 2009). According to former research results, we believe that this distribution pattern may have existed for a long time. Based on this, analysis suggests that transition zones of vertical movement changes or high gradient zones of vertical movement changes are advantaged areas for future earthquakes.

It can be clearly seen from the three isoline maps of average movement velocities of three components that the overall stress distribution in Xinjiang is strong in the southwest and weak in the northeast. The overall characteristics of stress distribution in this region results in a large amount of energy accumulating in the joint area of the southern Tianshan Mountains and Kunlun Mountains, forming an active zone where moderate-strong earthquakes or even major earthquakes with magnitudes over 7.0 occur. Since 1970, there have been two major earthquakes with magnitudes over 7.0 in this region, the northwest Akto MS7.3 earthquake on August 11, 1974 and the Wuqia MS7.1 earthquake on August 23, 1985. Therefore, it is concluded that there will be a great risk of M≥7.0 earthquakes in this region in the next 10 years or 15 years.

4 CONCLUSION AND DISCUSSION

(1) The above analysis shows from different aspects that from the southwest of Xinjiang to the northeast, the NS components of reference stations present a trend of gradual decrease, and the EW components show an arch-shaped change trend of progressive increase first, gradual decrease and then increase again. However, for the movement changes of UP components, although the results are different from former research results, the pan-like distribution(Liu Guangming et al., 2012), the trend remains broadly consistent.

(2) When a turning of trend occurs to the time series curves of three components of displacement of reference stations, there is the possibility of M≥6.0 earthquakes in Xinjiang and its adjacent areas.

(3) It can be seen that the NS, EW and UP components of displacement of reference stations all present a high gradient zone of velocity change along the Wuqia-Bulungkol-Taxkorgan region, therefore, it is believed that there is a rapid change of stress in the region from Wuqia to Taxkorgan, which is a forceful place for stress accumulation and release.

(4) There are differential distributions of movement velocities for both EW and UP components in the middle section of the Tianshan Mountains, and such variation features make us pay attention to the judgment of earthquake risk in the middle section of the Tianshan Mountains, especially around Bayanbulak.

ACKNOLODGEMENT

The authors extend their great thanks to the National Earthquake Infrastructure Service, China Earthquake Administration, for providing data, and to research professor Wang Xiaoqiang from the Earthquake Agency of Xinjiang Uygur Autonomous Region for his guidance on data analysis.

This paper has been published in Chinese in the journal of Inland Earthquake, Volume 30, Number 1, 2016.

REFERENCES
 Gu Guohua, Zhang Jing. Time series of displacements from GPS observation at fiducial station in the crustal movement observation network of China[J]. Journal of Geodesy and Geodynamics, 2002, 22(2): 61–67 (in Chinese with English abstract). Gu Guohua. Recent progress in researches on crustal movements through GNSS (GPS) observations[J]. Recent Developments in World Seismology, 2007(7): 9–15 (in Chinese with English abstract). Jiang Zaisen, Ma Zongjin, Niu Anfu, Zhang Xiaoliang, Wang Shuangxu, Chen Bing. Approaches and preliminary results of crust movement researches based on the GPS in China[J]. Earth Science Frontiers, 2003, 10(1): 71–79 (in Chinese with English abstract). Li Guirong, Ailixiati·Yushan, Wang Xiaoqiang, Zhu Zhiguo, Paerhati·Zainula, Buaijieer·Kuerban. Study on abnormal features in EW displacement before Minxian-Zhangxian MS6.0 earthquake in Gansu[J]. China Earthquake Engineering Journal, 2013, 35(3): 542–548 (in Chinese with English abstract). Li Jie, Wang Xiaoqiang, Tan Kai, Liu Daiqing, Paerhati·Zainula, Jiang Jinxiang, Fang Wei. Analysis of movement characters of present-day active tectonics of northern Tianshan region[J]. Journal of Geodesy and Geodynamics, 2010, 30(6): 1–5 (in Chinese with English abstract). Li Jie, Liu Daiqin, Wang Qi, Wang Xiaoqiang, Zhu Zhiguo. Research on the characteristics of recent crustal movement of the South Tianshan Area in Xinjiang through the strain field and the change rate of baseline[J]. Journal of Seismological Research, 2012, 35(1): 59–65 (in Chinese with English abstract). Li Jie, Wang Xiaoqiang, Liu Daiqin, Paerhati·Zainula, Chen Shujiang. Distribution characteristic of strain rate field in Tianshan region by least squares collocation[J]. Acta Seismologica Sinica, 2015, 37(1): 103–112 (in Chinese with English abstract). Liu Guangming, Tang Yingzhe, Wu Fumei, Qin Xianping, Zeng Anmin. Coordinates and velocities of crustal movement observation network in China[J]. Journal of Geodesy and Geodynamics, 2012, 32(S1): 53–56 (in Chinese with English abstract). Shi Chuang, Zhao Qile, Lou Yidong, Li Min, Geng Jianghui, Ge Maorong, Liu Jingnan. PANDA: comprehensive processing software for satellite navigation systems and its research progress[J]. Spacecraft Engineering, 2009, 18(4): 64–70 (in Chinese with English abstract). Wang Xiaoqiang, Lu Xing, Liu Bin, Li Jie, Liu Daiqin. Study on present-day vertical crust movement and seismic activity in Xinjiang area[J]. Journal of Geodesy and Geodynamics, 2009, 29(6): 22–27, 31 (in Chinese with English abstract). Wang Xiaoqiang, Song Heping, Cheng Ruizhong, Zheng Liming, Li Jie, Li Guirong, Parhati·Zainula, Zhu Zhiguo, Fang Wei, Zhou Jun, Buerjer, Chen Shujiang, Sun Xiaoxu. A preliminary study on crustal movement of Tianshan and seismic activity based on crustal movement observation network of China[J]. Inland Earthquake, 2012, 26(2): 97–107 (in Chinese with English abstract). Wu Zemin, Bian Shaofeng, Xiang Caibing, Ji Bing, Jiang Dongfang. A GNSS integer ambiguity resolution method based on partial search strategy[J]. Journal of Naval University of Engineering, 2015, 27(1): 31–34 (in Chinese with English abstract).