Paddy Field Expansion and Aggregation Since the Mid-1950 s in a Cold Region and Its Possible Causes

Over the last six decades, paddy fields on the Sanjiang Plain have experienced rapid expansion and aggregation. In our study, land use and land cover changes related to paddy fields were studied based on information acquired from topographic maps and remote-sensing images. Paddy field expansion and aggregation were investigated through landscape indices and trajectory codes. Furthermore, the possible causes of paddy field expansion and aggregation were explored. Results indicated that such fields have increased by approximately 42,704 ha·y−1 over the past six decades. Approximately 98% of paddy fields in 2015 were converted from other land use types. In general, the gravity center moved 254.51 km toward the northeast, at a rate of approximately 4.17 km·y−1. The cohesion index increased from 96.8208 in 1954 to 99.5656 in 2015, and the aggregation index grew from 91.3533 in 1954 to 93.4448 in 2015, indicating the apparent aggregation of paddy fields on the Sanjiang Plain. Trajectory analyses showed that the transformations from marsh as well as from grassland to dry farmland and then into paddy fields were predominant. Climate warming provided a favorable environment for rice planting. Meanwhile, population growth, technological progress, and government policies drove paddy field expansion and aggregation during the study period.


Introduction
Accounting for approximately 15% of the world's cultivated land, paddy fields serve more than 50% of the population of the world [1].Simultaneously, paddy fields are a type of man-made wetland [2,3].Sharing both characteristics of wetlands and cultivated land, paddy fields influence food security, greenhouse gas (GHG) emissions, water scarcity, and virus transmission [4][5][6].Growing on flooded soil, rice is the most water-consuming crop.Its water-use efficiency has a great influence on water scarcity [7,8].Accounting for over one-tenth of the entire methane emission into the atmosphere, rice planting is an important source of methane flux that cannot be ignored [9,10].In addition, methane emission caused by paddy planting is one of the dominant sources of anthropogenic GHG emissions [6].As an important habitat for waterfowl, paddy fields have a relationship with the transmission of diseases such as avian influenza virus [11].Therefore, it is vital to monitor paddy field distribution information.Accompanying global warming, a large number of newly reclaimed paddy fields have emerged in middle-high latitudes [11,12].However, paddy field distribution maps are lacking in these regions, especially on a long-time scale.
As a critical driving force of global ecological and environmental changes, land use and land cover change (LULCC) is still a key environmental challenge of global concern [7,[13][14][15][16].LULCC manifests in atmospheric interactions, hydrological changes, and biodiversity losses [8, [17][18][19].Approximately 39-50% of terrestrial ecosystems have been affected by anthropogenic activities [20].Population growth, agricultural activities, government policies, and economic development all influence LULCC [12,15,21].Many efforts have been made to monitor the spatiotemporal patterns of LULCC, but these studies have mostly adopted a bitemporal detection method [21][22][23].Trajectory analysis seeks the idealized signatures of LULCC [24,25].Trajectory analysis can recover the history of LULCC as well as discover trends in LULCC at multiple times [24][25][26].Some studies have attempted to analyze LULCC using trajectory analysis [25,27,28].In our study, bitemporal detection and trajectory analysis are combined to trace the LULCC patterns related to paddy field expansion for each location.
The Sanjiang Plain is a floodplain that includes the largest marshlands in China.It is located in Northeastern China.However, due to the large-scale reclamation of marshlands, the Sanjiang Plain has been converted from a "great northern wilderness" to a "great northern granary" [28,29].Much attention has been paid to marshland loss and its effect in this area or part of this area [2,[29][30][31][32][33][34][35], but studies regarding landscape change and paddy fields are lacking.Song et al. analyzed the dynamic changes of cultivated land in this area, but they did not separate paddy fields from cultivated land [33].Considering the rapid increase in paddy field area in Northeastern China [36], it is necessary to monitor the spatial-temporal extent of paddy fields in detail in order to trace greenhouse gas (GHG) emissions and analyze food security as well as water security.Based on topographic maps and remote-sensing images, this study analyzed the expansion and aggregation characteristics of paddy fields on the Sanjiang Plain over the last six decades.Additionally, possible causes of paddy field expansion were analyzed.
In this study, our objectives were to (1) characterize landscape changes of paddy fields during the study period by related landscape indices in the study area and (2) explore the LULCC related to paddy fields using transition probabilities index and trajectory analysis from 1954 to 2015.

Study Area
This study focuses on the Sanjiang Plain (Figure 1), an alluvial plain formed by the Heilongjiang, Ussuri, and Songhua rivers.Covering an area of approximately 10.89 million ha, the elevation of the Sanjiang Plain ranges from 27 to 968 m, with an average value less than 200 m.The Sanjiang Plain has a temperate continental monsoon climate with four seasons.The frost-free days range from 120 to 140 days.The annual accumulated air temperature above 10 • C is approximately 2300-2500 • C. The average annual temperature varies between −21 and −18 • C in January, the coldest month, while the maximum average temperature ranges from 21 to 22 • C in July.The mean yearly rainfall varies from 500 to 650 mm, mainly concentrated in June and October [34,35].Rich water resources of the Sanjiang Plain include the Heilongjiang River, the Ussuri River, the Songhua River, the Naoli River, and the Xingkai Lake, providing facilities for agricultural development, especially paddy rice planting.

Data Preprocessing
Based on our previous studies [12,23,28,[36][37][38][39][40], we developed land-use-land-cover (LULC) maps for four different years : 1954, 1976, 1986, and 2000.Topographic maps with a scale of 1:100,000, Landsat Multispectral Scanner (MSS) images with a spatial resolution of 80 m, and Landsat Thematic Mapper (TM) images with a spatial resolution of 30 m were applied to map LULC patterns in 1954, 1976, and 1986, and 2000.Meanwhile, using 14 Landsat-8 Operational Land Imager (OLI) images with a spatial resolution of 30 m, we updated the LULC map in 2000 to 2015.Considering the coarse spatial resolution of MSS images (80 m), all the calculations were done at the resolution of 80 m.
The detailed methods to map LULC distribution for different years have been introduced in our previous articles [12,28,[36][37][38][39][40].The standard methods are visual interpretation and digitization of remote-sensing images after these images are geo-referenced and ortho-rectified.Ground control points were used to do geometric rectification of remote sensing (RS) images and the root mean squared (RMS) error of geometric rectification was no more than 1.5 pixels (120 m).In our study, the minimum mapping unit is larger than 23.04 ha, and the smallest edge is no less than 320 m.Additionally, we outlined the LULC by comparing images of different periods to assess the real LULCC.Uniform coordinates and projection, the Beijing 1954 Krasovsky Albers projection, was used to integrate satellite images in different periods.The interpretation accuracy was verified using a large number of photos by field survey, historical data such as aerial photos, and interviews with local people [36,41,42].In 2015, we also applied unmanned aerial vehicle (UAV) images to correct and check the interpretation of paddy fields (Figure 2).The battery-powered quadrocopter of our UAV can fly for 23 min at a time.We usually flew the UAV to an altitude of 200 m and obtained images with a resolution of 5-6 cm [43,44].The first level of classification is the same as Liu et al. [37], which includes cultivated land, forest land, grassland, water body, settlement, and unused land.We separated paddy fields from cultivated land and marsh from unused land to analyze paddy field changes and the related LULCC.The overall accuracy of the first level of classification was no less

Data Preprocessing
Based on our previous studies [12,23,28,[36][37][38][39][40], we developed land-use-land-cover (LULC) maps for four different years : 1954, 1976, 1986, and 2000.Topographic maps with a scale of 1:100,000, Landsat Multispectral Scanner (MSS) images with a spatial resolution of 80 m, and Landsat Thematic Mapper (TM) images with a spatial resolution of 30 m were applied to map LULC patterns in 1954, 1976, and 1986, and 2000.Meanwhile, using 14 Landsat-8 Operational Land Imager (OLI) images with a spatial resolution of 30 m, we updated the LULC map in 2000 to 2015.Considering the coarse spatial resolution of MSS images (80 m), all the calculations were done at the resolution of 80 m.
The detailed methods to map LULC distribution for different years have been introduced in our previous articles [12,28,[36][37][38][39][40].The standard methods are visual interpretation and digitization of remote-sensing images after these images are geo-referenced and ortho-rectified.Ground control points were used to do geometric rectification of remote sensing (RS) images and the root mean squared (RMS) error of geometric rectification was no more than 1.5 pixels (120 m).In our study, the minimum mapping unit is larger than 23.04 ha, and the smallest edge is no less than 320 m.Additionally, we outlined the LULC by comparing images of different periods to assess the real LULCC.Uniform coordinates and projection, the Beijing 1954 Krasovsky Albers projection, was used to integrate satellite images in different periods.The interpretation accuracy was verified using a large number of photos by field survey, historical data such as aerial photos, and interviews with local people [36,41,42].In 2015, we also applied unmanned aerial vehicle (UAV) images to correct and check the interpretation of paddy fields (Figure 2).The battery-powered quadrocopter of our UAV can fly for 23 min at a time.We usually flew the UAV to an altitude of 200 m and obtained images with a resolution of 5-6 cm [43,44].The first level of classification is the same as Liu et al. [37], which includes cultivated land, forest land, grassland, water body, settlement, and unused land.We separated paddy fields from cultivated land and marsh from unused land to analyze paddy field changes and the related LULCC.The overall accuracy of the first level of classification was no less than 94.3%, and that of the 25 LULC subclasses was no less than 91.2%.Additionally, the LULCC accuracies were also checked between two periods.For example, the mean interpretation accuracy of LULCC from 1995 to 2000 was 97.6% [36].
Remote Sens. 2018, 10, x FOR PEER REVIEW 4 of 19 than 94.3%, and that of the 25 LULC subclasses was no less than 91.2%.Additionally, the LULCC accuracies were also checked between two periods.For example, the mean interpretation accuracy of LULCC from 1995 to 2000 was 97.6% [36].

Annual Change Rate
The annual change area (CA, ha•y −1 ) and change rate (CR, %•y −1 ) was calculated using the following equations [21]: where Aa and Ab represent the area of the paddy fields at Times A and B, respectively, and T is the interval years.

Gravity Center
The land resource centroid index (Equations ( 3) and ( 4)) can reflect the spatial changes of land use [2]: where X and Y represent the gravity center coordinate for paddy fields in the entire study area; C and n are the area of the ith patch and the number of patches of paddy fields, respectively; X and Y represent the ith patch's gravity center coordinates.

Trajectory Codes
Trajectory code Yi (Equation ( 5)) was calculated as follows [43]: where n represents the number of time periods and n equals 5 in this study.(G1)i, (G2)i, and (Gn)i are the codes for different LULC categories in polygon i at different times.In our study, Codes 1 through 8 represent paddy, dry farmland, forest, grassland, water, settlement, marsh, and other unused land, respectively.Thus, for example, Code 44211 represents grassland → grassland → dry farmland → paddy→paddy.

Annual Change Rate
The annual change area (CA, ha•y −1 ) and change rate (CR, %•y −1 ) was calculated using the following equations [21]: where A a and A b represent the area of the paddy fields at Times A and B, respectively, and T is the interval years.

Gravity Center
The land resource centroid index (Equations ( 3) and ( 4)) can reflect the spatial changes of land use [2]: where X and Y represent the gravity center coordinate for paddy fields in the entire study area; Cti and n are the area of the ith patch and the number of patches of paddy fields, respectively; X and Y represent the ith patch's gravity center coordinates.

Trajectory Codes
Trajectory code Y i (Equation ( 5)) was calculated as follows [43]: where n represents the number of time periods and n equals 5 in this study.(G1) i , (G2) where P ij is the cumulative transition probabilities index; i and j represent LULC types; S t ij is the area converted from category i to category j between time t and t + 1; S T represents the entire area.

Spatial aggregation Analysis
Global Moran's I index [45][46][47] was used to analyze the aggregation characteristic of the entire Sanjiang Plain.It is calculated by the following formula: where x i and x j are the values of x in adjacent paring spatial units; x stands for the mean value of variable x; W ij stands for the adjacent weight; n is the number of space units.Global Moran's I coefficients are between −1 and 1.A positive value implies that the spatial distribution tends to aggregation status, while a negative value indicates a fragmentation trend.Local Moran's I index [47][48][49] was adopted to illustrate the aggregation of various spatial unites (Equation ( 8)): Due to the complexity of the landscape pattern of patches, we also chose more than one index including mean area (MA), largest patch index (LPI), patch density (PD), edge density (ED), patch cohesion (COHESION), division (DIVISION), splitting (SPLIT), and aggregation (AI) index to qualify the spatial expansion and aggregation features of paddy fields by Fragstats 4.2 software [50][51][52][53][54][55].More detailed descriptions of related indices can be found in the Fragstats 4.2 help documentation [56].In our study, MA, LPI, PD, and ED were chosen to illustrate the expansion characteristics, while the other four indices were applied to analyze the aggregation features (Supplementary Materials, Equations ( 1)-( 8)).

Change Rate
The annual change area/rate of paddy fields during the 1954-2015 period is shown in Table 1         This was in relation to the adoption of modern machine agriculture in the 1978-1985 period that was promoted by "Agricultural Modernization" activities.The spatial configuration changes in paddy fields are illustrated using four landscape indices for the study period (Figure 5).In general, all four indices showed an increasing trend over the last 60 years, indicating an expansion trend of paddy fields during the study period.The LPI and ED showed a continued increasing trend, while PD and MA increased initially and then decreased.The PD increased in the 1954-2000 period, with a rapid increase in the 1986-2000 period, and then showed a decreased trend from 2000 to 2015.During the 1986-2000 period, the PD increased quickly, while MA showed a declining trend, which meant that more patches of paddy fields were cultivated but their size was not as large as before.During the 2000-2015 period, the PD showed a declining trend, while the other three indices of LPI, ED, and MA showed an increasing trend, indicating that paddy fields were forming on a large scale.The spatial configuration changes in paddy fields are illustrated using four landscape indices for the study period (Figure 5).In general, all four indices showed an increasing trend over the last 60 years, indicating an expansion trend of paddy fields during the study period.The LPI and ED showed a continued increasing trend, while PD and MA increased initially and then decreased.The PD increased in the 1954-2000 period, with a rapid increase in the 1986-2000 period, and then showed a decreased trend from 2000 to 2015.During the 1986-2000 period, the PD increased quickly, while MA showed a declining trend, which meant that more patches of paddy fields were cultivated but their size was not as large as before.During the 2000-2015 period, the PD showed a declining trend, while the other three indices of LPI, ED, and MA showed an increasing trend, indicating that paddy fields were forming on a large scale.The remaining four indices (COHESION, DIVISION, SPLIT, and AI in Table 2) were chosen to illustrate the aggregation characteristics during the study period.A greater value of COHESION means more aggregation, while a greater value of DIVISION means less aggregation.In our study, the COHESION increased from 96.8208 in 1954 to 99.5656 in 2015, while the DIVISION decreased from 1 to 0.9986 during the study period.A greater SPLIT can be considered less aggregated, a greater AI (%) considered more aggregated.On the Sanjiang Plain, the SPLIT showed a continuous decreasing trend from 1954 to 2015, with a rapid decrease from 1954 to 1976.Including several fluctuations, the AI showed an increasing trend, with a peak in 1986.These changes indicated that paddy fields on the Sanjiang Plain have become increasingly aggregated over the past 60 years.The gravity centers of the paddy fields at different times were calculated using ArcGIS software (Figure 6).The gravity center of the paddy fields on the Sanjiang Plain transformed from The remaining four indices (COHESION, DIVISION, SPLIT, and AI in Table 2) were chosen to illustrate the aggregation characteristics during the study period.A greater value of COHESION means more aggregation, while a greater value of DIVISION means less aggregation.In our study, the COHESION increased from 96.8208 in 1954 to 99.5656 in 2015, while the DIVISION decreased from 1 to 0.9986 during the study period.A greater SPLIT can be considered less aggregated, a greater AI (%) considered more aggregated.On the Sanjiang Plain, the SPLIT showed a continuous decreasing trend from 1954 to 2015, with a rapid decrease from 1954 to 1976.Including several fluctuations, the AI showed an increasing trend, with a peak in 1986.These changes indicated that paddy fields on the Sanjiang Plain have become increasingly aggregated over the past 60 years.The gravity centers of the paddy fields at different times were calculated using ArcGIS software (Figure 6).The gravity center of the paddy fields on the Sanjiang Plain transformed from (45.88°N,130.83°E) in 1954 to (47.48°N,133.2°E) in 2015.As is clearly shown in Figure 6, the spatial pattern of paddy fields in the study area has moved northeastward since the mid-1950s.In general, the gravity center moved 254.51 km toward the northeast at a rate of approximately 4.17 km•y −1 over the last six decades (Table 3).The gravity centers of paddy fields migrated 31.38   In general, the gravity center moved 254.51 km toward the northeast at a rate of approximately 4.17 km•y −1 over the last six decades (Table 3).The gravity centers of paddy fields migrated 31.38 7).There were only two types of aggregation (H-H and H-L) from 1954 to 1986, which had a strong relationship with the relatively small paddy field area during this period.As paddy fields became more widely distributed, the other two types of aggregation (L-H and L-L) occurred in 2000.H-H accounted for the largest area and has shown a significant upward trend in the past six decades.H-H was only distributed in the south of the Sanjiang Plain in 1954, while it was widely distributed in the entire study area, especially the north of the Sanjiang Plain, in 2015, indicating that paddy field aggregation was enhanced especially in the north of the study area.L-L, which means their surrounding locations have low values, indicates that there are few paddy fields.By combining this information with LULC data, we found that the L-L area was mainly forestland.H-L and L-H have high negative local Moran's I values.Both H-L and L-H implies that the locations have values that are different from their neighbors.H-L indicates a low value in high-value surroundings, while L-H is the opposite.The H-L area increased from 1954 to 1976, which was mainly caused by small-scale and decentralized rice cultivation in this period.The H-L area moved northward from 1976 to 2000, and it mostly located in the south of the Sanjiang Plain in 2015.These changes indicated that paddy planting moved northward and the small-scale planting became aggregated in the northern part.The L-H area was relatively small and only showed an increasing trend from 2000 to 2015.

Trajectory Computing
To better illustrate the trajectories of paddy field changes, the initials of different LULC types were used to represent numeric codes.For example, "P" and "D" were used to replace "1" and "2", respectively.Figure 8 illustrates the trajectories of paddy field changes.Over the past 60 years, the unchanged paddy field, namely paddy field→paddy field→paddy field→paddy field→paddy field (PPPPP), occupied just 1.39% of the paddy field area in 2015, while the percentage of two-step changes was 53.07%and was the largest.Both one-step and three-step occupied approximately 20% of the paddy field area in 2015, while the percentage of four-step changes was 4.63% (Table 5).

Trajectory Computing
To better illustrate the trajectories of paddy field changes, the initials of different LULC types were used to represent numeric codes.For example, "P" and "D" were used to replace "1" and "2", respectively.Figure 8 illustrates the trajectories of paddy field changes.Over the past 60 years, the unchanged paddy field, namely paddy field→paddy field→paddy field→paddy field→paddy field (PPPPP), occupied just 1.39% of the paddy field area in 2015, while the percentage of two-step changes was 53.07%and was the largest.Both one-step and three-step occupied approximately 20% of the paddy field area in 2015, while the percentage of four-step changes was 4.63% (Table 5).Seventeen categories were included in the one-step changes (Supplementary Materials, Table S1), which mostly included conversion from marsh to paddy field (such as marsh→paddy field→paddy field→paddy field→paddy field(MPPPP), marsh→marsh→marsh→paddy field→paddy field, (MMMPP), and marsh→marsh→paddy field→paddy field→paddy field (MMPPP)) and that from dry farmland to paddy (such as dry farmland→paddy field→paddy field→paddy field→paddy field (DPPPP) and dry farmland→dry farmland→paddy field→paddy field→paddy field (DDPPP)).In addition, MPPPP (marsh→paddy field) accounted for the largest proportion of these one-step changes, followed by MMMFF (marsh→forest, 2.11%), DPPPP (dry farmland→paddy field), and DDDDP (dry farmland→paddy field).Seventeen categories were included in the one-step changes (Supplementary Materials, Table S1), which mostly included conversion from marsh to paddy field (such as marsh→paddy field→paddy field→paddy field→paddy field(MPPPP), marsh→marsh→marsh→paddy field→paddy field, (MMMPP), and marsh→marsh→paddy field→paddy field→paddy field (MMPPP)) and that from dry farmland to paddy (such as dry farmland→paddy field→paddy field→paddy field→paddy field (DPPPP) and dry farmland→dry farmland→paddy field→paddy field→paddy field (DDPPP)).In addition, MPPPP (marsh→paddy field) accounted for the largest proportion of these one-step changes, followed by MMMFF (marsh→forest, 2.11%), DPPPP (dry farmland→paddy field), and DDDDP (dry farmland→paddy field).Two-step changes (84 types, Table S2 in Supplementary Materials) during the study period mainly included MDDPP (10.16%),MDDDP (8.7%), MMDDP (6.36%), MMDPP (6.05%), GDDPP (3.13%), and GDDDP (3.03%), indicating that the main two-step changes of paddy fields were marsh→dry farmland→paddy field and grassland→dry farmland→paddy field.There were numerous categories of three-step as well as four-step changes, all accounting for only a small proportion (Tables S3 and S4 in Supplementary Materials).The conversion from dry farmland reclaimed from marsh and grassland to paddy occupied a larger proportion of these trajectories.Furthermore, MPDPP (marsh→paddy field→dry farmland→paddy field→paddy field) and MDPDP (marsh→dry farmland→paddy field→dry farmland→paddy field) occupied the greatest percentage of three-step and four-step changes, respectively.These three-step and four-step trajectories illustrated that conversions between paddy field and dry farmland took place relatively frequently, particularly the transformation from dry farmland to paddy field.The conversion from rice planting to dry farmland was mostly a result of limited water resources, while the conversion from dry farmland to rice planting was due to the promotion of related polices and higher economic benefits of rice.For example, "promoting the conversion from dry farmland to paddy" policy and "controlling flooding by rice planting" policy in the northeastern study area greatly promoted the expansion of paddy fields.S2 in Supplementary Materials) during the study period mainly included MDDPP (10.16%),MDDDP (8.7%), MMDDP (6.36%), MMDPP (6.05%), GDDPP (3.13%), and GDDDP (3.03%), indicating that the main two-step changes of paddy fields were marsh→dry farmland→paddy field and grassland→dry farmland→paddy field.There were numerous categories of three-step as well as four-step changes, all accounting for only a small proportion (Tables S3 and S4 in Supplementary Materials).The conversion from dry farmland reclaimed from marsh and grassland to paddy occupied a larger proportion of these trajectories.Furthermore, MPDPP (marsh→paddy field→dry farmland→paddy field→paddy field) and MDPDP (marsh→dry farmland→paddy field→dry farmland→paddy field) occupied the greatest percentage of three-step and four-step changes, respectively.These three-step and four-step trajectories illustrated that conversions between paddy field and dry farmland took place relatively frequently, particularly the transformation from dry farmland to paddy field.The conversion from rice planting to dry farmland was mostly a result of limited water resources, while the conversion from dry farmland to rice planting was due to the promotion of related polices and higher economic benefits of rice.For example, "promoting the conversion from dry farmland to paddy" policy and "controlling flooding by rice planting" policy in the northeastern study area greatly promoted the expansion of paddy fields.

Cumulative Transition Probability (P ij )
LULCC is a mutual transformation process.Results indicated that the values of P ij between paddy field and other LULC types were 9.1% (from paddy field to other LULC types) and 32.6% (from other LULC types to paddy field) of the total area, respectively.The transformations between paddy field and other LULC types are demonstrated in Figure 9.The dominant contributor to paddy field expansion was dry farmland (69.7% of the conversion from other LULC types to paddy field), followed by marsh (18.5%), grassland (8.5%), and forest (1.9%).Dry farmland (82.3%), marsh (7.9%), grassland (3.5%), and built-up land (2.3) made the greatest contributions to paddy field decrease.The mutual transformation between paddy field and dry farmland was the main LULCC type of paddy field conversions.The transformation from dry farmland to paddy field was mostly promoted by the exploitation of groundwater for irrigation, with some policies including promoting the conversion from dry farmland to paddy field and controlling flooding by rice planting, while the conversion from paddy field to dry farmland was mostly due to limited water resources.
LULCC is a mutual transformation process.Results indicated that the values of Pij between paddy field and other LULC types were 9.1% (from paddy field to other LULC types) and 32.6% (from other LULC types to paddy field) of the total area, respectively.The transformations between paddy field and other LULC types are demonstrated in Figure 9.The dominant contributor to paddy field expansion was dry farmland (69.7% of the conversion from other LULC types to paddy field), followed by marsh (18.5%), grassland (8.5%), and forest (1.9%).Dry farmland (82.3%), marsh (7.9%), grassland (3.5%), and built-up land (2.3) made the greatest contributions to paddy field decrease.The mutual transformation between paddy field and dry farmland was the main LULCC type of paddy field conversions.The transformation from dry farmland to paddy field was mostly promoted by the exploitation of groundwater for irrigation, with some policies including promoting the conversion from dry farmland to paddy field and controlling flooding by rice planting, while the conversion from paddy field to dry farmland was mostly due to limited water resources.

Uncertainty Analysis
Integrating remote-sensing images and historical data is one of the most economically feasible methods to obtain long time series of LULC maps.One question that needs attention when illustrating LULCC is the interpretation accuracy for LULC categories and accuracy for change detection.To validate the accuracies of our LULC data, a large amount of photos that is taken by camera, UAV images, historical data such as Statistical Yearbook, and field site data, as well as interviews with local residents were used to check interpretation accuracy.Visual interpretation is the main method of obtaining LULCC data in our study, which may be labor-intensive and timeconsuming, but it guarantees relatively high accuracy.One problem that we should pay attention to is the amplification of individual classification errors in the change detection process.To reduce these errors, we outline LULC by comparing images of different periods to assess the real LULCC.Additionally, the accuracy for change detection was also validated by historical data and a field survey.In general, our LULCC maps can meet the requirement of user accuracy despite some uncertainties.

Climatic Warming Favorable for Paddy Field Expansion
Continuously rising temperatures over the past few decades in Northeastern China have been reported in previous studies [34,57].In our study, we calculated the average temperature change over the past six decades on the Sanjiang Plain (Figure 9).The annual temperature data was produced by

Uncertainty Analysis
Integrating remote-sensing images and historical data is one of the most economically feasible methods to obtain long time series of LULC maps.One question that needs attention when illustrating LULCC is the interpretation accuracy for LULC categories and accuracy for change detection.To validate the accuracies of our LULC data, a large amount of photos that is taken by camera, UAV images, historical data such as Statistical Yearbook, and field site data, as well as interviews with local residents were used to check interpretation accuracy.Visual interpretation is the main method of obtaining LULCC data in our study, which may be labor-intensive and time-consuming, but it guarantees relatively high accuracy.One problem that we should pay attention to is the amplification of individual classification errors in the change detection process.To reduce these errors, we outline LULC by comparing images of different periods to assess the real LULCC.Additionally, the accuracy for change detection was also validated by historical data and a field survey.In general, our LULCC maps can meet the requirement of user accuracy despite some uncertainties.

Climatic Warming Favorable for Paddy Field Expansion
Continuously rising temperatures over the past few decades in Northeastern China have been reported in previous studies [34,57].In our study, we calculated the average temperature change over the past six decades on the Sanjiang Plain (Figure 9).The annual temperature data was produced by data from seven meteorological stations (China Meteorological Administration, CMA) across the study area since the 1950s.
Results (Figure 10) indicated that yearly average temperature showed an increasing trend despite several fluctuations.The regression relationship between yearly average temperature and year can be described by the following formula: y = 0.021x + 2.6139 where y represents the annual average temperature, and x represents a year.
study area since the 1950s.Results (Figure 10) indicated that yearly average temperature showed an increasing trend despite several fluctuations.The regression relationship between yearly average temperature and year can be described by the following formula: where y represents the annual average temperature, and x represents a year.The annual average temperature has increased at a rate of 0.21 °C/10 y on the Sanjiang Plain.By the Kriging interpolation method, we obtained the spatial pattern of average temperature in two periods (Figure 11). Figure 11 indicated that the average temperature showed an obvious upward trend in the past six decades.For example, the orange color, which means the average temperature from 3.1 to 3.7 °C, moved northward obviously.In cold regions, agricultural activities can be limited by low temperatures.For example, rice cannot be grown in regions with lower yearly average temperature [2].On the one hand, climatic warming has promoted the conversion of marsh into cropland (paddy field and dry farmland), leading to paddy field expansion to some extent.On the other hand, climatic warming has provided a favorable environment for large-scale paddy planting in the northern part of the Sanjiang Plain, leading to the aggregation of paddy fields.Strong wind is also one of the limiting factors affecting agricultural production.Wind speed has shown an obvious decline trend in the past several decades in Northeastern China [58], which has benefited paddy growth.Rising temperature and declining wind speed have contributed to the conversion from wetland and grassland to cultivated land.The annual average temperature has increased at a rate of 0.21 • C/10 y on the Sanjiang Plain.By the Kriging interpolation method, we obtained the spatial pattern of average temperature in two periods (Figure 11). Figure 11 indicated that the average temperature showed an obvious upward trend in the past six decades.For example, the orange color, which means the average temperature from 3.1 to 3.7 • C, moved northward obviously.In cold regions, agricultural activities can be limited by low temperatures.For example, rice cannot be grown in regions with lower yearly average temperature [2].On the one hand, climatic warming has promoted the conversion of marsh into cropland (paddy field and dry farmland), leading to paddy field expansion to some extent.On the other hand, climatic warming has provided a favorable environment for large-scale paddy planting in the northern part of the Sanjiang Plain, leading to the aggregation of paddy fields.Strong wind is also one of the limiting factors affecting agricultural production.Wind speed has shown an obvious decline trend in the past several decades in Northeastern China [58], which has benefited paddy growth.Rising temperature and declining wind speed have contributed to the conversion from wetland and grassland to cultivated land.

Population Growth and Technological Progress
Previous studies have indicated that LULCC has a close relationship with human activities [14,18,24,59,60] including population growth and intensive reclamation.Population growth and technological progress have provided convenience for paddy field expansion.From 1956 to 1960, approximately 8.2 × 10 4 veterans migrated to the Sanjiang Plain to reclaim marsh as cultivated land [61].Then approximately 4.5 × 10 4 educated young people moved into the study area to work in agriculture from 1970 to 1972 in response to the call for educated youth to live in the countryside to accept poor farmers' re-education [62].This population growth led to large-scale reclamation as well as an increase in paddy field area.
During the early 1980s, "agricultural modernization" activities in China brought advanced agricultural machinery, leading to large-scale planting of crops including rice.Water resources were also an important limiting factor for rice planting.Rice needs abundant water to grow [2].Some water conservancy projects on the Sanjiang Plain provide convenient conditions for rice planting [30].As shown in Figure 12 (the photographs were taken during our field survey), a dam was built to irrigate paddy fields.Additionally, technological innovations such as the seeding of cold-tolerant rice, the use of fertilizer, and advances in irrigation have provided convenient conditions for paddy field expansion.Modern machine and advances in irrigation also provided convenient conditions for large-scale planting, leading to paddy field aggregation.Figure 10 shows the irrigation project and advanced germination and soaking facilities of the Suibin state-farm, located in the north of the Sanjiang Plain.The irrigation project introduced water from the Amur River for irrigation.

Role of Government Policies
On the Sanjiang Plain, government policy shifts have also played an essential role in the progress of paddy field expansion and aggregation.The First Five-Year Plan emphasized "food first" in agricultural production, promoting an increase in cultivated land from 1953 to 1957.In 1964, the government encouraged educated youth to participate in rural socialist construction, causing the large-scale movement of educated young people to the Sanjiang Plain for agricultural work [61,62].During this period, marsh was converted into farmland on a large-scale, promoting paddy field expansion.The introduction of modern machinery by the "Agricultural Modernization" [2,43] policy from 1978 to 1985 promoted large-scale reclamation as well as paddy field aggregation in the Sanjiang Plain.
In 1992, the planned economy was replaced by a market economy in China [60].As a result of the higher income of rice planting, as compared to that of dry farmland crops, farmers paid more attention to the conversion of dry farmland to paddy fields, leading to paddy field expansion in the Sanjiang Plain.In addition, large-scale dry farmland was converted to paddy fields as a result of the

Role of Government Policies
On the Sanjiang Plain, government policy shifts have also played an essential role in the progress of paddy field expansion and aggregation.The First Five-Year Plan emphasized "food first" in agricultural production, promoting an increase in cultivated land from 1953 to 1957.In 1964, the government encouraged educated youth to participate in rural socialist construction, causing the large-scale movement of educated young people to the Sanjiang Plain for agricultural work [61,62].During this period, marsh was converted into farmland on a large-scale, promoting paddy field expansion.
The introduction of modern machinery by the "Agricultural Modernization" [2,43] policy from 1978 to 1985 promoted large-scale reclamation as well as paddy field aggregation in the Sanjiang Plain.
In 1992, the planned economy was replaced by a market economy in China [60].As a result of the higher income of rice planting, as compared to that of dry farmland crops, farmers paid more attention to the conversion of dry farmland to paddy fields, leading to paddy field expansion in the Sanjiang Plain.In addition, large-scale dry farmland was converted to paddy fields as a result of the "promoting the conversion from dry farmland to paddy field" policy from 1992 to 1995, especially in the state farms of the Sanjiang Plain, which promoted paddy field expansion and aggregation.As an important commodity grain base, the implementation of the "Layout Planning for the Advantageous Regions of Rice (2008-2015)" policy issued by the Ministry of Agriculture in 2009 also greatly promoted the expansion of paddy field area on the Sanjiang Plain.The "Building High-Standard Basic Farmland" policy in the Twelfth Five-Year Plan (2011-2015) promoted agricultural modernization on a large scale, leading to paddy field aggregation in the study area.

Future Study
The planting of paddy fields on the Sanjiang Plain plays an important role in the development of a regional agricultural economy.However, the expansion of paddy rice fields will inevitably influence biogeochemical processes, water resources utilization, and grain yield.The rapid expansion and aggregation of paddy fields can lead to serious over-exploitation of groundwater, increasing the contradiction between supply and demand of water resources as well as the danger of groundwater depletion to some extent.Therefore, during paddy field development, based on a comprehensive evaluation of regional water resources carrying capacity, one should heed an early scientific warning regarding the utilization of water and soil resources and their coupling and coordinating relations.One should scientifically and reasonably implement a matching of water and land resources as well as healthy and sustainable development of the social economy and ecological environment.Paddy fields are also an important source of methane emissions.Zhang et al. found that methane emissions from paddy fields on a global scale in the 2000s ranged from 18.3 ± 0.1 Tg under intermittent irrigation to 38.8 ± 1.0 Tg under continuous flooding [6].Therefore, the effect of paddy field expansion on the water cycle and the regional agricultural climate should receive substantial attention.

Conclusions
Data analyses of topographic and remote-sensing images over the last 60 years have indicated that paddy field area has increased from 59,325 ha in 1954 to 2,604,946 ha in 2015.The increased areas were mainly converted from dry farmland (69.7%), marsh (18.5%), and grassland (8.5%).Trajectory analyses showed that the transformations from marsh and grassland into dry farmland, and subsequently into paddy fields, were the main categories to be considered.The lost paddy field areas were notably smaller than those of increased area, and the lost paddy field areas were mostly transformed to dry farmland.Paddy field expansion and aggregation has been driven by population growth, technological progress, climate change, and governmental policies.Paddy field expansion will affect food security, GHG emissions, and water security, among others.Future studies should pay increasing attention to these issues.

Figure 1 .
Figure 1.Study area location and elevation.

Figure 1 .
Figure 1.Study area location and elevation.

Figure 2 .
Figure 2. Field survey by unmanned aerial vehicles (UAVs).(left) UAV images; (right) the UAV that used in our field survey.

Figure 2 .
Figure 2. Field survey by unmanned aerial vehicles (UAVs).(left) UAV images; (right) the UAV that used in our field survey.

where s 2
represents the variance of x.A high positive I means that the locations have values similar to those of their surrounding locations.High-high clusters (H-H) mean high values in high neighbors while low-low clusters (L-L) are the opposite.High negative I values imply that values of these locations are quite different from values of their neighbors.High negative I values contain a high-low outlier (H-L, a high value surrounded by low values) and a low-high outlier (L-H, a low value surrounded by high values).

3. 1 . 1 .
Percentages of Paddy FieldsThe area percentage of paddy fields during the 1954-2015 period on the Sanjiang Plain is illustrated in Figure3.Paddy fields have expanded from 59,325 ha in 1954 to 2,604,946 ha in 2015, increasing by 43 times (approximately 2.5 million ha).The area percentage of paddy fields dramatically increased from 0.54 in 1954 to 23.92% in 2015.As is shown in Figure3, the area of paddy fields has continued to grow since the mid-1950s, especially after 1986.Remote Sens. 2018, 10, x FOR PEER REVIEW 6 of 19 dramatically increased from 0.54 in 1954 to 23.92% in 2015.As is shown in Figure 3, the area of paddy fields has continued to grow since the mid-1950s, especially after 1986.

Figure 3 .
Figure 3.The percentages of paddy from 1954 to 2015.

Figure 4
Figure 4 demonstrates the spatial distribution of paddy fields in 1954, 1986, and 2015 on the Sanjiang Plain, showing an obvious trend of area expansion.In 1954, paddy fields were mainly distributed in the southern part of the Sanjiang Plain and were of small area.In 1986, paddy fields were cultivated northward, with large-scale paddy fields on the central plain.Paddy planting continued to move northward, with a large distribution on the northern part of the Sanjiang Plain in 2015.

Figure 4 .
Figure 4. Spatial distribution of paddy fields in different years.
. The paddy field area increased by 15,808.86ha•y −1 with an annual change rate of 26.65%•y −1 , from 1954 to 1976.Since 1976, the annual change area of paddy fields has continued to increase, reaching a maximum of 78,945.39ha•y −1 in the 2000-2015 period.The great growth in the annual change area in the 1986-2000 period, a number approximately 3.5 times that of the period from 1976 to 1986, is

Figure 3 .
Figure 3.The percentages of paddy from 1954 to 2015.

Figure 4
Figure 4 demonstrates the spatial distribution of paddy fields in 1954, 1986, and 2015 on the Sanjiang Plain, showing an obvious trend of area expansion.In 1954, paddy fields were mainly distributed in the southern part of the Sanjiang Plain and were of small area.In 1986, paddy fields were cultivated northward, with large-scale paddy fields on the central plain.Paddy planting continued to move northward, with a large distribution on the northern part of the Sanjiang Plain in 2015.

Figure 3 .
Figure 3.The percentages of paddy from 1954 to 2015.

Figure 4
Figure 4 demonstrates the spatial distribution of paddy fields in 1954, 1986, and 2015 on the Sanjiang Plain, showing an obvious trend of area expansion.In 1954, paddy fields were mainly distributed in the southern part of the Sanjiang Plain and were of small area.In 1986, paddy fields were cultivated northward, with large-scale paddy fields on the central plain.Paddy planting continued to move northward, with a large distribution on the northern part of the Sanjiang Plain in 2015.

Figure 4 .
Figure 4. Spatial distribution of paddy fields in different years.
3.1.2.Change RateThe annual change area/rate of paddy fields during the 1954-2015 period is shown in Table1.The paddy field area increased by 15,808.86ha•y −1 with an annual change rate of 26.65%•y −1 , from 1954 to 1976.Since 1976, the annual change area of paddy fields has continued to increase, reaching a maximum of 78,945.39ha•y −1 in the 2000-2015 period.The great growth in the annual change area in the 1986-2000 period, a number approximately 3.5 times that of the period from 1976 to 1986, is

Figure 4 .
Figure 4. Spatial distribution of paddy fields in different years.
3.1.2.Change RateThe annual change area/rate of paddy fields during the 1954-2015 period is shown in Table1.The paddy field area increased by 15,808.86ha•y −1 with an annual change rate of 26.65%•y −1 , from 1954 to 1976.Since 1976, the annual change area of paddy fields has continued to increase, reaching a maximum of 78,945.39ha•y −1 in the 2000-2015 period.The great growth in the annual change area in the 1986-2000 period, a number approximately 3.5 times that of the period from 1976 to 1986, is notable.

Figure 5 .
Figure 5. Spatial configuration in different time intervals.(a) PD index; (b) LPI index; (c) ED index and (d) MA index.

Figure 5 .
Figure 5. Spatial configuration in different time intervals.(a) PD index; (b) LPI index; (c) ED index and (d) MA index.
(45.88 • N,130.83• E) in 1954 to (47.48 • N,133.2 • E) in 2015.As is clearly shown in Figure 6, the spatial pattern of paddy fields in the study area has moved northeastward since the mid-1950s.Remote Sens. 2018, 10, x FOR PEER REVIEW 8 of 19

Figure 6 .
Figure 6.The gravity centers of paddy fields during different times.
km northwestward and 75.43 km northeastward from 1954 to 1976 and from 1976 to 1986, respectively.

Figure 6 .
Figure 6.The gravity centers of paddy fields during different times.

3. 1
.5.Spatial Aggregation Global Moran's I coefficients for paddy fields in the Sanjiang Plain from 1954 to 2015 are shown in Table 4.In general, Global Moran's I coefficients increased from 0.0583 in 1954 to 0.1833 in 2015, indicating that paddy fields have become more aggregated over the past 60 years.Specifically, Global Moran's I displayed a slight downward trend during the 1954-1976 and 1976-1986 periods, while it showed an upward tendency during the 1986-2000 and 2000-2015 periods.From 1954 to 1986, more patches were reclaimed as paddy fields, but their distribution was less aggregated than before, leading to the decline in Global Moran's I.The increase in Global Moran's I from 1986 to 2015 was caused by large-scale planting of paddy rice.

Figure 8 .
Figure 8.(a) Trajectories with different steps and (b) main trajectories of paddy field changes.P represents paddy field, D represents dry farmland, F represents forest, G represents grassland, and M represents marsh, respectively.

Figure 8 .
Figure 8.(a) Trajectories with different steps and (b) main trajectories of paddy field changes.P represents paddy field, D represents dry farmland, F represents forest, G represents grassland, and M represents marsh, respectively.

Figure 9 .
Figure 9. Area contribution of transformation between paddy field and other LULC types from 1954 to 2015.(Left) Transformation from paddy field to other LULC types; (right) transformation from other LULC types to paddy field.

Figure 9 .
Figure 9. Area contribution of transformation between paddy field and other LULC types from 1954 to 2015.(Left) Transformation from paddy field to other LULC types; (right) transformation from other LULC types to paddy field.

Figure 11 .Figure 11 .
Figure 11.(a) Average temperature from 1954 to 1964 and (b) average temperature from 2004 to 2014.4.3.Population Growth and Technological Progress Previous studies have indicated that LULCC has a close relationship with human activities [14,18,24,59,60] including population growth and intensive reclamation.Population growth and technological progress have provided convenience for paddy field expansion.From 1956 to 1960, approximately 8.2 × 10 4 veterans migrated to the Sanjiang Plain to reclaim marsh as cultivated land Remote Sens. 2018, 10, x FOR PEER REVIEW 15 of 19

Table 1 .
Rate of change of paddy fields during different time periods.

Table 1 .
Remote Sens. 2018, 10, x FOR PEER REVIEW 7 of 19 notable.This was in relation to the adoption of modern machine agriculture in the 1978-1985 period that was promoted by "Agricultural Modernization" activities.Rate of change of paddy fields during different time periods.
km northwestward and 75.43 km northeastward from 1954 to 1976 and from 1976 to 1986, respectively.The migration rate was 1.43 km•y −1 from 1954 to 1986, while it increased to 7.54 km•y −1 during the period of 1976-1986.From 1986 to 2000, it migrated 33.07 km northeastward at a rate of 2.36 km•y −1 , then again 146.19 km northeastward at a rate of 9.75 km•y −1 during the 2000-2015 period.

Table 3 .
The migration of the gravity center of paddy fields in different time periods.

Table 4 .
Global Moran's I for paddy fields at different times.

Table 5 .
The percentage of paddy field change trajectories with different steps.

Table 5 .
The percentage of paddy field change trajectories with different steps.