Evaluation of riparian condition of Songhua River by integration of remote sensing and field measurements

Riparian zone is crucial to the health of streams and their surrounding environment. Evaluation of riparian condition is essential to achieve and maintain good stream health, as well as to sustain ecological functions that riparian areas provide. This manuscript is aimed to evaluate riparian conditions of Songhua River, the fifth longest river in China, using physical structural integrality (PSI) values derived from remote sensing and validated by field measurements. The variation and clusters of PSI values were discriminated by the spatial statistics to quantify variation of riparian condition in each measurement section. Evaluation results derived from 13 measurement sections indicated that over 60% of the riparian zones have been disturbed by human activities. Analysis of land use patterns of riparian zone in the cold and hot spots found that land-use patterns had an important effect on riparian condition. The build-up and farmland areas had been the main human disturbances to the riparian condition, which were increased from 1976 to 2013. The low-low clusters (low PSI values with low neighbors) of PSI values can be implemented to identify the vulnerability of the riparian zone.

width, riparian vegetation and bank stability were important and feasible to extract from remote sensing data for condition assessment of riparian zones. However, evaluate riparian conditions using those indicators still remain to be examined.
This study evaluated the riparian condition of the Songhua River across the Northeastern Plain of China, using a series of indicators developed from 2,081 basic evaluation units (BEUs). The specific objectives of this paper are to: (1) demonstrate the feasibility of the multi-metric approach through comparisons with field measurements; (2) discriminate the variation and clusters of riparian condition to identify the vulnerability and stability of the riparian zones; and (3) explore the change of landscape of the riparian zone using the land use data from 1976 to 2013 to understand the effects of human impacts on riparian conditions.

Materials and Methods
Study area. Songhua River is the largest tributary of the Heilongjiang (Amur River) in the Northeast China.
With the length of 1,897 kilometers long and the drainage area of 545,600 square kilometers, Songhua is ranked the fifth longest river in China. It is originated from Changbai Mountain, passing through the Northeastern Plain, and injected into Amur River at about 48° north latitude. This study focused on the 1,679 kilometers riparian zone of Songhua River started from the Fengman hydrologic dam passing through major cities such as Jilin, Harbin, and Jiamusi, and ended at the Amur River, the border river between China and Russia (Fig. 1). The elevation of the river basin is between 54 to 2,735 meters above mean sea level. The basin has temperate continental climate with four clearly distinct seasons. The mean annual rainfall in the area is about 550-800 millimeter. The river reaches its maximum flow in summer and is frozen from November to April next year.
Songhua River basin is a major agricultural center and commodity grain production base of China. Large scale operations of state farms were established after the establishment of the People's Republic of China in 1949 and land reclamation projects begun since. Songhua River is an important transportation artery for agricultural products of the region. The river flows through major port cities and serves as the source of drinking water for millions of people. The ecological functions and integrality of riparian zones of the Songhua River have experienced much degradation under the influence of anthropogenic activities such as overgrazing, construction of transportation infrastructure, urban sprawl, sand mining, tourism development, reclaimed wetland and other human activities.
Land cover data and ancillary data. This study adopted the land-use data at 1:100,000 scale ( Fig. 2) for 1976,1986,1995,2000,2005 and 2013 from the Resources and Environment Science Data Center, Chinese Academy of Sciences. The land-use datasets included information of forest, grassland, farmland, wetland, water body, barren land and build-up (Table 1). The land use dataset of 2013 and 2015 derived from updating the land use map of 2005 with interpretation keys from field measurement sites using Landsat 8 Operational Land Imager (OLI) images acquired in 2013 and 2015, respectively. Table 2 describes selected multi-spectral OLI images that covered the study area with 30 m spatial resolution.
Other datasets adopted including 1:500,000 geomorphic map (Fig. 3) and 1:100,000 topographic map developed by the Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, and Shuttle Radar Topographic Mission (SRTM) generated Digital Elevation Model (DEM) data at 30 m spatial resolution. Referencing the 1:100,000 topographic map and 1:500,000 geomorphic map, the elevation range of the riparian zone was between 115 and 370 meters.
Field data collection. The studied riparian zone was divided into 13 measurement sections according to the Chinese Ministry of Water Resources technical protocol 11 . Within each section three or four measurement sites were chosen, and each site consisted of three transections (Fig. 4). Each transection was located in an area with a consistent management regime and with similar vegetation and stream characteristics. The entire riparian zone consisted of 45 measurement sites (Fig. 1) and 135 transections. Each transection laid two 30 m width × 50 m length measurement area on both side of banks. The positions of measurement sites and transections were guided by Global Positioning System (GPS). The sub-indicators of field-based PSI include bank slope and slope length of river bank, Canopy cover, bank sediment size (clay, sand, gravel, cobbles or bedrock), bank erosion and anthropogenic activities ( Table 3).
The bank slope and slope length were directly measured by a laser rangefinder, and used to calculate slope height. Canopy cover was measured by 1 × 1 m quadrat within the 30 × 50 m evaluation area. The 1 × 1 m vegetation quadrat was photographed with a camera mounted on a vertical shooting platform. These photos were subsequently processed to calculate vegetation coverage in every quadrat with the binary segmentation method by image-processing, and then calculated all Canopy cover in the evaluation area. Anthropogenic activities, bank sediment size, erosion and river connectivity were recorded by field photos. All field measurements were executed during the normal river flow period, in Septembers of 2013 and 2015, respectively.
Description of BEUs. The BEUs is defined as the basic evaluation units with various sizes in the measurement sections of the riparian zone (Fig. 1). To obtain the BEUs, the riparian zone was retrieved as a single featured GIS polygon using the 1:500,000 geomorphic map (Fig. 3). The river centerline was extracted from the river mouth to the source for the main channel based on the 1:100,000 topographic map and the DEM. The 600 m length was adopted following the Chinese Ministry of Water Resources technical protocol 11 as a splitting criterion to divide the river centerline. Then, the polygon of the riparian zone was separated using splitting lines perpendicular to the river centerline. This process generated 2,081 BEUs in the study area. The minimum, average and maximum size of the BEUs were 0.28 km 2 , 2.89 km 2 and 19.81 km 2 , respectively. All indicators and sub-indicators were calculated and scored based on the grid of BEUs.
Scientific RepoRts | 7: 2565 | DOI:10.1038/s41598-017-02772-3 PSI Indicators derived from remote sensing data. Adapting field-based indicators to remote sensing has been done through a combination of direct conversion and substitution (Table 4). Bank slope was calculated by DEM. Vegetation percent cover (VPC) and human disturbance (RD) were calculated by the Equation (1) and (2) using the updated 1:100,000 land use data. The width of riparian zone is affected synthetically by bank sediment size, flow erosion and deposition 23 . Previous studies utilized water level width and riparian zone area as sub-indicators to assess bank stability [24][25][26] . In addition, these sub-indicators could be accurately calculated by remote sensing data. In this paper, the area of riparian zone was calculated based on the 1:500,000 geomorphic map (Fig. 3) and 1:100,000 topographic map. The water level width of transection was extracted using OLI images. The two sub-indicators were normalized for converting the value to 0-100 (Equation (3)). Natural wetland conservation (NWC) was calculated using the 1976, 2013 and 2015 1:100,000 land-use data (Equation (4)). Vegetation canopy coverage (BVC) was calculated by Equation (5)  Where, Area BEU is the area of a BEU; Area i is the area of each vegetation type in the BEU; i is the vegetation types, including arboreal forest land, open forest land and shrub land, high-, mid-and low-cover grassland.     Where, Area is the area of natural wetland with closely hydraulic connection with river; C is the year of evaluating PSI, and R is the historical year; Ns is the number of wetland.
Where, Normalized Difference Vegetation Index (NDVI) was calculated using the September OLI images; NDVI soil and NDVI veg represent the NDVI in the vegetated and non-vegetated area, respectively. In this study, NDVI soil was the NDVI of 5% cumulative frequency, and NDVI veg was the NDVI of 95% cumulative frequency in the study area 28,29 .
Calculated PSI of the riparian zone. Sub-indicators and indicators of PSI based on remote sensing method were scored by modification of the score sheets of Chinese Ministry of Water Resources technical protocol (Table A, Supplement Materials). Each indicator was given a score between 0 and 100, with higher number implying better condition. The indicators were calculated with associated sub-indicators that were scored multiplying by their corresponding weights. The indicator of bank stability was calculated by the Equation (6). The PSI of riparian zone was calculated by the Equation (7).

Statistical validation of evaluation results. To validate evaluation results, an independent samples t-test
and F-test were performed to examine the statistically significant difference in PSI scores between field-evaluation and remote sensing observation. A test of normality was first performed using the Kolmogorov-Smirnov statistic and Shapiro-Wilk statistic to examine the possibility that the PSI scores derived from 45 measurement sites conform to a normal distribution. Under the circumstance of no significant difference, the study sought to understand the similarity between the remote sensing and field derived PSI, which led to the use of coefficients used for understanding model outputs. We used the common hydrological modeling coefficient of Nashe and Sutcliffe (NS) to evaluate similarity of riparian PSI between the field and remote sensing observation. NS is a method for quantitatively analyzing modeled outputs as they compare to observed values 30,31 . The method was appropriate as we compared the observed field data and the modeled remote sensing approach.
Spatial statistic and analysis of riparian PSI. The Getis-Ord General Gi and the Anselin local Moran's I statistics were calculated to analyze the spatial distribution of riparian PSI. The General Gi indicator was used to identify concentrations of high or low PSI values, while Moran's I indicator was applied for calculating spatial autocorrelation of PSI between each BEU and its neighboring BEUs, and identifying spatial outliers of PSI value [32][33][34] . A search radius or threshold distance is essential for the two local statistics. In this study, the 'Incremental Spatial Autocorrelation' Toolbox in ArcGIS was used to determine the optimal distance radius for spatial statistics analysis. This tool measures spatial autocorrelation for a series of distances and optionally creates a line graph of those   distances and their corresponding z-score. The statistically significant peak z-scores indicate distance where spatial processes promoting clustering are most pronounced. The distance corresponding to the first peak z-score was adopted for spatial statistics analysis by the default recommendation. The Getis-Ord General G statistic was calculated using the fixed Euclidean distance method with an optimal threshold distance by the 'Hot-Spot Analysis' Toolbox in ArcGIS. The Anselin local Moran's I was calculated using the inverse distance method with the 'Cluster and Outlier Analysis' Toolbox. The existence of a dispersed or a clustered distribution of riparian PSI based on a z-score (the standard deviation about the mean) where the p-value of statistical significance is typically set at 0.05 or 0.01 (a 95% or 99% confidence interval (CI)). Positive z-scores indicated clustering of large PSI values and negative z-scores indicated clustering of small PSI values (e.g., a z-score of 2.58 indicated at the 99% CI BEUs with large neighbors, and −2.58 indicated BEUs with small neighbors).

Results
Riparian PSI. The BEUs-based PSI results of the entire riparian zone were illustrated in Fig. 5. The mean and standard deviation of PSI, and number of BEUs for each measurement section were summarized in Table 5. The riparian condition had a significant difference among 13 measurement sections due to the range of PSI values. The dominated land cover types of riparian zones in section A and F is build-up category, and the riparian zones lost its ecological functions, resulted in PSI values below 45, which was lower than other measurement sections. In addition, the two sections with the smallest standard deviation (i.e., less than 2) indicated that the riparian zones were completely developed by human induced build-up structures. Section D has been less disturbed by human activities, and resulted in a better riparian condition, with the highest PSI value of over 70. However, the section also possessed the highest standard deviation up to 10.62, indicating a wide difference within the section. The PSI in section I and J was over 65 and the standard deviation was less than 3, indicating that the riparian zones have been in a homogeneous and stable condition. Section B with the 56.74 PSI value but a high standard deviation (5.20), indicating a weak stability of the riparian zone.
Accuracy validation. The comparison between PSI derived from field measurements and remote sensing observation was implemented both in measurement sites and measurement sections to validate the evaluating results (Fig. 6). The trend of the riparian PSI was consistent in both measurement sites and measurement sections, and the riparian PSI values were between 38 and 85. Specifically, there were a little differences in some measurement sites and sections. The PSI differences between field measurements and remote sensing observations for all sites were mainly distributed between 0-10, with the exception of the eleventh measurement site, which produced the highest difference (15.84). Measurement sites with difference less than 5 account for over 44% of all measurement sites. The PSI differences between field measurements and remote sensing observations were less than 8 for all measurement sections. The measurement sections with the difference less than 3 account for over 46.2%. The tenth measurement section produced the highest difference of 8.52.
The normality test indicated that the PSI scores derived from 45 measurement sites represented a normal distribution (Shapiro-Wilk statistic, p = 0.51). The field and remote sensing derived measurements were compared with Students t-test and found that there was no significant difference between the two measurements (p = 0.16) ( Table 6). Spatial statistic and analysis of riparian PSI The Anselin local Moran's I and Getis-Ord General Gi statistics were calculated to discriminate structural difference and look for evidence of clustering throughout the riparian zone (Figs 7 and 8). Anselin local Moran's I statistic found that the high-high clusters, i.e., high PSI values with high neighbors, appeared in section C, D, F, I and L. However, the PSI value in the high cluster of section F was below 50, which was significant less than other sections. Section A and B with dominated urban and farmland area had the lowest PSI values, and did not present strong high-high clusters. The low-low clusters, i.e., low PSI values with low neighbors, presented in all measurement sections. A few cases presented outliers, in the form of low-high or high-low clusters for each sections. Within the section D, the riparian PSI value changed from low-low clusters to high-high clusters.
For the Getis-Ord Gi statistic, Negative z-scores (z < −1.65, p < 0.1) represent statistically significant clusters of low-low PSI values, and become cold spots at the 90-95% confidence interval. Positive z-scores (z > 1.65, p < 0.1) indicates statistically significant clusters of high-high PSI values, and become hot spots at the 90-95% confidence interval. Clusters were non-significant (−1.65 < z < 1.65) in neighborhoods with a wide variation in PSI values, indicating those riparian conditions were randomly distributed.
Land use patterns in riparian zone. According to the evaluation results of the riparian PSI and its spatial statistics, the land use pattern was analyzed and found that the riparian landscape in the study area experienced a significant change from 1976 to 2013 (Tables 7 and 8). The areas of build-up and farmland categories accounted for the total section had increased in every measurement section from 1976 to 2013. In addition, the average annual change rate of build-up and farmland in each section had a positive growth, with the exception of section K, which the average annual change rate was −1.4‰ per year. The combination of riparian PSI values in each measurement section found that for section A and F with low PSI values, the build-up and farmland categories comprised over 80% of the total area from 1976 to 2013. However, in the section D with the highest PSI value, forest, grassland and wetland accounted for over 60% of the total area, especially over 80% from 1976 to 1986. For section H, I, J and L with the moderate PSI values, forest, grassland and wetland comprised over 45% of the total area. The average annual change rate of build-up and farmland of section C, J and M was over 14‰, more than other sections, which indicated that these sections were seriously disturbed during the 37 years. Section A and F with the lowest average annual change rate of build-up and farmland indicated that those riparian zones were changed into the urban and farmland area.
In addition, quantification of land use pattern of the BEUs in the high-high (hot spots) and low-low (cold spots) clusters found that build-up and farmland areas in sections A, B and F both accounted for over 70% of the total area in the hot spots and cold spots, while wetland, forest and grassland of the section D comprised over 70% of the total area (Table 9). In BEUs of the hot spots, wetland of section E, H and J accounted for over 47% of the total area, and forest and grassland of section C also accounted for over 47%. Build-up and farmland of section G, I and K comprised over 50% of the total area. In BEUs of the cold spots, build-up and farmland accounted for over 50% of the total area in all sections, with the exception of section D, E, H, I and J. Wetland, forest and grassland of section H, I and J accounted for over 50% of the total area.

Discussion
The multi-metric method derived a single indicator of riparian PSI from the combination of a number of sub-indicators to evaluate riparian condition. The consistency of BEUs-based evaluation results with field measurements was validated by statistical t-test and the NS model coefficient. No significant difference between two evaluation results and the NS value of 0.63 demonstrated that remote sensing derived riparian PSI and evaluation of riparian condition was efficient, reliable and comparable. The BEUs-based method with sub-indicators calculated from remote sensing data was able to ensure objectiveness and reliability of evaluation results through reducing the concerns of accessibility in field-based investigation. As anthropogenic and natural factors drive constant change, repeated measurements over time is essential for monitor and manage the riparian zone. In this study, remote-sensing based method was able to monitor and evaluate the riparian condition repeatedly over time. BEUs-based evaluation results were able to provide a guideline for field measurements to select the typical and suitable measurement sites, and riparian protection. The BEUs-based evaluation approach developed from this representative plain river of the northeast China can be extend to the other plain rivers. As for mountainous rivers, this study provided a generalized approach for combining the indicators and sub-indicators to derive the   PSI values. Modifications of scoring criteria and weights of sub-indicators and addition of sub-indicators may be needed for specific regions. The PSI value is able to indicate riparian condition. A high PSI score indicated the land-use pattern was appropriate, and riparian ecosystem was in healthy condition. But it is necessary to consider the sub-indicators when interpreting the final PSI sores. That is because although two sites or BEUs may have the same PSI score, they may have very different sub-indicator scores for indicating environment problems. For example, low PSI score for sections A and F indicated that those sections of the frail riparian ecosystem were undermined, but did not provide information regarding specific problems. After analysis of sub-indicators scores, the high score of sub-indicator of human disturbance revealed that these section of riparian zone were disturbed by the land-use pattern. In those sections, the build-up and farmland accounted for over 80% of the total area. The forest land and wetland only made up 4-7.08%. Comparatively, sections D, I and J with high PSI scores, the wetland and forest covered over 50-60% of the total riparian zone, and farmland comprised about 30-40% of the landscape. These proportions indicated that it was essential to find a balance between riparian development and fluvial ecosystem preservation.

Number of BEUs
Field measurements were designed for site-scale assessments of the current condition of a riparian zone. With the limitation of the number of measurement sites, field-based method should be more appropriate to provide section-scale evaluation results. This study provided the BEUs-based evaluation results of riparian condition. The spatial statistics revealed variation of riparian PSI. Anselin local Moran's I and Getis-Ord Gi statistics both can identify different kinds of PSI assemblages. The low-low clusters of riparian PSI discriminated the vulnerable area of the riparian zone. The information was essential for riparian protection to identify current or potential threats to an area, and adopt corresponding protective measures. Meanwhile, the high-high clusters were equally important for riparian management to identify the parts of riparian zone with a high PSI value in each section. However, the clusters were local comparisons of PSI values. It did not indicate the BEUs with the high-high clusters had a good riparian condition. For example, section A and F both had a high-high clusters of PSI values, but analysis of respective evaluation results derived from BEUs found that the PSI values were below 50 in all BEUs. The sections of riparian zone were disturbed by human activities, and in a poor condition. This demonstrated that the combination of spatial cluster analysis and BEU-based evaluation results was able to make a correct judgment for riparian condition.
Analysis of land use pattern of the BEUs in the cold and hot spots found that the sections with dominated urban and rural area have a low PSI value due to disturbance by human activities. The sections dominated wetland and forest possessed a high PSI values. The conclusions were consistent with reported studies that land-use types and patterns in a riparian zone have a great influence on riparian condition 13,[35][36][37][38][39] . Analysis of land use data from 1978 to 2013 found the riparian landscape has a great change. The build-up and farmland categories accounted for over 35% of the total area in each section, more than other land cover types indicated that build-up and farmland have been the main human disturbance to the riparian zone. In addition, the positive average annual change rate of build-up and farmland indicated that the riparian zone had been developed by the human activities from 1976 to 2013. An appropriate land use plan is needed with emphasis on protecting the riparian ecosystem and maintaining a balance between economic development and ecosystem health in the future.
Analysis of the sub-indicator scores of bank stability found that section C, D, L and M with the highest PSI values possessed the highest scores of over 60. But section A, F and G with the lower PSI values processed higher scores (about 50) of bank stability than sections B, I and J (less than 50). This indicated that the trend of the riparian PSI values was inconsistent with the scores of bank stability. This difference revealed that river connectivity and riparian vegetation had a great influence on riparian condition. River banks in the section A, F and G had realigned or modified by human-built structures or activities (e.g., bridges, weirs, dam, walking tracks, hard footpaths and recreation access), which created hard concrete banks and weakened bank erosion, resulted in high scores. Section B, I and J in the rural area could be influenced by agricultural practices or overgrazing, and caused soils degradation and accelerated bank erosion. While human-built structures precluded the lateral connectivity between the river channel and its floodplain, changed the longitudinal continuity of natural riparian vegetation, and decreased riparian vegetation coverage. In this study, section A, and F with below 10% of vegetation cover and a number of weirs resulted with the low scores of river connectivity and percent canopy cover. These demonstrated that bank stability, river connectivity and vegetation structure are closely related with riparian quality.

Conclusions
A BEUs method of remote sensing assessment was used to calculate the riparian PSI of Songhua River, Northeast China. No significant difference between evaluation results was found when compared with field measurements.   Table 7. The main land-use types comprise percent of the total area of the section in the riparian zone from 1976 to 2013. "a": build-up and farmland; "b": forest and grassland; "c": wetland.    Table 9. The main land-use types comprise percent of the total area of the BEUs in the hot spots and cold spots. "a": build-up and farmland; "b": forest and grassland; "c": wetland.

Sections
riparian zone had an important effect on riparian condition. The dominated build-up and farmland area had a poor riparian condition and the dominated wetland and forest area had a good riparian condition. The build-up and farmland had been the main human disturbance to the riparian condition. In addition, the two land cover types increased over the past 37 years. The area of build-up and farmland had increased by 787.9 km 2 from 3476.0 km 2 in 1976 to 4263.9 km 2 in 2013.