Superimposed tectonic movement and karst erosion resulted in a combination of fractures and irregular caves in deep/ultradeep carbonate rocks, typically along major fault swarms. Outlining these fault-karst reservoirs mainly depends on recognizing the strong reflection in seismic profiles; however, characterizing their internal structures is still difficult, which are represented as weak amplitude in seismic profiles. This study intended to propose a method to dissect the internal structure of fault-karst reservoirs, which contains four steps: (1) elimination of the signal interference by the covering bed with strong energy and recognition of internal reservoirs with low energy based on seismic data conversion, frequency division, and inversion; (2) gradient structure tensor analysis based on an anisotropic Gaussian filter for fault-karst reservoir outlining; (3) internal faults and cave recognition on the basis of wave-based inversion; and (4) reevaluation of the number and scale of these faults and caves based on seismic recognition and well test results and verification of their volumes and hydrocarbon reserves. The method was used in the evaluation of the fault-karst reservoir in the Halahatang (HLHT) oilfield, which is located in the north of Tarim Basin. The results show that the fault-karst reservoirs along major faults and their internal structures are effectively recognized, and the error of the predicted depth of the reservoirs decreases from more than 20 m before to less than 4 m now; the drilling success ratio increases from 70% to 90%. In addition, the method supports the recognition of untapped fault-karst reservoirs around shut-in wells, which provides guidance for sidetracking plans. Further, by comparing the geophysical volume of fault-karst reservoirs and the reserve predicted by production performance, the untapped reserve in a certain reservoir can be evaluated; on this basis, producing wells received high yields by targeted acid fracturing. In summary, the method effectively improves the prediction accuracy and the recovery efficiency of fault-karst reservoirs.
Marine oil and gas fields account for one-third of giant hydrocarbon fields worldwide, whereas marine carbonate reservoirs play a leading role in oil and gas production. They are widespread in Tarim, Sichuan, and Erdos Basin and are the next major oil and gas exploration targets after continental clastic reservoirs. In the marine carbonate reservoirs in Paleozoic strata alone, the oil equivalent reserve reaches approximately 4.5 billion tons. Marine carbonate reservoirs are typically deeply buried (~7000 m-8000 m in depth) in China, especially the Ordovician carbonate reservoirs in Tarim Basin. Most of these reservoirs were untapped. In the reservoirs that have already been exploited, the recoverable reserve and the recovery efficiency are low in the first stage of exploration and exploitation. The cause of formation of these reservoirs was thought to be buried hills [1, 2], paleokarst cave systems, karst along unconformity , or interlayer karst  before 2015, whereas the reservoir patterns contained isolated caves, fault caves, and fracture pores. Vertical wells were prevalent in reservoir exploration, appraisal, and exploration, with main targets being “string beads” [5–7] in seismic profile. The drilling success rate in this stage was as low as approximately 70%.
Exploration and exploitation results in fracture-cave reservoirs these years show that oil and gas accumulation is closely related to deep faults. In most cases, these faults connect with oil sources and provide space for karst water flowing. Deep fault-fracture swarms and related karst caves and pores are therefore the main targets for exploration now . Xinbian et al. . named this complex system as “fault-karst reservoir” in 2015 based on abundant practices and interpreted its cause of formation as karst erosion along the fracture swarms by the subsurface water near deep fault zones. The whole system is therefore a typical fault-controlled karst system, and as the oil and gas migrate through these faults, they tend to accumulate at both the faults and karst zones under the seal of marlstone above or compact limestone around. The concept of “fault-karst reservoir” is a great breakthrough in geological theory of fracture-cave reservoirs, which provides new ideas for exploration and exploitation in marine carbonates. The concept is more and more accepted now, leading to several successful exploration [10–15]. In actual practice, the fault-karst reservoirs are divided into different types of trap for oil production. Geologists and engineers found that total yield varies greatly between them, among which the type of banded fault-karst reservoirs receives the highest yield . Based on the progress of related geological theory, production technologies such as targeted water injection and steam injection are applied to enhance oil recovery [16–18]. Nevertheless, further effective exploitation of these reservoirs is still changeling, since their great depth and strong heterogeneity make it difficult to depict their internal structures. An elaborate characterization of their internal structure and a reasonable and effective exploitation plan therefore become the major problem to be solved in the oil production.
Geophysical interpretation is the key to solve these problems in fault-karst reservoir quantitative characterization and evaluation. In previous studies, geologists and engineers failed to understand that the “string beads” feature in seismic profile represented a complex system of faults, fractures, and caves and pores [10, 19]. As the degree of exploitation goes higher, it is gradually realized that separate belts (no fractures or caves develop in these belts) exist in the fault-karst reservoirs and lead to complex structures inside the reservoirs [20, 21]. However, there is no detailed solution for the description of the internal structure characteristics of fault-karst reservoirs in the published literature. At present, the research mainly focuses on trap-type classification and evaluation, fault-karst reservoir edge characterization, physical simulation, and seismic response characteristic analysis. On the other hand, there are few reports on the characterization of the internal structure of fault-karst reservoirs. Vague characterization of internal structure of the reservoir is the major problem that constricts further deployment of high-yield wells and remaining oil extraction [8, 22, 23]. Difficulties encountered in a clear characterization of internal structure of the reservoir include (a) a clear outlining of caves with weak seismic reflection or near faults in seismic profiles, (b) accurate characterization of fault and cave systems (fault-karst reservoirs) near slip fault of small displacement , (c) reasonable explanation for the huge difference in production performance between wells in different fault swarms or even in one swarm, and (d) a feasible combination of production data and seismic interpretation to characterize the internal structure of fault-karst reservoirs. Technological breakthrough is imperative around quantitative characterization and elaborate evaluation of internal structure of reservoirs.
Seismic reflection feature is a basic tool for characterizing internal structures of fault-karst reservoirs. Based on seismic waves, related geological patterns, and production behavior, this paper discussed the feasibility of several technologies to outline and dissect fault-karst reservoirs, including covering bed signal elimination, boundary outlining using orthogonal tensor gradient, tensor gradient constrained FWI (full wave inversion) to recognize internal structure of reservoirs, and second iteration recalibration of the reservoir structures as well as reserve calculation based on FLCP (full lifecycle performance). These technologies provide new ideas for well infill and remaining oil recovery in carbonate reservoirs with complex oil-water relationship.
2. Information of the Study Area
2.1. Geological Background
HLHT (Halahatang) oilfield is located in the middle part of North Tarim Uplift of Tarim Basin  (Figure 1). The oilfield borders Luntai Bulge at north, North Sag at south, Yingmaili Low Bulge at west, and Lunnan Low Bulge at east, with a total area of 4000 km2. It is in the middle of south flank of Great North Tarim anticline, and the Ordovician carbonate strata constitute the west inclined wall layer of the Ordovician buried hills on Lunnan Low Bulge. The whole reservoir shows a giant nose structure and is further complicated by two sets of conjugated fracture striking NE and NW.
Middle Ordovician YJF Fm (Yijianfang Formation) and bottom YS Fm (Yingshan Formation) constitute the main reservoir of HLHT oilfield. Karst erosion happened along faults and related fracture swarms in the carbonate, which leads to reservoir development. Core sampling shows poor matrix porosity and permeability of carbonate rocks in this area, and secondary karst fractures and caves are therefore the main reservoir space. As recorded by well logs, the porosity of the two formations ranges from 1.8% to 90% with an average of 4.2%, and the permeability ranges from 0.01 to μm2 with an average of μm2. No apparent relationship exists between the porosity and permeability .
Crude oil in HLHT oilfield generally has low freezing point, density, and sulphur content; however, its density varies largely in different zones. From southwest to northeast, the oil gradually gets heavier, and the whole reservoir can be divided into 2 parts by oil density: light oil zone in the south and heavy oil zone in the northeast and east (Figure 2).
2.2. Production History
Since the beginning of exploitation at the north of Tarim River in 2009, the production in HLHT oilfield increased in seven consecutive years, then stabilized for a short time and fell in 2016. The composite water cut in Ordovician reservoirs reaches 36% currently (Figure 3).
3. New Methods and Workflows
3.1. Recognition of Low-Energy Fault-Karst Reservoirs
Strong reflection in seismic profiles has long been thought as response of effective fault-karst reservoirs; however, in actual practice, high-yield oil flows were obtained in some reservoirs with weak reflection [13, 27–30]. How to accurately outline and characterize the internal structural of these reservoirs was still challenging. The Ordovician reservoirs in HLHT are covered by thick clastic layers, and the interface between them forms strong continuous reflection in seismic profiles. These strong signals shield (weaken) the seismic response of the reservoir below. When the fractures and pores of the reservoirs are smaller, the shield effect is stronger. Similar situation also happens in the central Tarim Basin, the Tahe oilfield, and the Permian Maokou Formation in Sichuan Basin. These fault-karst reservoirs are thus difficult to be recognized in seismic profiles, unless the response of the interface is eliminated.
The process of the elimination contains three steps. First, seismic data volume is transited from time domain to wheeler domain, and the new data volume is named as Volume A. In this step, responses of the tectonic information in the seismic volume and other attribute volumes are removed, and the syndepositional distribution of the formation is restored. The occurrence and phase of the seismic events in Volume A therefore directly reflect the stratigraphy and sedimentary information. Particularly, the continuous strong reflection above the fault-karst reservoirs is flattened, and the seismic events have better phase congruency now. Second, wave packets with good phase congruency are extracted by PCA (principal component analysis), which constitute Volume B. Volume B is then subtracted from A to get volume C, which represents the weak reflections in the volume in wheeler domain. C is further inverse transited to time domain. Third, through frequency division and impedance inversion of Volume C in time domain, characteristic of reservoirs with weak reflection is finally featured.
The workflow above is a seismic filtering process, in which the low-frequency component is removed, and the response of small faults is featured. More importantly, strong reflection of the covering clastic layers is also eliminated. Then, the original weak reflection of fault-karst reservoirs is enhanced, and the fault swarms become clearer in seismic profile (Figure 4).
3.2. Fault-Karst Reservoir Boundary Outlining by Gradient Structure Tensor
GST (gradient structure tensor) analysis is a new seismic attribute analysis method, it brings the concepts of image texture processing into seismic interpretation. Seismic data is processed as images in GST analysis. Textures and structures in these seismic images are further recognized by combined parameters  and the scale of tensor eigenvalue (Figure 5). These textures (laminated or chaotic) and structures actually represent anomalies in rocks (e.g., faults and caves) and can be used as signals for internal structure interpretation and automatic boundary outlining  for fault-karst reservoirs.
This study proposed a GST analysis method using an anisotropic Gaussian filter. The nature of the method is texture recognition and analysis in seismic images by scaling the amplitude gradient. The advantage of this method is that firstly, it recognizes reservoir boundaries of various scales and generates chaotic volume that, compared to normal structure tensor, reflects sedimentary textures with stronger noise immunity. Secondly, it generates planar structure gradient that features horizontal discontinuity, which is favorable for boundary outlining of complex deposition. Thirdly, it supports more accurate inclination calculation and, combined with III-generation coherence analysis method, generates high-precision coherence and curvature volume.
is a diagonal matrix, whose eigenvalue is .
Seismic response of fault-karst reservoirs is distinct from that of neighboring rocks, and the in three directions are therefore different from each other. The tensor gradient property is used to describe the boundary of fracture solution. The larger value is, the stronger the chaos is in the three directions, and the more likely it is to reflect the boundary of the fracture solution. The chaotic values represent the boundary difference between the fault solution and the surrounding rock. Their sum is defined as the result of GST. Higher value means stronger heterogeneity, which is more likely to be the boundary of fault-karst reservoirs, as shown in Figure 6.
3.3. Internal Structure Characterization of Fault-Karst Reservoirs
Wave impedance inversion is a commonly used technology in (semi-) quantitative characterization of reservoirs. However, it is capable only of recognizing laterally developed cave reservoirs. For the vertically developed fault-karst reservoirs, this technology is usually ineffective [5, 33]. In this study, an improved waveform inversion method is proposed based on strata-signal elimination (see Section 3.1 for details).
Strong reflection of the covering clastic beds is already eliminated in the workflows above. On this basis, GST is crossplotted with wireline-calculated impedance to determine the mathematic relationship between them. The function is then used to transform GST volume to wave impedance volume, which is further used in inversion as original volume and constraint. This method generates data volume that is capable of both defining the boundary and characterizing the internal structure of the fault-karst reservoirs.
Based on removing the strong reflection at the interface of the overburden strata of the fault solution, the intersection analysis of the gradient structure tensor property and the well side impedance of the known well is carried out, and the relation formula of converting the attribute body to wave impedance body is obtained. The gradient structure tensor is transformed into wave array antibody, which is used as the initial model of wave impedance inversion. Then, the initial model is applied to the waveform indication inversion process as a constraint condition, and the inversion data which can characterize both the contour characteristics of the fault solution and the internal structure characteristics are obtained.
The results show that faults directly control the development of karst caves, and the reservoirs extend largely in a vertical direction. Horizontally, the reservoirs can be divided into fault zone, break zone, and compact zone (Figure 7).
3.4. Second Iteration Recalibration for Reserve Evaluation
Two methods are generalized for reserve calculation and evaluation in the previous studies. One is the reserve calculation by geological analysis, whose result is named as “geological reserve.” Geological models are used in this method to count the total volume of the caves, fractures, and pores, and the reserve is calculated on this basis. The other method is the reserve evaluation by analyzing the production behavior, whose result is called the “behavior reserve.” With this method, the reserve is evaluated using actual production data and performance. This is an effective way to evaluate actual reserve which the producing wells control. Great difference between the results of the two method is another problem commonly encountered in fault-karst reservoir characterization. Commonly, geological reserve is much larger than behavior reserve, which leads to wrong decisions in development mode planning and remaining oil exploring. This study proposed a second iteration recalibration method, which combines production behavior, testing results, and production decline laws. The internal structures and volumes of the reservoirs interpreted from seismic profiles are adjusted in this method, which lays foundations for more precise development of fault-karst reservoirs.
Second iteration recalibration contains 3 steps: (a) internal structure characterization using well test results, (b) behavior reserve calculation, and (c) seismic interpretation calibration by the results of (b). The relationship between geological reserve and behavior reserve is established from the above steps, which can be used in reserve evaluation for undrilled fault-karst reservoirs.
3.4.1. Internal Structure Characterization Using Well Test Results
Geologists refer to well testing as the “eyes for the reservoirs.” Well test provides the geological and production information of the well-controlled area in the reservoirs. Interpretation model for the results of well test in fault-karst reservoirs is more complicated [34–37] when compared to that in conventional reservoirs. Commonly, the log-log diagnostic testing curves are used for the recognition of the number of reservoir bodies, calculation of the equivalent porosity and permeability, and speculation of the production behavior of fault-karst reservoirs (Figure 8).
3.4.2. Calculation of Behavior Reserve
The volume of fault-karst reservoirs can be evaluated by GST and other geological technologies such as seismic interpretation, but the reserve calculation and recalibration also need consideration of production behavior. Behavior reserve (reserve that is speculated or reflected by production performance) calculation methods mainly includes material balance method and production transit testing [38–46]. Strong heterogeneity, various well behaviors, and small aquifer pressure are distinct characteristics of fault-karst reservoirs. In this study, a comprehensive method that combines Blasingame, AG (Agarwal-Gardner), NPI (Normalized Pressure Integral), and FMB (flowing material balance) [45, 47–52] was used to calculate behavior reserve of the reservoirs. The advantage of this comprehensive method is its ability to calculate reserves by production behavior monitoring without shutting in the producing wells.
Take FMB as an example. In most cases, fault-karst reservoirs are unsaturated (by natural gas) and are exploited above the bubble point pressure. The original pore volume remains the same before and after crude oil exploitation. The reserve of crude oil, however, changes with exploitation and pressure variation, as in shown Figure 9. If the original reserve is referred as , the accumulative production is , the remaining reserve is , the original and present formation pressure is and , and the compressibility of the rock is “” then in a certain time of exploitation.
The above equation clearly illustrates the relationship between initial pressure, flowing pressure, and rate.
As in Figure 10, it is apparent that the pseudo-steady-state equation is the sum of two distinct pressure loss components:
A graph of normalized pressure vs. time produces a straight line (Figure 11), the slope of which provides the original oil in place. This is the most basic form of flowing material balance and is directly applicable to undersaturated oil reservoirs. In well testing, this procedure is commonly referred to as a Reservoir Limit Test.
Plotting normalized rate vs. normalized cumulative (Figure 12), the following plot is derived, where the -intercept is equal to the original oil in place.
The above equation and associated figure describe the flowing oil material balance analysis for a volumetric circular reservoir.
Standard chart in Harmony® is used for the calibration of behavior reserve. The reservoir volume determined by well testing (see Section 3.4.3) is also used to give the final calibrated reserve of the fault-karst reservoir. Figure 13 shows the charts that are used in the calculation of the behavior reserve in well P5, which reaches m3.
3.4.3. Workflow of Second Iteration Recalibration for Reservoir Volume
There are several factors that lead to an “amplifying effect” of the geological features in seismic profiles: (1) seismic wave attributes, such as horizontal-vertical proportion, resolution, and SNR (signal to noise ration); (2) interpretation methods, such as the migration method; and (3) rock characteristics, such as dimension, depth, and lithology (especially the lithology of neighboring rocks). These characteristics are different in various districts, which results in different extent of amplification. Exclusive calibration curve, mainly established by seismic forward modeling, is therefore necessary for each study area. Important elements for establishing these curves include (a) simplified geological model that reflects the geological settings of the area; (b) acquisition parameters for forward modeling should be the same as the acquisition parameters for the seismic data; and (c) the processing methods and parameters for forward modeling results should be the same as those for seismic data.
Seismic forward modeling is an effective way to establish the recalibration curve. Through seismic forward modeling of fault-karst reservoirs with different volumes, the “pseudoseismic responses” of the reservoirs can be obtained. The attributes of these responses such as volume or amplitude were further plotted vs. reservoir volumes calculated from geological models, and thus, the relationship between them is established. This relationship can be used to recalibrate the actual seismic response and correct its attributes. For example, in actual practice, it is concluded that small caves would be amplified horizontally, and seismic forward modeling tells us to what extent the width of caves was amplified. The relationship (Figure 14) therefore helped us to correct the “pseudoreservoir volume” to “actual volume,” with errors small enough to accurately calculate the reserve. As shown in Figure 14, where the seismic response width is beyond 250 m, the correction coefficient approaches gradually to 1: apparently, when the reservoir is large enough, the amplification effect is negligible. However, where the width of response is smaller than 250 m, the correction coefficient increases as the width decreases. The way of correction from “pseudoreservoir volume” to “actual volume” is to divide the width of the beads in seismic profiles by corresponding correction coefficient.
Second iteration recalibration of seismic interpretation and reserve calculation is a critical technology in oilfield evaluation and exploitation. It is directly used in reserve calculation, potential evaluation, and remaining oil development and contributes to higher yield. Based on the chapters above, the workflow of this innovative technology is summarized as follows:
Calibrate the seismic-interpreted (by GST) reservoir internal structures by well test results, and accurately characterize the number of caves and their interconnectivity in the reservoir
Calculate the bulk volume of caves and fractures of the reservoir, i.e., the reservoir space volume, by the recalibrated seismic data (in step a), and determine a more precise geological reserve
Recalibrate the geological reserve (in step b) by behavior reserve, and establish the relationship between the two reserves. Plot the curves of the relationship in a certain fault swarm, which is further used in evaluating the reserve in the whole area
4. Application and Discussion
Using the recognition of low-energy fault-karst reservoir technology, the problem of seismic prediction of weak energy formation inside fault solution in P6 well area of HLHT oilfield is effectively solved, the prediction accuracy of reservoir is improved, and the drilling success rate is significantly improved. Based on the fine identification technology of the internal structure of the fault-karst reservoirs, the reservoir around the abandoned well was reevaluated, and the high oil production was obtained after sidetracking well P10. By using the second iteration recalibration for reserve evaluation technology, the reservoir scale near the low-producing well can be reestablished, which provides geological basis for acid fracturing and other stimulation measures and significantly improves the effect of stimulation measures.
4.1. Low-Energy Fault-Karst Reservoir Prediction in P6
Using the innovative and comprehensive technology packet mentioned above, this study characterized the internal structure of the Ordovician fault-karst reservoirs in the P6 block of HLHT oilfield. The result is shown in Figures 15 and 16, which provides foundation for fault segment characterization and explains the discrepancy in the scale of reserves between different fault swarms and different structures in the same reservoir. The result of seismic interpretation of internal structures is verified by layer thickness comparison: first, multiply the time duration of the reservoirs in each layer by average seismic velocity to get the thickness of each reservoir, and add them up to get the total thickness () of the reservoirs; is then compared with the reservoir thickness drilled (). The error () is less than 4 m, as shown in Table 1.
The accurate characterization and evaluation of the fault-karst reservoirs effectively guide the exploration and well employment. The drilling success rate of the wells drilled in recent years increased from 70% to 90% now.
4.2. The Internal Structure Characterization of P10
The location of well P10 is shown in Figure 17. Well log recorded poor porosity and permeability of the reservoirs. After acid fracturing, the well test showed that the formation drilled was dry layer, which indicated drilling failure.
Recent interpretation of the reservoir by the new methods found that the borehole went through the mud filling of the fault-karst, and acid fracturing was not so effective to channel the neighboring reservoirs. A more precise characterization of the internal structure of the reservoir indicated the existence of internal unfilled caves, which was proved by later side tracking. The side-tracking well received good testing result and high yield later, as shown in Figure 18.
4.3. Volume and Reserve Recalibration of P11
Main target formation in the north Tarim Basin is the deposits of the middle Ordovician carbonate, which is covered by clastic deposits. Seismic wave velocity is quite different in the two formations (average velocity in compact carbonate is 6000 m/s, and that in clastic rocks is just 4950 m/s), which leads to great difference in their wave impedance. Consequently, the interface of the two formations forms continuous strong response in seismic profile, which is recorded as a peak with anomalous large amplitude. Previous exploration and exploitation mainly focused on the middle Ordovician reef flat reservoirs in the carbonate platform margin, instead of the fault-karst reservoirs in the platform center. The latter ones mainly develop fractures and corrosion cavities, whose impedance is similar to that of the compact carbonates and form weak reflection in seismic profiles. Therefore, the strong response of the interface shields the response of the fracture-cave carbonate reservoirs, resulting in difficulties in predicting these weak-reflection reservoirs.
A case study of well P11 shows the effectiveness of strong response elimination: in the original seismic profile (Figure 19(a)), two fracture-karst reservoirs (recorded as ① and ②) were recognized, and in the enhanced profile (Figure 19(b), with strong response removed and weak reflection featured), a third reservoir (recorded as ③) was found. To further verify the accuracy of the characterization, the cumulative production vs. pressure drops of well P11 was plotted (Figure 20). The three segments in the diagram represent three reservoirs, which is identical to the result of seismic interpretation.
On the basis of strong reflection elimination, this study recognized 10 wells (out of 12 wells in P11 block) that drilled through caves with weak reflection, and the rate of cave drilling reached 83%. These caves are small corrosion cavities with a height less than 1.5 m. In the drilling process, circulation lost happened rather than drilling break. Application of this technology in the whole P11 block increased t new oil reserve, which strongly supports the next round of well employment.
The method has a good application prospect in the fault-karst reservoir recognition in ultradeep carbonate formation and their outline depiction, internal structures characterization, and reserve evaluation.
On the basis of seismic data conversion, restructuring, frequency division, and inversion, the internal structures of the reservoirs with low energy were effectively characterized; based on the technique of gradient structure tensor analysis, the boundaries of the fault-karst reservoirs were depicted, which improves the prediction accuracy of such kind of reservoirs in ultradeep carbonate formation. With the application of this method in HLHT oilfield, the error of the predicted depth of the target formation decreased from more than 20 m to less than 4 m, and drilling success ratio in the past three years increased from 70% to more than 90%
By combining the interference signal elimination and wave-based inversion, the internal structures of the fault-karst reservoirs and untapped reserves were accurately characterized and evaluated. The application of this technique in well P10 in HLHT supported an effective sidetracking plan towards the untapped reservoir, and high yield was received later
On the basis of well test interpretation, the number and scale of the fault-karst reservoirs were compared with reserves predicted by production performance, and thus, the untapped reserves can be evaluated in a certain reservoir. With the application of this method in well P11 in HLHT, three separate caves with different depths were recognized in a fault-karst reservoir, and P11 was back on production later by targeted acid fracturing
In summary, the internal structure characterization of fault-karst reservoir is an integrated method for fault-karst reservoir evaluation and shows its effectiveness in the application of HLHT oilfield. This method can be effectively used in reservoir prediction, trajectory optimization, shut-in well opening, and EOR.
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare no conflict of interest.
Conceptualization was done by Peng Cao and Shaoying Chang; methodology was done by Peng Cao; validation was done by Peng Cao and Yongjin Zhu; writing the original draft preparation was done by Peng Cao and Shaoying Chang; writing the review and editing were done by Peng Cao, Shaoying Chang, Yongjin Zhu, Jinlong Shen, Zhanfeng Qiao, Guanming Shao, and Xiaowei Sun.
This research was funded by the National Science and Technology Major Project (grant number 2017ZX05008005).