#### Abstract

Quantifying the condition of large cultural relics, such as marine archaeological shipwreck, is important to verify stability and reliability. Deformation monitoring system plays a key role in the preservation and long-term conservation of cultural relics. Two-dimensional digital image correlation (2D-DIC) method has proven its efficiency in being able to provide accurate quantitative information of structural deformations. In this study, a deformation monitoring system with four cameras based on 2D-DIC is developed to perform noncontact, optically based measurement to monitor the deformation of shipwreck in museum environment with low and varying illumination. Because the consistency of the accuracy of 2D-DIC measurements for different locations is the most basic requirement in the application of structural deformation monitoring, selecting the appropriate exposure time and quantifying the bias errors on 2D-DIC measurements should be helpful to the optimal use of this optical nondestructive testing technique. A theoretical criterion is deduced to quantitatively characterize the dependence of interpolation bias upon natural patterns and illuminations. Then, an exposure adjustment scheme is built based on the aforementioned criterion. Numerical experiments reveal that the exposure adjustment scheme is able to provide consistency interpolation error for different natural patterns even though the environmental illumination is different as well. The deformation monitoring system with the proposed exposure adjustment scheme is promising for developing flexible and robust in situ structural health monitoring for use in museum environment with low and varying illumination, making 2D-DIC technique a really useful tool for in situ long-term monitoring of large cultural relics.

#### 1. Introduction

Structurally intact wooden shipwrecks on the seabed constitute time capsules carrying significant historic, archaeological, social, and scientific information [1]. In 1974, after more than 700 years buried in marine clay sediments off the Quanzhou Bay (24°37′N, 118°37′E), the Song Dynasty shipwreck was raised from the seabed with the aim of restoration, preservation, and long-term conservation and has since been kept at the Quanzhou Maritime Museum. The Song Dynasty shipwreck (Figure 1) provides a unique environment for research into microbial wood deterioration, strength and stiffness degradation, biogeochemical cycling of iron and sulfur, and the medium to long-term impacts of conservation of marine archaeological wood. From a conservation or structural perspective waterlogged archaeological wood differs from modern wood because it has been degraded by various irreversible changes of physical, chemical, and biological processes in the marine environment [2]. Environmental degradation [3] of marine archaeological wood results in loss of cellulose [4, 5], which provides longitudinal strength and stiffness. These changes result in a very porous material, due to the formation of cavities and cracks, and the result is a significant loss in mechanical strength [6]. The major challenges faced during conservation of marine archaeological wood include the nonuniform degradation of timbers and the seasonal variation of relative humidity (RH) in museum. The marine archaeological wood substance strives for equilibrium with the ambient relative humidity. In the environment of high relative humidity, the cavities within the wood will fill with water, which if descended can result in severe shrinkage of the wood due to an inability to resist the strong capillary forces during the water-desorption process. Uncontrolled water-desorption process can lead to collapse, shrinkage, distorted shape and surfaces (twists, warping, cracks, and splits), disintegration, precipitation of salts, and corrosion products in the marine archaeological wood [7]. Therefore, it must be achieved while avoiding the destructive forces of water that can result in shrinking, warping, and cracking of the wood in the current preservation environment. In addition, the process such as shrinking and warping could result in further local deformation and subsidence of the wooden shipwreck. An easy-to-use, effective deformation monitoring system for shipwreck condition assessment, would allow the detection of potentially dangerous situations at an early stage and providing an efficient basis to control preservation environment or reinforce the shipwreck, while not interfering with exhibition.

Digital image correlation (DIC) among other optical nondestructive testing (NDT) methods is becoming more widely accepted for materials and structural monitoring. The two-dimensional DIC is applied to a series of images of a deforming planar object acquired by a single camera [8]. The surface of the object is speckled to present a randomly distributed intensity pattern. The principle of the 2D-DIC method lies in determining the local correlation of two speckle images. The randomness of the speckle pattern is essential for obtaining a unique solution in the correlation process. Reference subsets in the reference image are compared to subsets in a deformed image to find the target subset, which is the subset in the deformed image that shows the maximum similarity with the reference subset.

To achieve high-accuracy measurements, various aspects involved in the implementation of 2D-DIC must be given careful consideration. As an optical NDT method based on image matching, the measurement of 2D-DIC relies heavily on the quality of the acquired images. Various researchers have demonstrated that the speckle pattern has an important influence on DIC measurements. A number of speckle image assessment criteria have been presented in recent literature [9–18]. Among them, the mean intensity gradient (MIG) [10] and the sinusoidal approximate formula for noise-induced bias [17] are typical representatives. In addition, efficient speckle pattern optimized algorithms and corresponding bias error estimation and correction approaches have been proposed [19–23]. In practical application of DIC measuring and/or monitoring, a continuous and even illumination on the test object with stable and controllable white light source is an important means to ensure sufficient and constant contrast for reliable and accurate speckle image matching and to produce a series of quality speckle images for subsequent DIC analysis. However, in an actual museum environment, especially in situ long-term monitoring of cultural relics in the exhibition hall where medium and/or high light-sensitivity objects are on display, the requirement of continuous and even illumination is unable to be achieved, as ambient light adopts general lighting and has requirements for total amount of illumination (annual lighting exposure). In 2009, China issued a National Standard for the Design of Museum Lighting (GB/T 23863-2009). This code articulates the different deterioration phenomena caused by light and establishes the highest levels of illumination to limit damage [24]. China has divided museum objects into three categories, low, medium, and high sensitivity. For Song Dynasty shipwreck that belongs to high sensitivity objects, the total amount of illumination should be less than 360 000 lux-hours per year, which limits the number of display days to less than 300 per year when the illuminance is 150 lux for eight hours.

On-site investigation and measurement result show that the environmental illumination of Quanzhou Maritime Museum below 100 lux with adoption of natural lighting. On the other hand, the environmental illuminations fluctuate from the medium illumination of daytime to low illumination of dark night and vary from location to location on shipwreck. Speckle patterns markedly affect the accuracy of the correlation results because they are associated with the similarity before and after deformation. When the ambient light illumination changes significantly, the contrast of speckle patterns may change accordingly. In practice, the most detrimental lighting conditions are the ones in which the illumination gradient change over the image. As a result, the similarity between the deformed images and the originally recorded reference image may decrease notably, resulting in a failure of 2D-DIC analysis. In order to suppress influences of varying ambient light on speckle patterns and acquire quality speckle images suitable for 2D-DIC analysis, Pan et al. [25] proposed a monochromatic light illuminated active imaging DIC method for high fidelity deformation measurement. Simončič et al. [26] modified the inverse compositional Gauss-Newton algorithm to achieve high-accuracy measurement in image sequences with varying lighting intensity by DIC. Xu et al. [27] proposed a weighted normalized gradient-based algorithm to decrease the influence of illumination and got better results by DIC in case of nonlinear gray-value-based illumination variation.

The speckle pattern serving as the information carrier is a key issue in relation to the accuracy in using DIC. The speckle patterns cannot be fabricated artificially for protection reasons, so natural wood texture distribution inherent in shipwreck surface was used as natural pattern to perform correlation calculation. Monitoring of multiple locations is necessary for the deformation monitoring system. It should be noted that the natural patterns on different locations of the shipwreck may lead to distinctly different bias errors. The bias error of 2D-DIC includes the random error, interpolation bias caused by interpolation, and the noise-induced bias caused by coupling of interpolation and noise. Among them, interpolation error is the main bias error form. In order to synthetically evaluate the result of deformation monitoring, the same bias errors are required for different monitoring regions. This requires that the deformation monitoring system based on 2D-DIC must have the capability to estimate bias errors for natural patterns and function to adjust the quality of natural patterns in order to meet the consistency requirements of bias errors. In fact, to ensure constant bias errors of different natural patterns captured by ordinary industrial cameras, an alternative approach is to adjust exposure time. The industrial camera's exposure time determine the amount of light that reaches the image sensor. Changing exposure time can change the brightness and contrast of speckle patterns. A good speckle pattern on the test sample surface should have high contrast. The decrease in illuminance causes the decrease in the grayscale gradient of the speckle pattern, and then it will lead to lower measurement accuracy significantly [28]. In order to acquire effective speckle patterns with the constant bias error suitable for DIC monitoring in low and varying illumination, how to choose appropriate exposure time is a problem to solve.

In this work, a deformation monitoring system based on 2D-DIC was established to investigate the deformation of the Song Dynasty shipwreck. A quantitative criterion for bias error estimation of natural pattern is proposed based on theoretical analysis and experimental verification. Next, we apply the quantitative criterion as the basis for exposure adjustments. The remainder of present paper is organized as follows. Section 2 briefly explains the 2D-DIC measurement principle and the influence of different natural patterns with different illumination conditions on measuring accuracy. In Section 3, a quantitative criterion based on the theoretical analysis of interpolation bias is presented, and adjustment scheme of exposure is experimented and analyzed. In Section 4, the influence of image noise is discussed. Finally, conclusions are summarized in Section 5.

#### 2. The Influence of Speckle Pattern on 2D-DIC Measurements

##### 2.1. Principle of 2D-DIC

2D-DIC uses image registration algorithms to retrieve full-field relative displacements between undeformed (reference) image and deformed images [29]. Each pixel of these images stores a grayscale value matrix because of a pattern at the surface of the specimen called speckle pattern [30]. DIC relies on the speckle pattern to obtain surface displacement fields [10]. In subset-based 2D-DIC algorithms, the reference image is divided into smaller regions referred to as subsets. Reference subsets are initially a contiguous rectangular group of points that are on integer pixel locations. The deformation is assumed to be homogeneous inside each subset, and the deformed subsets are then tracked in the deformed image using the correlation matching algorithm. To find the deformation of a subset, 2D-DIC algorithms find the extremum of a correlation function. The practical and robust zero-mean normalized sum of squared difference (ZNSSD) correlation criterion is used to evaluate the similarity between the reference and deformed subsets [31, 32]. The grayscale intensity functions of the speckle pattern in each subset is represented by a function pair which refers to the reference configuration and to the deformed incremented state. Here are the coordinates of a reference subset point and the coordinates of a deformed subset point. The correlation function is the subset DIC residual, which is the zero-mean normalized sum of squared difference between the gray levels in subsets extracted in the deformed and reference images. For each subset, we have

Functions and correspond to the mean grayscale values of the reference and current subset with denoting the total number of points within subset .

The transformation of the coordinates from the reference subset to the deformed configuration is constrained to a linear, first-order transformation referred to as subset shape function , which is utilized to depict the shape of the deformed subset relative to the reference subset, is the local coordinates of the pixel point in each subset, and is the displacement parameter vector. It can be written in the following form:Here with denoting the displacement components, and the other four parameters represent the displacement gradients. If the linear intensity variation is assumed, the ZNSSD correlation criterion together with an affine shape function and an incremental shape function is defined by the following expression: where and denote the grayscale levels at of reference image and the deformed image.

To accurately match the two subsets, the roles of the deformed subset and reference subset are reversed in inverse compositional Gauss-Newton (IC-GN) algorithm [32–34]. A backward matching strategy is used by exerting an incremental shape function to the reference subset and comparing it with the deformed subset, where the affine shape function is defined. The incremental shape function can be defined in the same manner as the affine shape function .where denotes the incremental deformation parameter vector, which is calculated by the minimization procedure. In order to solve for , a first-order Taylor expansion of (3) with respect to is performed. We obtainwhere is the gradient in the x- and y-directions of the reference subset, which can be calculated by interpolation. In this work, quintic B-spline interpolation is employed. In (5), the term represents the Jacobian of the shape function, which can be written as Minimization of with respect to gives the least-squares solution of :where and .

Based on the minimization procedure, the incremental shape function is inverted and composed with the current estimated shape function. The obtained result is the updated shape function of the deformed subset, from which the incremental parameter vector can be determined. The operation can be expressed as follows:

The iteration is repeated until the preset convergence condition is satisfied. In this work, the convergence condition, , is predefined.

##### 2.2. The Influence of Speckle Patterns on 2D-DIC Measurements

Speckle patterns have an important influence on the accuracy in the displacements measurement of DIC [23]. To demonstrate the influence of speckle patterns on 2D-DIC measurements, eight reference speckle images obtained from real experiments were utilized in the following numerical experiments. The image acquisition operation is done in a darkroom with an illuminance of 0. In this environment, the ambient illuminance can be adjusted by using a brightness-adjustable white artificial light, and the reading of the illuminance value is read by a digital illuminometer. We slowly increased the brightness of the artificial light in the condition that the room is completely dark, continuously read the ambient illuminance values by digital illuminometer under various brightness conditions, and captured image until the illuminometer reading is approximately equal to 60 lux or other values.

Figure 2 illustrates the reference images and their corresponding power spectrums. The pattern A was acquired by randomly spraying black and white paints on flat specimen surface under the condition of the illumination of 60 lux and exposure time of 100 ms. The patterns B, C, and D were obtained from the same specimen of the speckle pattern A but different illumination or exposure time. The pattern E was acquired another flat specimen surface under the same condition of the pattern A. The patterns F, G, and H were obtained from the same specimen of the pattern E but unequal illumination or exposure time. All of the speckle patterns were acquired by a CCD camera with a fixed aperture lens (f = 5.6). The specific parameters for different speckle patterns are listed in Table 1. The illuminance is adjusted by changing the brightness of the white artificial light. The exposure time is adjusted by the camera software.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

For each speckle pattern, 20 subpixel translated images are generated by applying appropriate shifts in the Fourier domain of the reference speckle pattern according to the shift theorem [35]. The subpixel displacements are applied in the x-direction, ranging from 0 to 1 pixel with a spacing of 0.05 pixels between two successive speckle patterns. The displacements of each translated speckle pattern were computed by using IC-GN subpixel registration algorithm.

To quantitatively evaluate the accuracy of different natural patterns adopted in this numerical experiment, the systematic error of the computed displacements associated with speckle pattern is represented as mean bias error [36]. The mean bias error of the measured displacement is defined as follows: where represents the mean of the estimated displacements and denotes the exact imposed displacement.

Figure 3 shows the mean bias error as a function of preassigned subpixel displacements for the eight speckle patterns (solid lines). The sinusoidal-shaped mean bias error can be attributed to the interpolation error. This observation agrees well with existing results reported in [37]. By inspection, it is seen that speckle pattern has a significant influence on the bias error; the mean bias errors of different speckle patterns with the same illumination and exposure are different (compare the experimental results by DIC of speckle patterns A and E); lighting condition has a significant influence on the bias error; the mean bias errors of different speckle patterns with different illumination and the same exposure are different (compare the experimental results by DIC of speckle patterns A and G); for given speckle patterns, whether or not the patterns and illuminations differ, it is always possible to find appropriate exposure time to produce approximate bias errors.

#### 3. Deformation Monitoring System

##### 3.1. Systematic Error Estimation

Sinusoidal-shaped systematic error in 2D-DIC was periodic with a period of 1 pixel. Schreier et al. attributed this periodic error to imperfect interpolation [35]. To obtain subpixel accuracies, grayscale must be evaluated at noninteger locations in 2D-DIC. Therefore, gray values and gray-value derivatives must be interpolated. On the other hand, the IC-GN algorithm must estimate the gradients of the reference image. Gradients can be estimated by interpolation. Nonideal interpolation will lead to this systematic error, which is also called interpolation bias [38]. The interpolation bias is inevitable because the ideal interpolation cannot be implemented in practice. If the dependence of interpolation bias upon the interpolation algorithm and the speckle pattern can be determined analytically, the quantitative criterion for interpolation bias estimation can be predicted, assessment of speckle patterns will be benefited, and appropriate exposure time for different patterns under different illumination will be gained based on this quantitative criterion.

In the absence of noise, the measured displacement is not equal to the actual displacement due to subpixel interpolation and gradient estimation. The imperfect interpolation will introduce systematic errors, which are referred to as interpolation bias [35]. Reference function with domain is sampled to produce a reference sequence . For two-dimensional situation, translate towards the positive direction of by units and towards the positive direction of by units to obtain the deformed function , which is sampled to produce a deformed sequence . If convolution-based interpolation is employed, the reconstruction of the reference function and deformed function are where and denote the interpolation basis for the reference sequence and deformed sequence, respectively. For convenient description, suppose a continuous function with continuous-time Fourier transform and its translated version . The functions and are, respectively, sampled to produce the reference sequence and the deformed sequence ; (10) and (11) can been rewritten as

The corresponding continuous-time Fourier transform is

For translational situation, the following assumptions are used:

According to (7), the measured displacement satisfies where is the interpolation error sequence, is the gradients of reference sequence, and their discrete-time Fourier transforms are and , respectively. and are the discrete-time Fourier transforms of and , respectively. Applying Parseval’s theorem, the frequency representation of (18) can be written as The Poisson summation formulae are employed to present in terms of and in terms of :Therefore, the form of is given byThe continuous Fourier transform of is denoted as provided that . An application of the Poisson summation equation yields in terms of :

With the expressions of and , (19) can be rewritten as

Since in practical situations , it is reasonable to employ a linear approximation

Substituting (25) into (24), utilizing the band-limited hypothesis of , and recognizing that and are approximation of an ideal low-pass filter, for simplicity, let actual displacement along the -axis be ; it is evident that there is no interpolation bias along the -axis due to symmetry. Hence, can be approximatively derived as follows:

Equation (26) becomes identical to (41) given by Su et al. in [39]. Similarly, (26) explains the well-known sinusoidal-shaped curves of the interpolation bias with respect to the prescribed subpixel displacement; it also explicitly presents the dependence of the interpolation bias upon the reference function and interpolation basis and in the frequency domain. Within (26), is the power spectrum of is exclusively determined by the subpixel and gradient interpolation algorithm, and it has been called the interpolation bias kernel by Su [38, 39]. In this work both the deformed image and the gradient of reference image were interpolated quintic B-spline interpolation (known as a convolution-based interpolation). The uniform quintic B-spline basis function can be expressed aswhere .

Equation (26) implies that the interpolation bias is mainly determined by the interpolation bias kernel and the power spectrum. The amplitude of interpolation bias plays a central role in this work and is of importance because it quantifies the bias error response for different patterns, thus providing a criterion of exposure adjustments of the deformation monitoring system.

The interpolation biases of the eight speckle patterns utilized in Section 2.2 were estimated by theory (see (26)). The interpolation biases by DIC and by theory are illustrated in Figure 3. Figure 3 demonstrates that the interpolation biases by DIC and by theory (dashed lines) show good agreement. The correctness of theoretical estimations of interpolation bias is verified. Figure 4 indicates the standard deviation errors of computed displacements for these eight speckle patterns. It is observed from Figure 4 that the standard deviation errors of DIC are approximately stable and do not depend on the imposed subpixel displacement.

##### 3.2. Exposure Adjustment Scheme and Experimental Verification

In order to acquire appropriate speckle patterns with constant systematic interpolation basis for multiple monitoring locations, the exposure time of cameras must be adjustable. Based on this consideration, a deformation monitoring system with multicamera based on 2D-DIC is established in this work. The deformation monitoring system consists of four illuminometers, four digital CMOS cameras, and four prime lenses. To minimize the negative impact of natural patterns and varying lighting condition and ensure the same accuracy at different location of the shipwreck, we develop an exposure adjustment scheme, which can adjust exposure time of the digital CMOS cameras according to the interpolation bias amplitude of patterns captured at different locations. For this purpose, an initial value of exposure time and interpolation bias amplitude should first be selected, and then the illuminance of regions of interest of the shipwreck is monitored in real time by illuminometers. Speckle images of the surface of the shipwreck are acquired by CMOS cameras. Afterwards, the amplitudes of interpolation bias of each image are firstly calculated and then compared with the initial amplitude ; if approximately equals , it means that the interpolation bias estimations of this pattern are identical with the predefined values, and the captured natural pattern can be used as a reference image. If is not equal to , it means that the interpolation bias estimations of this pattern are different from the predefined values and the exposure time should be adjusted. Both image collection and the value of interpolation bias kernel will be reperformed in the same way as stated in the foregoing step. The process is ended until approximately equals . During above process, if the illumination changes more than 50%, the exposure time should be adjusted as well. For clarity, Figure 5 presents a simple flowchart describing the fundamental principle of the proposed exposure adjustment scheme.

Actual natural patterns of the Song Dynasty shipwreck are utilized to verify the correctness of proposed exposure adjustment scheme for deformation monitoring system. Four different locations are chosen as regions of interest, where the illuminations of surfaces are different. These four exact locations of the region of interest are marked in Figure 6; the corresponding illuminations are listed in Table 2. Considering the characteristics of low contrast and anisotropy of natural grains of shipwreck, the initial value of exposure time is set to 200ms, and is set to the initial value of 1.3e-2. According to the exposure adjustment scheme (Figure 5) mentioned above, four images have been recorded by the CMOS cameras with different exposures. The exposure time of each appropriate reference natural pattern is shown in the final column of Table 2. The four natural patterns shown in Figure 7 are utilized for further analysis. Figure 8 illustrates their corresponding power spectrums.

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

Firstly, the interpolation bias of DIC was evaluated. Because it is difficult to perform regular translation experiments for shipwreck, 20 deformed images were also produced by shifting in Fourier space, which is suggested to be used since it will not introduce extra bias. For each pattern, quintic B-spline was employed, respectively, to estimate the gradients of the reference image and to interpolate the subpixel registration of the deformed images. Then, the interpolation bias was estimated by (26). The value of for each pattern is calculated and given in Table 3.

The interpolation biases by DIC and by (26) are illustrated in Figure 9. Figure 9 indicates that the theoretical estimations show good agreement with the DIC results; the interpolation biases of the four natural patterns almost coincide with each other, which satisfy the requirements of the deformation monitoring system to the consistency of bias errors. Standard deviation errors of measured displacement for the four natural speckle patterns were shown in Figure 10.

#### 4. Discussion

Actual natural patterns obtained from CMOS cameras are corrupted by noise unavoidably. The interpolation bias estimation is derived in the absence of noise. Thus further study needs to verify whether the adopted criterion is still valid in the presence of noise. To address this issue, numerical experiment based on natural patterns was conducted. In most cases, both the reference image and the target images contain noise. Natural pattern C shown in Figure 7 was utilized in this numerical experiment. A series of deformed patterns is generated by shifting in Fourier space along the-axis byunits, whereranges from 0 to 1, with an increase of 0.05 pixels. The noise is assumed to be additive Gaussian white noise with zero-mean and standard deviation, where standard deviations are 0 (noise-free), 0.01, 0.02, and 0.04, respectively. Both the interpolation biases by DIC and by theoretical analysis when natural patterns are noisy are shown in Figure 11. Figure 11 indicates that the interpolation bias increases as noise level increases; the numerical experiment results show good agreement with the theoretical estimations; the theoretical criterion adopted can successfully estimate the interpolation biases in the presence of image noise.

During the superimposition of the noise, the gray values of the speckle images are rescaled to prevent the truncation of the intensity at some pixels. This rescaling process produces changes in the strength and distribution of the noise signal as well as the original speckle image. Noise can change gray distribution of image, so the power spectrums of image contained noise are altered. As mentioned above, the variations of power spectrum caused the changes of interpolation bias.

#### 5. Conclusion

A deformation monitoring system based on 2D-DIC for cultural relics protection is established in this work. Compared with regular DIC measurement system, the deformation monitoring system consists of four or more industrial cameras for simultaneously monitoring multiple locations of the cultural relics. In order to synthetically evaluate the results of deformation monitoring, the same bias errors are required for different monitoring locations. Considering the fact that only natural patterns can be used for correlation calculation and the particular illumination requirements in museum, it requires that the deformation monitoring system must have the capability to quantify bias errors and to adjust the quality of natural patterns in order to meet the consistency requirements of bias errors. To tackle this issue, a sinusoidal approximation for the interpolation bias of the IC-GN based DIC method is derived, and a speckle pattern assessment criterion is presented. Based on these theoretical analyses, exposure adjustment scheme of the deformation monitoring system is proposed. The correctness of these approaches is verified by numerical experiments using actual natural patterns. Furthermore, the theoretical criterion remains valid in the presence of noise. Thus, it is believed that the exposure adjustment scheme is promising for developing flexible and robust in situ deformation monitoring systems for use in museum environment, making 2D-DIC technique a really useful tool for in situ long-term monitoring of cultural relics. In addition to the monitoring of shipwreck in the museum environment in our study, the proposed monitoring system also is applicable to measure and monitor sample in other environment with low and varying illuminance.

#### Data Availability

The source codes and images used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This work was supported by Fundamental Research Fund for the Central Universities (no. 2018ZY08) and the National Natural Science Foundation of China (Grant no. 11502022).