**2.2.1 Reconstruction of NDVI and Albedo Data Annual Scale **

The change of vegetation coverage is the most intuitive manifestation of desertification. It is necessary to evaluate the degree of desertification under the most lush vegetation time in a year. Therefore, based on the annual desertification index constructed by the annual NDVI and the annual Albedo, it is necessary to use the annual maximum of NDVI and the annual minimum of Albedo as the basic data. The calculation process includes the following aspects:

（1）Albedo Data Calculation

The conversion of MODIS narrow-band albedo to wide-band albedo is based on the algorithm of Liang

^{[8]}^{[9]} . The short-wave albedo calculation method is used. The formula is as follows:

\({\alpha }_{short}=0.16{\alpha }_{1}+0.291{\alpha }_{2}+0.243{\alpha }_{3}+0.116{\alpha }_{4}+0.112{\alpha }_{5}+0.081{\alpha }_{6}-0.0015\) (1)

Among them, αshort is short-wave albedo, α1–α6 represents bands 1–6 in MCD43A3, respectively. In this study, white space albedo is used to calculate Albedo, because white space albedo is the integral of each incident angle, which is closer to the surface albedo in general meaning.

（2）Calculation of Annual NDVI Data

Maximum Value Compositing (MVC) is used. This method is mainly used for remote sensing data processing, and is mainly used for data analysis and reconstruction of pixels in a certain period of time. Calculation formula:

\({NDVI}_{i}=\underset{1\le \mathrm{j}\le \mathrm{n}}{\mathrm{max}}{NDVI}_{ij}\) (2)

Where i is the pixel name, j is the time point of the [1, n] time interval, and NDVIij refers to the NDVI value of the pixel i at the time j.

（3）Calculation of Albedo Annual Data

The annual reconstruction of Albedo data uses the minimum synthesis method. The method selects the minimum value of the pixel as the new pixel value within a specified period of time. The calculation formula of the annual Albedo is:

\({Albedo}_{i}=\underset{1\le \mathrm{j}\le \mathrm{n}}{\mathrm{min}}{Albedo}_{ij}\) (3)

Where i is the pixel name, j is the time point of the time interval [1, n], and Albedoij refers to the Albedo value of the pixel i at the time j.

**2.2.2 Constructing Desertification Index based on Albedo and NDVI **

NDVI values are directly proportional to vegetation coverage, while Albedo is inversely proportional to vegetation coverage. According to the positive correlation between NDVI and vegetation coverage, and the negative correlation between NDVI and Albedo, the two-dimensional （2D）spatial feature map of Albedo-NDVI can be obtained. In this paper, taking 2010 data as an example, the NDVI and Albedo data are normalized firstly. The normalization processing formula is as follows:

\(\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}=\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}/10000\) (4)

\(\mathrm{A}\mathrm{l}\mathrm{b}\mathrm{e}\mathrm{d}\mathrm{o}=\mathrm{A}\mathrm{l}\mathrm{b}\mathrm{e}\mathrm{d}\mathrm{o}/1000\) (5)

Then, 1500 sample points are evenly selected in the study zone, and NDVI is used as the X axis, and Albedo is used as the Y axis to construct the NDVI-Albedo feature space. The feature space scatter diagram is shown in Figure 2. According to the results of linear statistical regression analysis, the negative correlation between Albedo and NDVI can be expressed by the following linear relationship:

\(Albedo=-0.2303\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}+0.3258\) (6)

**Figure 2
**Two-dimensional (2D) feature space map of Albedo-NDVI in China-Pakistan Economic Corridor 2010As shown in Fig. 2, in the feature space constructed by Albedo and NDVI, the different positions of the negative correlation between Albedo and NDVI represent the state and degree of different stages of desertification, and the degree of desertification increases with the decrease of NDVI value and is up with the rise of Albedo value, that is, the negative correlation linear expression of Albedo and NDVI can reflect the trend of desertification. According to the conclusions of Verstraete and Pinty

^{[10]}, if the Albedo-NDVI feature space is divided in the vertical direction representing the trend of desertification, different desertification land can be effectively distinguished as shown in Figure 3, and the desertification difference index model DDI is used to represent:

\(DDI=\alpha ×NDVI-Albedo\) (6)

**Fig.3
**Graphical expression of desertification difference index (DDI) in NDVI-Albedo spaceIn a specific application, the constant α can be determined from the slope in the Albedo-NDVI feature space, where k is the slope of the characteristic equation.

\(\alpha ×k=-1\) (7)

According to the above calculation results, the slope of the negative correlation of Albedo-NDVI is k=-0.2303, then α=4.3422, and the expression of DDI is as follows:

\(DDI=4.3422×NDVI-Albedo\) (8)

**2.2.3 ****Desertification Categorization**

In order to meet the needs of desertification evaluation and mapping, in 1984, the United Nations Food and Agriculture Organization (FAO) and the United Nations Environmental Programme (UNEP) in the “Desertification Assessment and Mapping Program” came up with the specific quantitative vegetation of desertification status, development rate and intrinsic risk assessment from vegetation degradation, wind erosion, water erosion and salinization as four aspects, and desertification was divided into mild, moderate, severe and extremely severe (level) according to the degree of development.

Through the DDI model constructed in the previous section, the dedicated data on the desertification index of the China-Pakistan Economic Corridor from 2000 to 2017 was obtained. When making the desertification categorization based on DDI, most researchers use the natural fracture method

^{[7]}^{[11]}^{[12]} in ArcGIS reclassification. The natural fracture method is a categorization point based on the Jank optimal method of statistics

^{[13]}. This classification minimizes the difference within the category and maximize the difference among the categorized. In this paper, the natural fracture method and the desertification degree quartering method are also used. The DDI value is divided into five intervals, and the five-level categorization index is determined (Table 2). The user can also combine with the field survey data to further fine-tune the DDI categorization table.

**Table 2
**DDI categorization table**Type** | **DDI Values** |

**Snowy and icy waters** | DDI≤−0.26 |

**Severe desertification** | −0.26＜DDI≤0.12 |

**Moderate desertification** | 0.12＜DDI≤0.55 |

**Mild desertification** | 0.55＜DDI≤1.6 |

**Non desertification** | 1.6＜DDI≤4.2 |

**2.2.4 Simulation Process of Dedicated Data for Desertification Index **

Firstly, the NDVI and Albedo data are reconstructed, the annual NDVI value is calculated by the maximum synthesis method, and the annual Albedo value is calculated by the minimum synthesis method. Then, the Albedo-NDVI feature space is constructed by fitting the Albedo and NDVI data to obtain the fitting slope k. Finally, the desertification index is constructed, the desertification grade is divided, and the desertification classification map based on the desertification index is completed. The data processing flow is shown in Figure 4.

**Fig. 4
**Simulation process for desertification index’s dedicated data