Detecting and Assessing Nondominant Farmland Area with Long-Term MODIS Time Series Images

: While most land use and land cover (LULC) studies have focused on modeling, change detection and driving forces at the class or categorical level, few have focused on the subclass level, especially regarding the quality change within a class such as farmland. The concept of nondominant farmland area (NAF) is proposed in this study to assess within class variability and quantify farmland areas where poor environmental conditions, unsuitable natural factors, natural disasters or unsustainable management practices lead to poor crop growth and thus low yield. A 17-year (2000–2016) time series of the Normalized Di ﬀ erence Vegetation Index (NDVI) was used to develop a NAF extraction model with abnormal features in the NDVI curves and subsequently applied to Heilongjiang province in China. The NAF model was analyzed and assessed from three aspects: agricultural disasters, soil types and medium- and low-yield ﬁelds, to determine dominant factors of the NAF patterns. The results suggested that: (1) the NAF model was able to extract a variety of NAF types with an overall accuracy of ~80%. The NAF area accumulated more than 8 years in 17 years is 6.20 thousand km 2 in Heilongjiang Province, accounting for 3.75% of the total cultivated land area; (2) the NAF had signiﬁcant spatial clustering characteristics and temporal variability. 53.24% of the NAF accumulated more than 8 years in 17 years is mainly concentrated in the west of Heilongjiang Province. The inter-annual NAF variability was related with meteorological variations, topography and soil properties; and (3) the spatial and temporal NAF patterns seem to reﬂect a cumulative impact of meteorological disasters, poor farmland quality, and soil degradation on crop growth. The determinant factors of the observed NAF patterns di ﬀ ered across regions, and must be interpreted in the local context of topography, soil properties and meteorological environment. Spatial and temporal NAF variability could provide useful, diagnostic information for precision farmland management.


Introduction
The world's population has reached 7.7 billion (United Nations [1]), and there is a substantial demand for food but there is a limited amount of reclaimable farmland. Therefore, increasing the yields of middle-and low-yield farmland is probably the most feasible means to grow more food to meet the increasing demand. Low-yield farmland is defined as lands with a lower or much lower 3 a certain growth period or the whole growing period regarding the growth, canopy structure, biomass, and leaf area index that differ from those of normal crops. These differences could be extracted by remote sensing monitoring, and the extraction range can be analyzed to identify the limiting factors of the NAF.
Hence, the overarching goal of this study was to build a NAF model based on remotely sensed imagery to better understand the spatial and temporal dynamics of farmlands. The constructed NAF model was then validated to assess its accuracy from three aspects: agricultural disasters, soil types and medium-and low-yield fields. The model was applied in a case study in Heilongjiang Province, China, to determine the driving factors by observing the NAF patterns in the region.

Study Area
Heilongjiang Province (Figure 1), located in the northeastern China, is an important grainproducing area and commodity grain base in China. Heilongjiang Province has a total area of 42.66 million hectares, accounting for 4.7% of the total national land area, and 3378.74 million hectares of agricultural land, accounting for 83.5% of the total land area of the province. In 2016, the sown area and crop yield of Heilongjiang Province accounted for 10.44% and 9.8% of the total in China, respectively. Songnen Plain within the province is formed by the alluviation of the Songhua River and Nenjiang in the western portion of Heilongjiang Province, and the midwestern portion of Songnen Plain belongs to the aeolian sandy soil zone. The region has experienced many years of drought. The northeastern portion of Heilongjiang Province is Sanjiang In 2016, the sown area and crop yield of Heilongjiang Province accounted for 10.44% and 9.8% of the total in China, respectively. Songnen Plain within the province is formed by the alluviation of the Songhua River and Nenjiang in the western portion of Heilongjiang Province, and the midwestern portion of Songnen Plain belongs to the aeolian sandy soil zone. The region has experienced many years of drought. The northeastern portion of Heilongjiang Province is Sanjiang Plain, which is surrounded by mountains on three sides. Xiao Hinggan Ling branch of the Green Montenegro is to the west, the Wandashan branch vein of the Fen Shui Gang is to the south and the Nadanhada Ridge of the Wandashan main vein reaches to the east. Heilongjiang Province has a temperate continental monsoon climate, four distinct seasons, a short summer, a long winter, and a large seasonal temperature difference. The annual average temperature is between −5 and 5 • C, the annual average precipitation is 400-600 mm, and the crop growth period is typically within the range of 122-274 days. The main crops grown in Heilongjiang Province are corn, soybean, rice, potato, wheat and sugar beet.

Data Acquisition
The remote sensing data used in this study are 250 m MODIS 8-day synthetic data (MOD09Q1 and MYD09Q1) from 2000 to 2016 and 250 m MODIS 16-day synthetic data (MOD13Q1). (Four MODIS 8-day data images that were of poor quality were replaced by the 16-day data.) There were two phases per year (day of year [DOY] 185 and 193) for a total of 34 phases of remote sensing images (see http://reverb.echo.nasa.gov/reverb). The Net Primary Productivity (NPP) and spatial distribution of the soil texture in China were obtained from 2000 to 2010. These datum set was provided by Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) (http://www.resdc.cn). Global 30-meter land cover data (see http://www.globallandcover.com) and digital elevation model (DEM) data at a 90-meter resolution were obtained for Heilongjiang Province (see http://www.gscloud.cn). The MODIS 8-day reflectance products and 16-day NDVI (Normalized Difference Vegetation Index) products were converted to projections using MODIS Reprojection Tool (MRT), and then MODIS 8-day reflectance products were computed in ArcGIS (geographic information system) toolbox to obtain the NDVI.
The meteorological data (temperature, precipitation) were obtained from 11 agrometeorological observation stations in Heilongjiang Province from 2000 to 2016, the meteorological data of agricultural meteorological disasters from 2000-2011, the crop growth development and the farmland soil moisture in Heilongjiang Province from 2000 to 2013 were obtained from the meteorological data center of China Meteorological Administration (see http://data.cma.cn/user/toLogin.html). The China meteorological spatial interpolation dataset for 20002015 were obtained from the Chinese Academy of Sciences, Resources and Environmental Science Data Center (see http://www.resdc.cn).
Agricultural insurance companies provided insurance data for 182,913, 180,208, 18,359, and 221,362 cropland patches from 2013-2016 through the digital imaging of high-resolution spatial images. The land survey confirmed that the accuracy of the land insurance data was higher than 93%, based on the data for 24 crop types. A total of 216 and 779 disaster points were obtained in 2015 and 2016, respectively ( Figure A1a). In 2017, 58 aeolian sandy soil samples were collected [27].

Determination of the Critical Period for the Remote Sensing Images
The agricultural meteorological datasets of Heilongjiang Province contain 30 coerced sites from 2000 to 2011 (Table A1). The results showed that 87.8% of the crops were stressed before 19 July (DOY 201) in 10 typical counties distributed throughout the province.
Climatic differences lead to different crop growing determinations in Heilongjiang Province. There are differences when sowing in the morning and evening and between early-and late-maturing varieties, and phenology zoning could be responsible for this difference [28]. The main crop phenology (Table 1) shows that the late-sowing crop was in the early growth stage on DOY 177, and the early maturing variety was mature on DOY 201. Therefore, the accuracy of crop growth will be limited by monitoring the status of crop growth before DOY 185 or after DOY 201. At the same time, the disasters mainly occurred in the early period, and early monitoring could better characterize the farmland quality. Therefore, it is reasonable to use the remote sensing data of DOY 185 and 193 to monitor the crop growth.  Province. SW, SD, HD, MT, FL, PD, JT,  TS, ML, TP and TL represent Sowing

Phenology Zoning of Farmland
Crop phenology refers to the regular growth and change in crops in a certain environment and shows the temporal variation in crop growth during the year. Crops changes in the timing of leaf out, leaf senescence, and flowering according to environmental cues [29,30], significantly impact ecosystem structure and function [31]. The phenological characteristics of the crops in different regions are diverse, and the crop growth is discrepant [32].
The geography and climatic conditions in the study area are quite different. It is difficult to accurately determine the NAF through unified monitoring of the entire study area; therefore, it was necessary to divide the study area into several regions with a uniform phenology. According to the variations in the time series NDVI curves, the phenology of the relatively homogeneous region of Heilongjiang Province was extracted using the phenological parameters by nearest neighbor classification method [33].

NAF Extraction Modeling
NAF refers to an area with poor environmental conditions, unsuitable crop environmental factors, a low yield, disaster-prone farmland or land that is not suitable for crop growth. NAF not only refers to the quality of farmland problems but also incorporates climate suitability issues, including vulnerability to meteorological disasters, poor farmland quality, low-yield fields, soil degradation and other zones. NAF formation is mainly affected by meteorology, soil, topography, crop variety selection, tillage methods and management measures.
The NAF was extracted using the vegetation index spatial statistical analysis method. Many researchers have used the anomaly index to monitor crop growth and disasters. This method extracts poor growth regions by comparing the same pixel in different years. If a region has experienced poor growth for several years, it is difficult to extract the NAF. NDVI is the most widely used index for crop growth monitoring [34]. According to the zoning results of the farmland phenology in Heilongjiang Province, the NDVI of each phenological region was statistically analyzed, and the NAF was extracted by the NDVI mean value and standard deviation. Figure 2 is a flowchart in this study.
NDVI values in the phenological region are normally distributed ( Figure A2), and the average NDVI value in the phenological region is calculated. The ratio of NDVI value to NDVI mean value in any pixel is calculated, and the standard deviation is graded on the comparison value. The pixel of the ratio less than the mean value minus the standard deviation belong to the NAF: where G is the crop growth index and NDVI M is the average NDVI in each phenological region.
Remote Sens. 2020, 12, 2441 6 of 22 where is the threshold of the NAF, which is less than the minimum value of the normal growth crop of in the phenology region, is the average of this phenological region, and is the standard deviation of this phenological region.

Vector Angle Method
The vector angle method was also applied to an n-dimensional vector. Assuming that and are two n-dimensional vectors, is [ , , ..., ] and is [ , , ..., ], the cosine of angle between and is equal to the following expression: where is the interannual variation in the NAF area in a county, is the interannual variation in the NAF area in another county, = 1 is the year 2000, and n is the cumulative year of the study.
When the cosine of the angle between the two curves is equal to 1, the two curves are completely coincident. When the cosine of the angle is close to 1, the two curves are similar and can be classified into one class. The smaller the angle cosine is, the more irrelevant the two curves are. This paper assumes that the interannual variation in the NAF area of a county is a vector (such as Nehe 2000Nehe -2001Nehe , 2001Nehe -2002. for the Nehe county vector). The vector angle method was used to classify the region with the NAF change trend in 80 counties of Heilongjiang Province from 2000-2016. The G MIN values of normal growth regions (NGR) in each phenological region (Equation (2)) were higher than the threshold of NAF, indicating that monitoring the NAF using DOY 185 and 193 data could better exclude the monitoring differences caused by the different crops: where G MIN is the threshold of the NAF, which is less than the minimum value of the normal growth crop of G in the phenology region, G MEAN is the average G of this phenological region, and G STD is the G standard deviation of this phenological region.

Vector Angle Method
The vector angle method was also applied to an n-dimensional vector. Assuming that A and B are two n-dimensional vectors, A is [A 1 , A 2 , . . . , A n ] and B is [B 1 , B 2 , . . . , B n ], the cosine of angle θ between A and B is equal to the following expression: where A i is the interannual variation in the NAF area in a county, B i is the interannual variation in the NAF area in another county, i = 1 is the year 2000, and n is the cumulative year of the study.
When the cosine of the angle between the two curves is equal to 1, the two curves are completely coincident. When the cosine of the angle is close to 1, the two curves are similar and can be classified into one class. The smaller the angle cosine is, the more irrelevant the two curves are. This paper assumes that the interannual variation in the NAF area of a county is a vector (such as Nehe 2000Nehe -2001Nehe , 2001Nehe -2002. for the Nehe county vector). The vector angle method was used to classify the region with the NAF change trend in 80 counties of Heilongjiang Province from 2000-2016.

Correlation Analysis
Crop growth is closely related to weather conditions factors such as temperature and precipitation that affect all crop stages from planting to harvest [35]. The average temperature and precipitation in the main growing period of each year (early April to early July) were compared with the NAF area in the typical counties in Heilongjiang Province, and their correlation was obtained.
NPP represents the accumulated organic matter of plants per unit area and time and is the total primary production minus the remaining organic or energy of plant respiration. Therefore, NPP can represent crop growth and verify the extraction accuracy of the NAF. The correlation between the NAF and NPP was used to verify the extraction accuracy of the NAF. According to the definitions of NAF and NPP, if the NAF is large, then the NPP value should be small. Therefore, Corn belt of Songnen Plain was obtained according to the planting structure of Heilongjiang Province extracted by Liu et al. [33].

Accuracy Assessment
Based on the field-based field disaster acquisition, the agricultural insurance companies provided 216, 779 disaster samples in 2015 and 2016, and the overall NAF accuracy was verified. Yan et al. [2] extracted a 5.61 ten thousand km 2 low-yield field in Heilongjiang using the NPP of farmland from 2000 to 2008. The NAF area was extracted in this paper and compared with that of Yan et al. Fifty-eight aeolian sandy soil samples collected in the field in 2016 were used to verify the NAF caused by soil-limiting factors: where P is the overall accuracy, P C is the number of correct categories and N is the total number of sample points.

Model Validation
To validate the accuracy of the threshold, the sample points of three crops (soybean, corn and rice) on three typical farms (Keshan, Yaluhe and Hailin farm) were collected in 2016, which were 56, 209 and 42, respectively. The average and standard deviation of G on DOY 185 and 193 NGR of the three crops were counted to calculated G MIN (Equation (2)

Model Validation
To validate the accuracy of the threshold, the sample points of three crops (soybean, corn and rice) on three typical farms (Keshan, Yaluhe and Hailin farm) were collected in 2016, which were 56, 209 and 42, respectively. The average and standard deviation of G on DOY 185 and 193 NGR of the three crops were counted to calculated GMIN (Equation (2) Figure 4 is spatial distribution of NAF in Heilongjiang Province from 2000 to 2016. The NAF area accumulated more than 8 years in 17 years is 6.20 thousand km 2 in Heilongjiang Province, accounting for 3.75% of the total cultivated land area. The NAF area from 2000 to 2008 extracted in this study was 5.55 ten thousand km 2 in Heilongjiang Province, which is almost consistent with 5.61 ten thousand km 2 low-yield field extracted by Yan et al. [2]. In addition, according to the NPP threshold of low-yielding fields set by Yan et al., 2.59 ten thousand km 2 of low-yielding fields were extracted by using 2009 and 2010 NPP data, which is roughly consistent with the 2.41 ten thousand km 2 NAF extracted in this paper. The correlation coefficient between the NAF area of the corn belt on Songnen Plain and the NPP was −0.57. At the same time, the overall accuracy of the disaster verification points in 2015 and 2016 were respectively 82.87% and 79.33% and that of aeolian sandy soil was 91.38%. The above results demonstrated that the NAF model can extract low-yielding fields commendably.

NAF Patterns
The spatial distribution of the NAF in Heilongjiang Province was showed according to the NAF accounting for the proportion of cultivated land from 2000-2016 ( Figure 5). The NAF proportions on the southwest of Songnen Plain, including Durbert, Lindian, Daqing, Tailai, Anda, Gannan and Fuyu county, were more than 15%. 53.24% of the NAF accumulated more than 8 years in 17 years is mainly concentrated in the west of Heilongjiang Province, due to mainly located in the aeolian sandy soil region that have a poor resistance to drought and flooding. However, the NAF proportions in the southwest of Songnen Plain, including Suihua, Qing'an, Hulan, Mingshui Suileng and Wangkui county, were less than 3%, because of being mainly located in the black soil region with good-quality farmland. Moreover, the NAF proportions of Sanjiang Plain was lower overall except Jixi county being located among mountains and prone to floods during the rainy season ( Figure A4).
Songnen Plain is one of the two major grain production bases in Heilongjiang Province. By superposing the NAF of similar years on Songnen Plain, we obtained three spatial distribution patterns of NAF (  Figure A5c), which was more than other patterns. Although the precipitation ranges of Pattern 1 and Pattern 3 are roughly the same, the spatial distribution is significantly different. The precipitation of pattern 1 gradually decreases from southwest to northeast ( Figure A5a), but the precipitation of pattern 3 gradually decreases from northwest to southeast ( Figure A5c). 8 thousand km 2 low-yield field extracted by Yan et al. [2]. In addition, according to the NPP threshold of low-yielding fields set by Yan et al., 2.59 ten thousand km 2 of low-yielding fields were extracted by using 2009 and 2010 NPP data, which is roughly consistent with the 2.41 ten thousand km 2 NAF extracted in this paper. The correlation coefficient between the NAF area of the corn belt on Songnen Plain and the NPP was −0.57. At the same time, the overall accuracy of the disaster verification points in 2015 and 2016 were respectively 82.87% and 79.33% and that of aeolian sandy soil was 91.38%. The above results demonstrated that the NAF model can extract low-yielding fields commendably.  concentrated in the west of Heilongjiang Province, due to mainly located in the aeolian sandy soil region that have a poor resistance to drought and flooding. However, the NAF proportions in the southwest of Songnen Plain, including Suihua, Qing'an, Hulan, Mingshui Suileng and Wangkui county, were less than 3%, because of being mainly located in the black soil region with good-quality farmland. Moreover, the NAF proportions of Sanjiang Plain was lower overall except Jixi county being located among mountains and prone to floods during the rainy season ( Figure A4).  10 gradually decreases from southwest to northeast ( Figure A5a), but the precipitation of pattern 3 gradually decreases from northwest to southeast ( Figure A5c).

NAF Classification
The vector angle method was used to classify the regions where the variation tendency in the NAF was consistent from 2000 to 2016 in 80 counties of Heilongjiang Province. when the cosine values between two counties are greater than 0.88, they belong to the same category. Therefore, 80 counties were divided into 11 categories (Figure 8), clearly demonstrating that almost all adjacent counties in 10 gradually decreases from southwest to northeast ( Figure A5a), but the precipitation of pattern 3 gradually decreases from northwest to southeast ( Figure A5c).

NAF Classification
The vector angle method was used to classify the regions where the variation tendency in the NAF was consistent from 2000 to 2016 in 80 counties of Heilongjiang Province. when the cosine values between two counties are greater than 0.88, they belong to the same category. Therefore, 80 counties were divided into 11 categories (Figure 8), clearly demonstrating that almost all adjacent counties in

NAF Classification
The vector angle method was used to classify the regions where the variation tendency in the NAF was consistent from 2000 to 2016 in 80 counties of Heilongjiang Province. when the cosine values between two counties are greater than 0.88, they belong to the same category. Therefore, 80 counties were divided into 11 categories (Figure 8), clearly demonstrating that almost all adjacent counties in the same category. For instance, the five counties of Qiqihar, Tailai, Durbert, Daqing and Anda on Songnen Plain were classified into the same category, and the three counties of Tongjiang, Fujin and Raohe on Sanjiang Plain were classified into the same category as well. Moreover, four categories were selected from 11 categories to map the interannual variation in the NAF in similar regions from 2000 to 2016 ( Figure A7). The curve trends of Nehe, Gannan, Yi'an, Fuyu and Baiquan county were the most consistent, as seen in Figure A7a, and the curves in Figure A7b-d showed roughly the same trend.
Remote Sens. 2020, 04, x FOR PEER REVIEW 11 of 24 11 the same category. For instance, the five counties of Qiqihar, Tailai, Durbert, Daqing and Anda on Songnen Plain were classified into the same category, and the three counties of Tongjiang, Fujin and Raohe on Sanjiang Plain were classified into the same category as well. Moreover, four categories were selected from 11 categories to map the interannual variation in the NAF in similar regions from 2000 to 2016 ( Figure A7). The curve trends of Nehe, Gannan, Yi'an, Fuyu and Baiquan county were the most consistent, as seen in Figure A7a, and the curves in Figure A7b, A7c, and A7d showed roughly the same trend.

NAF Pattern Formed Factors
According to our knowledge, we analyzed the formed factors of NAF patterns from soil properties and meteorological conditions. Figure 9a clearly showed the NAF of Songnen Plain distribution was mainly in the southwest. The NAF area accumulated more than 8 years in 17 years is 4.22 thousand km 2 on Songnen Plain. Combining the field samples and the soil-sand distribution (Figure 9b), most of the soil types in the NAF of the southwestern Songnen Plain were saline-alkali soil, aeolian sandy soil and sandy soil, which are low organic matter and weak water holding capacity [27], leading to the long-term poor crop growth in these areas. At the same time, the NAF accounted for more than 10% of the total farmland area ( Figure 5), indicating that the soil texture has a greater influence on the NAF in Duerbote, Daqing, Tailai, Gannan, and Longjiang county than the meteorological conditions.

NAF Pattern Formed Factors
According to our knowledge, we analyzed the formed factors of NAF patterns from soil properties and meteorological conditions. Figure 9a clearly showed the NAF of Songnen Plain distribution was mainly in the southwest. The NAF area accumulated more than 8 years in 17 years is 4.22 thousand km 2 on Songnen Plain. Combining the field samples and the soil-sand distribution (Figure 9b), most of the soil types in the NAF of the southwestern Songnen Plain were saline-alkali soil, aeolian sandy soil and sandy soil, which are low organic matter and weak water holding capacity [27], leading to the long-term poor crop growth in these areas. At the same time, the NAF accounted for more than 10% of the total farmland area ( Figure 5), indicating that the soil texture has a greater influence on the NAF in Duerbote, Daqing, Tailai, Gannan, and Longjiang county than the meteorological conditions. In order to understand the effect of precipitation and temperature on NAF, the correlation coefficients between the NAF and the meteorological data of 9 typical counties from 2000 to 2016 shown in Table 2.
The NAF of three typical counties (Yilan, Baoqing and Fuyuan county) on the Sanjiang Plain was more related with precipitation, yet it was not significant correlation with temperature, demonstrating precipitation was the main factor of NAF formation. For example, the average precipitation of Yilan county was 194.1 mm during the main growing period from 2000 to 2016. The precipitation was lower in 2001 (73.7 mm) than that in 2002 (299.3 mm), but the NAF in 2001 was higher than that in 2002 (Figure 10a,b).
Compared with Sanjiang Plain, the NAF of three typical counties (Gannan, Longjiang and Lanxi county) on Songnen Plain was affected by different factors. For instance, the correlation between precipitation and temperature with NAF in Longjiang county was 0.57 and 0.5 (Table 2), respectively, however, that in Gannan and Lanxi county were −0.38 and 0.07, and −0.1 and 0.5 (Table 2), respectively.
The NAF were more related with temperature on northwest of Xiao Hinggan Ling. The correlation between NAF and temperature in Nenjiang and Heihe county was 0.38 and 0.24 (Table  2), respectively. Taking Nenjiang county as an example, the average temperature of Nenjiang county during the main growing period from 2000 to 2016 was 13.63 °C. The temperature was higher (14.32 °C) in 2002 than that in 2013 (12.95 °C), but the NAF in 2002 was lower than that in 2013 ( Figure  11a,b).
Nevertheless, southeast of Xiao Hinggan Ling is unlike any other place in Heilongjiang Province. the correlation between precipitation and temperature with NAF in Ning'an county was −0.14 and −0.06 (Table 2), respectively, indicating the NAF of Ning'an county was neither affected by precipitation nor by temperature. As expected, the elevation of profile fluctuates greatly in Ning'an county ( Figure A4b). Hence, topography was the main factor of NAF formation for southeast of Xiao Hinggan Ling. In order to understand the effect of precipitation and temperature on NAF, the correlation coefficients between the NAF and the meteorological data of 9 typical counties from 2000 to 2016 shown in Table 2. Table 2. Correlation coefficients between NAF of typical sites and meteorological data. NW and SE respectively represent northwest and southeast of Xiao Hinggan Ling. P: Precipitation from April to early July; T: average temperature from April to early July; "*" and "**" represent a significant difference at the 0.05 level and that at the 0.01 level, respectively. The NAF of three typical counties (Yilan, Baoqing and Fuyuan county) on the Sanjiang Plain was more related with precipitation, yet it was not significant correlation with temperature, demonstrating precipitation was the main factor of NAF formation. For example, the average precipitation of Yilan county was 194.1 mm during the main growing period from 2000 to 2016. The precipitation was lower in 2001 (73.7 mm) than that in 2002 (299.3 mm), but the NAF in 2001 was higher than that in 2002 (Figure 10a,b).

Sanjiang Plain
Compared with Sanjiang Plain, the NAF of three typical counties (Gannan, Longjiang and Lanxi county) on Songnen Plain was affected by different factors. For instance, the correlation between precipitation and temperature with NAF in Longjiang county was 0.57 and 0.5 (Table 2), respectively, however, that in Gannan and Lanxi county were −0.38 and 0.07, and −0.1 and 0.5 (Table 2), respectively.
Remote Sens. 2020, 12, 2441 13 of 22 13 early July; T: average temperature from April to early July; " * " and " ** " represent a significant difference at the 0.05 level and that at the 0.01 level, respectively.    The NAF were more related with temperature on northwest of Xiao Hinggan Ling. The correlation between NAF and temperature in Nenjiang and Heihe county was 0.38 and 0.24 (Table 2), respectively. Taking Nenjiang county as an example, the average temperature of Nenjiang county during the main growing period from 2000 to 2016 was 13.63 • C. The temperature was higher (14.32 • C) in 2002 than that in 2013 (12.95 • C), but the NAF in 2002 was lower than that in 2013 (Figure 11a,b).
Remote Sens. 2020, 04, x FOR PEER REVIEW 13 of 24 13 Table 2. Correlation coefficients between NAF of typical sites and meteorological data. NW and SE respectively represent northwest and southeast of Xiao Hinggan Ling. P: Precipitation from April to early July; T: average temperature from April to early July; " * " and " ** " represent a significant difference at the 0.05 level and that at the 0.01 level, respectively.   Nevertheless, southeast of Xiao Hinggan Ling is unlike any other place in Heilongjiang Province. the correlation between precipitation and temperature with NAF in Ning'an county was −0.14 and −0.06 (Table 2), respectively, indicating the NAF of Ning'an county was neither affected by precipitation nor by temperature. As expected, the elevation of profile fluctuates greatly in Ning'an county ( Figure A4b). Hence, topography was the main factor of NAF formation for southeast of Xiao Hinggan Ling.

Discussion
In this paper, remote sensing images according to the main periods of crop stress and the main period of crop phenology were selected to build an NDVI ratio model in the different phenology regions to evaluate the growing conditions of crops, effectively avoiding phenological differences. At the same time, The G MIN values of NGR in each phenological region were higher than the threshold of NAF, indicating that monitoring the NAF using the model could better exclude the monitoring differences caused by the different crops. Previous researchers have mainly used LUE, temperature, moisture and leaf phenology to build productivity models and then estimate NPP, which is a laborious and time-consuming process [36]. In this paper, a variety of unfavorable factors were considered synthetically, eliminating the influence of crop types, phenology and topography. This method provides a quantitative way to evaluate the internal variations in the same land use type and the spatial and temporal patterns of internal farmland quality.
The spatial distribution of NAF in Heilongjiang Province was mostly aggregated (Figures 5 and 8), indicating that soil conditions, terrain and climatic factors play a decisive role in the spatial distribution. For example, these regions of classes 1, 3, and 7 were mainly caused by the soils, and the NAF proportions of them were 20-30%, 10-20%, 7-10% ( Figure 5), respectively; the NAF caused by terrain included classes 8 and 9, and the NAF proportions of them were 0-20% and 5-10% respectively; classes 2, 10 and 11 were mainly caused by precipitation, and the NAF proportions of them account for 0-10%, 5-20%, 7-10%, respectively; both classes 5 and 6 were caused by temperature, and the NAF proportions of them were 0-7%. It is worth mentioning that the range of pattern 1 (Figure 6a) was the common part of the three patterns ( Figure 6) due to the poor soil texture in the southwestern of Songnen Plain. Additionally, the typical regions of Heilongjiang Province were selected to analyze the reasons for the NAF formation. NAF in the southwestern of Songnen Plain is formed because of the soil texture, such as Lindian, Tailai, Gannan and Longjiang county ( Figure 9). Therefore, we could adjust the planting structure according to the limiting factors and characteristics of the local region to carry out the planting and breeding structure by combining agriculture and animal husbandry. Moreover, temperature seriously impacts on crop growth in Nenjiang and Heihe county, however, the temperature is more suitable for soybean growth. Thus, the rotation between grain and bean should be implemented in eco-friendly farming systems in Nenjiang and Heihe. Unlike Songnen Plain, Sanjiang Plain, such as Baoqing, Fuyuan and Yilan county, is seriously affected by precipitation. However, the soil is more fertile and the precipitation is more adequate in this region. Hence, flood-resistant and drought-resistant crops can be planted in flood-prone and drought-prone regions. In addition, Ning'an county is surrounded by mountains on both sides, causing the topography seriously effects on crop growth. Therefore, this region could be changed to the cultivation of medicinal herbs.
Remote sensing technology has advantages to eliminate the effects of crop type, phenology and topography by directly monitoring phenology zoning of the crop growth [37,38]. According to the analysis of the influencing factors of the NAF in the typical counties (Section 3.4), the dominant factors in the different regions were discrepant for the large-scale NAF and the middle-and low-yield fields. Currently, the existing methods generally constructed index systems, and the weights of the different index systems were fixed. In fact, the contribution rate of the factors in the different regions differed. For example, the northwestern region was mainly affected by temperature, the southwestern region was mainly affected by precipitation, soil desertification, and the salinization degree, and the eastern region was affected by topography and precipitation. If the evaluation model is constructed according to a uniform weight, the result of the extraction will deviate from the actual result. Agricultural experts can analyze the formation causes (such as topography, soil and hydrology) [39][40][41] according to the different spatial positions of the NAF, to make targeted improvements and adjust the planting structure. NAF is not only affected by topography, soil, and meteorological and other environmental factors, but the adjustment of crop planting structure also leads to the changes of NAF pattern. The NAF pattern on Sanjiang Plain shown in Figure 7c was due to the large number of dry lands changed to paddy fields. It would be difficult to monitor this change by using the traditional evaluation index system considering only the role of environmental factors.
The crop planting structure could be adjusted and agricultural production management policy could be drafted by accurately extracting the spatial and temporal distribution of the NAF. The annual precipitation of Songnen Plain ( Figure A5) showed the overall precipitation pattern, i.e., Figure A5c > Figure A5a > Figure A5b, on Songnen Plain. And thus, the range of pattern 2 was an arid area. Combined with Figure 6a-c, it was found that excessive or low precipitation on Songnen Plain could increase the NAF area. Therefore, we can predict the spatial and temporal distribution of NAF in this region according to the forecast of future precipitation by meteorological experts, which is helpful for early disaster warnings and provides a scientific basis for disaster prevention and mitigation and policy adjustment. China proposed the adjustment of the planting structure in the "Sickle Bay" region to reduce corn planting areas [42]. At present, the government adjusts the planting structure mainly based on statistical data, which are one-sided. The NAF extracted in this paper not only includes the nondominant area of corn but also the nondominant area of other crops, which can be used to pertinently adjust the crop planting structure.
The NAF extracted using lower spatial resolution of MODIS data limits the final extraction accuracy. For example, the mixed pixel problem of farmland and construction land in the central portion of Harbin was monitored as the NAF ( Figure 5). Further research is needed to consider the fusion of different time and space resolution data to construct a high time and space resolution dataset to improve the accuracy. The model proposed in this study is suitable only for one crop a year, however, the multiple cropping region in one year needs to be classified for different phenological zoning, and then the nondominant area of different crops needs to be extracted.

Conclusions
In this paper, a 17-year (2000-2016) time series of NDVI was used as the basis of a NAF extraction model that was used to extract the NAF with abnormal time series NDVI curves. The model accuracy was validated from three aspects: agricultural disasters, soil types and medium-and low-yield fields. Then, the driving factors of the NAF spatial and temporal pattern were analyzed based on meteorological, soil and topographic data. The following conclusion were obtained: (i) Extracting the NAF using time series remote sensing data had high accuracy, eliminated strong subjectivity in statistics and obtained a more accurate spatial and temporal NAF distribution. Based on the meteorological conditions, the NAF can be predicted ahead of time, which can result in early warnings and provide scientific guidance for the adjustment of planting structure, disaster monitoring, agricultural production management and macro decision making.
(ii) The NAF has significant spatial and temporal variations and regional clustering characteristics. The interannual variation trend in the NAF was consistent with similar topographic, soil and meteorological conditions. Based on the phenological division to monitor crop growth, the NAF extraction model eliminates the effects of crop type, phenology and topography.
(iii) The influence factors of NAF in different regions are discrepant. The reasons of NAF formation were not only topography, soil, meteorological factors and other environmental factors, but also the adjustment of crop planting structure, which leads to the NAF change patterns.

Acknowledgments:
We gratefully acknowledge the support from Sunlight Agricultural Mutual Insurance Company for their help in the field work. We further acknowledge our group members for their help. Da Hinggan Ling (Figure 1c). Hence, Sanjiang Plain, Songnen Plain and Xiaoxing'anling each selected a farm as a sampling area. Moreover, the main crops in Heilongjiang Province are corn, soybean and rice. In 2016, corn, soybean and rice grown occupy of the total account for 53.90%, 26.82% and 16.44%, respectively. Therefore, the sample points of corn, soybean and rice 209, 56 and 42, respectively. Figure A1a is the spatial distribution map with a total of 216 and 779 disaster points obtained in 2015 and 2016, respectively. Heilongjiang Province is mainly divided into four parts: Songnen Plain, Sanjiang Plain, Xiao Hinggan Ling and Da Hinggan Ling (Figure 1b), but there is very little cultivated land on Da Hinggan Ling (Figure 1c). Hence, Sanjiang Plain, Songnen Plain and Xiaoxing'anling each selected a farm as a sampling area. Moreover, the main crops in Heilongjiang Province are corn, soybean and rice. In 2016, corn, soybean and rice grown occupy of the total account for 53.90%, 26.82% and 16.44%, respectively. Therefore, the sample points of corn, soybean and rice 209, 56 and 42, respectively. Figure A1a is the spatial distribution map with a total of 216 and 779 disaster points obtained in 2015 and 2016, respectively.     Figure A2. NDVI histogram in a certain phenological region. Figure A2. NDVI histogram in a certain phenological region.

Conflicts of
Appendix A. 5

. The Elevation Map in Study Area
The elevations of Xiao Hinggan Ling and Da Hinggan Ling are higher than other regions. There is almost no farmland on Da Hinggan Ling. Figure A4b is the red line ( Figure A4a The spatial precipitation data showed three patterns on Songnen Plain ( Figure A5). The spatial precipitation data showed three patterns on Sanjiang Plain ( Figure A6). Figure A6b    Appendix A. 6

. Spatial Precipitation Map
The spatial precipitation data showed three patterns on Songnen Plain ( Figure A5). The spatial precipitation data showed three patterns on Sanjiang Plain ( Figure A6). Figure A6b,c is precipitation of pattern 2 (Figure 7b). The spatial precipitation data showed three patterns on Songnen Plain ( Figure A5). The spatial precipitation data showed three patterns on Sanjiang Plain ( Figure A6). Figure A6b    The interannual NAF curves of four similar regions were plotted ( Figure A7). The curve trends of Nehe, Gannan, Yi'an, Fuyu and Baiquan were the most consistent, as seen in A7a, and the curves in Figure A7b, A7c, and A7d showed roughly the same trend. (a) is category 1 in Figure 8, (b) is category 7, (c) is category 5, and (d) is category 8).

Appendix A.7. NAF Interannual Variation Curve
The interannual NAF curves of four similar regions were plotted ( Figure A7). The curve trends of Nehe, Gannan, Yi'an, Fuyu and Baiquan were the most consistent, as seen in A7a, and the curves in Figure A7b-d showed roughly the same trend. The interannual NAF curves of four similar regions were plotted ( Figure A7). The curve trends of Nehe, Gannan, Yi'an, Fuyu and Baiquan were the most consistent, as seen in A7a, and the curves in Figure A7b