Optimal Planting Density and Nutrient Application of Soybeans: A Case Study in Northeastern China

: In the context of the Chinese government’s policy guidance, there is black soil protection and ecological environment protection. The purpose of this paper is to solve the problem that the soil ecology of the black soil in Northeast China is changing year by year, and it is necessary to explore the sowing and fertilization strategy under the new situation; most Chinese growers rely excessively on their personal experience in the process of soybean sowing and fertilization. In this study, we used “Heihe 43” soybeans and used regression experimental design methods to analyze the effects of planting density, nitrogen, phosphorus, and potassium fertilizer application on soybean yield and to determine the optimal planting density and fertilizer ratios. The study reveals that the optimal soybean planting density in Northeast China is 45.37 × 10 4 plants/ha, with nitrogen at 98.4 kg/ha, phosphorus at 218.96 kg/ha, and potash at 47.62 kg/ha. Under these conditions, soybean yields can reach 3816.67 kg/ha. This study can provide a theoretical method for decision-making to obtain the optimal planting density and fertilizer ratio for different regions of the farming system.


Introduction
China's soybean self-sufficiency is insufficient, and it is overly dependent on imports.China's dominance in the global soybean trade, accounting for approximately 66% of it, has placed substantial pressure on global food supply chains [1].In recent years, China has significantly bolstered its agricultural support, accompanied by a continuous evolution of relevant policies.China's Central Document No. 1 of 2022 proposes to vigorously implement the soybean and oilseed production capacity enhancement project.The agricultural sector has increasingly prioritized collaboration with agricultural institutions for experiments related to crop cultivation and fertilization [2,3], with a particular focus on major grain crops like soybeans [4], maize [5], rice [6,7], and wheat [8], especially within the context of planting and fertilization experiments.According to the 2020 statistics from the Food and Agriculture Organization (FAO) of the United Nations, global agricultural fields received a total application of 85 million metric tons of nitrogen (N), 7 million metric tons of phosphorus (P), and 12 million metric tons of potassium (K).Nitrogen utilization efficiency ranged from 50% to 62% of the applied nitrogen.However, China significantly surpasses other agricultural powerhouses in fertilizer use.Excessive fertilizer application in China has escalated production costs, posing challenges to carbon peaking and carbon neutrality goals.Crop fertilizers, especially nitrogen (N), phosphorus (P), and potassium (K) fertilizers, exhibit low utilization rates.In China, nitrogen utilization efficiency (NUE) for crop production is notably lower, approximately 30%, compared to the global average [9].China faces the highest nitrogen surplus worldwide, with over 50% of nutrient substances Agronomy 2023, 13, 2902 2 of 21 applied to land as fertilizers going to waste.This wastage leads to soil compaction [10] and disrupts the native soil ecological equilibrium [11].The black soil of Northeast China is one of only four black soils in the world and one of the major grain-producing areas in China, but it is facing degradation problems caused by wind and water erosion, organic matter decline, and soil sloughing.Optimizing fertilizer application strategies can reduce the transfer of chemical fertilizers into the soil, reduce environmental pollution, and protect the soil ecology of the Northeast black soil.Implementing appropriate planting densities can enhance resource utilization [12] (e.g., light, moisture, and nutrients) by crops, thereby increasing fertilizer efficiency and reducing waste.This, in turn, can significantly enhance soybean yields [13].
Chinese agriculture grapples with challenges due to its widely dispersed land, making it difficult to establish specific fertilization standards as guidelines.Furthermore, a lack of technical support personnel leaves farmers relying solely on their individual experiences for planting and fertilizing.As a result, a contradiction exists between theoretical concepts of planting and fertilization in China and the practical realities of local agricultural practices, as well as the level of scientific and technological advancement in the region.Innovation in planting and fertilization theories is imperative to promote the healthy and sustainable development of eco-friendly agriculture.This innovation must be grounded in the local context, considering factors such as soil conditions, crop varieties, and production costs [7].
In 2008, Oz, M. conducted a field experiment to investigate the impact of soybean planting density and nitrogen application rates on soybean yield.The experiment featured the soybean variety A-3935 and was conducted at four different planting densities and four distinct nitrogen fertilizer application levels.The study assessed parameters such as plant height, minimum pod height, per-plant branching count, per-plant pod count, per-plant yield, harvest index, hundred-seed weight, and soybean production [14].In 2009, Goncalves, R. J. D. and Abreu, A. D. B. conducted field experiments to evaluate the performance of various common soybean cultivars under varying growth habits, different fertilizer application rates, and various planting densities.Subsequently, they developed corresponding fertilization strategies [15].In 2016, Ferreira, A. S. et al. conducted two field experiments using a randomized complete block design to assess the influence of nitrogen on soybean grain yield, yield components, as well as oil and protein concentrations across various planting densities and two nitrogen fertilizer levels [16].In 2018, a field experiment was conducted to evaluate the interactive effects of seeding rate and phosphorus (P) and potassium (K) fertilization on growth, grain yield, and protein and oil content in soybeans [17].In 2021, Xu, C. L. et al. conducted a two-year field experiment to investigate the impact of planting density on soybean plant growth and yield components.They analyzed the relationship between photosynthetic rates, dry matter accumulation, and yield under different planting densities and plant distributions [18].In 2021, Salvagiotti F and Magnano L curated an extensive dataset spanning from 2009 to 2018, encompassing critical agricultural parameters, including soybean seed yield, total biomass at physiological maturity, and data on the absorption of essential nutrients such as nitrogen (N), phosphorus (P), potassium (K), and sulfur (S) [19].
Some scholars have primarily focused their research on the impact of soybean yield from various perspectives, including cultivar improvement [20,21], crop rotation patterns, and tillage methods [22][23][24][25], fertilization and nutrient management [26,27], pest and disease control, as well as temperature and precipitation [28,29].In recent years, China's black soil ecosystems have deteriorated due to soil erosion and land degradation, particularly soil environmental pollution resulting from the use of chemical substances such as fertilizers and pesticides in agricultural production [30,31].To investigate optimal planting density and fertilizer application rates under diverse soil conditions, this study addresses the prevalent practice among the majority of Chinese growers, who traditionally rely on empirical methods for both fertilization and seeding.Often, these practices exceed optimal thresholds for fertilizer application, resulting in seed and fertilizer waste, soil compaction [10], environmental pollution, increased production costs, and adverse effects on both the quality and yield of harvests.This research aims to tackle the issues impacting the progress of sustainable agricultural development in China.In the experimental zone of Shanhe Farm in Nenjiang County, Heihe City, Heilongjiang Province, located in the northeastern region of China, foundational data were collected using experimental methods.Employing nonlinear regression analysis, we constructed a seeding and fertilization model with soybean yield as the optimization objective and three types of fertilizer application rates and planting density as experimental factors.Subsequently, this model underwent optimization to determine the optimal planting density, nitrogen, phosphorus, and potassium fertilizer application rates that maximize soybean yield.This research endeavor aims to provide decision-making support for seeding and fertilization practices, guiding judicious fertilizer application at the local level and furnishing a framework for achieving yield maximization.Rational seeding and fertilization not only unlock the full potential of crops and reduce costs but also contribute to the equilibrium of the soil ecosystem [31].

Experimental Design
In this study, the experimental site was carefully chosen within the Shanhe Farming Experimental Zone, located in Nenjiang County, Heihe City, Heilongjiang Province, China.The coordinates of this site are approximately 49 • 15 54 N latitude and 125 • 44 53 E longitude.This site boasts annual average precipitation ranging from 480 to 490 millimeters and enjoys a solar radiation duration of 2500 to 2800 h per year.Moreover, the annual average accumulated temperature above 0 • C, known as effective growing degree days, hovers around 1900 • C, with a frost-free period extending approximately 110 to 120 days.Notably, the Shanhe Farming Experimental Zone is characterized by flat terrain and primarily comprises meadow black calcareous soil.The soil organic matter content within this zone falls within the range of 4% to 6%.On 25 October 2019, the experimental field was plowed and harrowed.On 20 April 2020, plots were divided in the experimental field.The dimensions of the experimental plot measure 139 m in length and 33 m in width, with specific parameter settings outlined in Table 1.Throughout this research, "Heihe 43" soybean seeds were employed-a cultivar consistently cultivated within the Shanhe Farming area.The germination period of soybean "Heihe 43" lasts 5-7 days, the growing period is 80-120 days, the flowering period is 20-30 days, and the fruiting period is 30-40 days.In this study, all plots in the same block group were sown on one day, and the field was sown on 28 April 2020.Nitrogen fertilizer was administered in the form of urea (46% N), phosphorus fertilizer as ammonium phosphate (18% N and 46% P), and potassium fertilizer as potassium sulfate (60% K).The germination rate of the seeds reached 93%.Cultivation management in this field trial started on 4 May 2020 and ended on 28 September 2020, which was about 5 months.The harvesting date of this trial was from 28 September to 1 October 2020, and the soybeans were harvested by hand.This experiment adhered to international standards and best practices, with the selected site being ideally suited to fulfill the study's objectives.The experimental area was partitioned into 26 primary zones aligned with the plowing direction.Starting from the outermost perimeter, every three ridges were designated Agronomy 2023, 13, 2902 4 of 21 as subzones, each surrounded by protective rows.As a result, each primary zone was further subdivided into five subzones, resulting in a total of 130 subzones.Two subzones were intentionally left vacant at one end of the experimental site, while the remaining 128 subzones were meticulously utilized for implementing the 128 experimental treatments.This allocation also determined the placement of the aisles.The precise layout of the experimental field is depicted in Figure 1.The experimental area was partitioned into 26 primary zones aligned with the plowing direction.Starting from the outermost perimeter, every three ridges were designated as subzones, each surrounded by protective rows.As a result, each primary zone was further subdivided into five subzones, resulting in a total of 130 subzones.Two subzones were intentionally left vacant at one end of the experimental site, while the remaining 128 subzones were meticulously utilized for implementing the 128 experimental treatments.This allocation also determined the placement of the aisles.The precise layout of the experimental field is depicted in Figure 1.

Field Experimental Design
When harvesting a specific crop, it holds the potential to impact or even disrupt neighboring crops, thereby exerting an influence on the final yield [32].This research employed a block sampling approach to minimize field disturbance during crop field trials.Specifically, contiguous blocks within the field were selected for indoor seed examinations and other assessments to minimize overall crop disturbance.A comparison with nonblock sampling methods, exemplified in Figure 2a, clearly highlights that the depicted technique resulted in the highest degree of plant disruption, yielding a maximum disruption ratio of 9 between the disturbed plant area and the sample size.Conversely, utilization of the block sampling method with a sample shape approximating a square led to the least overall plant disturbance.For instance, in the case of grid-based systematic sampling, as illustrated in Figure 2b, the disruption ratio was reduced to 2.78 [33].Therefore, this study used the sampling method in Figure 2b.

Field Experimental Design
When harvesting a specific crop, it holds the potential to impact or even disrupt neighboring crops, thereby exerting an influence on the final yield [32].This research employed a block sampling approach to minimize field disturbance during crop field trials.Specifically, contiguous blocks within the field were selected for indoor seed examinations and other assessments to minimize overall crop disturbance.A comparison with non-block sampling methods, exemplified in Figure 2a, clearly highlights that the depicted technique resulted in the highest degree of plant disruption, yielding a maximum disruption ratio of 9 between the disturbed plant area and the sample size.Conversely, utilization of the block sampling method with a sample shape approximating a square led to the least overall plant disturbance.For instance, in the case of grid-based systematic sampling, as illustrated in Figure 2b, the disruption ratio was reduced to 2.78 [33].Therefore, this study used the sampling method in Figure 2b.
Diverse factors, encompassing climate, crop types, and geographical and historical environments, contribute to variations in soil nutrient content and ecosystem characteristics [34,35].To minimize experimental error, the entire test field in this study was subdivided into multiple large plot groups.In the preceding year, maize was cultivated in the experimental field.The experimental design and implementation plan for the selected plots can be found in Table 2. Climate variability exerts a profound influence on both the magnitude and stability of crop yields [36].Factors such as sunlight exposure [37-40], temperature [28], and precipitation [29,41] significantly impact soybean production.In September 2019 (one year before the official trial), this study pre-treated the trial plots using farm machinery and equipment, including deep plowing, fine harrowing, and leveling, to make the condition of the soil in each plot at the time of the experiment essentially the same.As a result of these treatments, it was assumed in this study that the organic matter content of the soil in each plot was uniform.The organic matter content of the soil (including sunshine, temperature, precipitation, etc., in each plot in the experimental field) was used as the same experimental conditions in this study, in this context, to study soybean yield about the number of inorganic fertilizers (urea, diamine, potassium sulfate) applied and the planting density.
The experiment took place in 2020 at the Shanhe Agricultural Experimental Zone of Nenjiang County, Heihe City, Heilongjiang Province, China.Precipitation in Nenjiang County is predominantly concentrated from June to September, accounting for approximately 80% of the annual precipitation.Notably, the months of July and August experience the highest levels of precipitation.Conversely, the period from November to February exhibits cold climatic conditions, rendering it unsuitable for agricultural production activities.Temperature patterns, including average high and low temperatures, as well as total precipitation for Nenjiang County from March to October, are illustrated in Figure 3.
tics [34,35].To minimize experimental error, the entire test field in this study was subdivided into multiple large plot groups.In the preceding year, maize was cultivated in the experimental field.The experimental design and implementation plan for the selected plots can be found in Table 2. Climate variability exerts a profound influence on both the magnitude and stability of crop yields [36].Factors such as sunlight exposure [37][38][39][40], temperature [28], and precipitation [29,41] significantly impact soybean production.
In September 2019 (one year before the official trial), this study pre-treated the trial plots using farm machinery and equipment, including deep plowing, fine harrowing, and leveling, to make the condition of the soil in each plot at the time of the experiment essentially the same.As a result of these treatments, it was assumed in this study that the organic matter content of the soil in each plot was uniform.The organic matter content of the soil (including sunshine, temperature, precipitation, etc., in each plot in the experimental field) was used as the same experimental conditions in this study, in this context, to study soybean yield about the number of inorganic fertilizers (urea, diamine, potassium sulfate) applied and the planting density.
The experiment took place in 2020 at the Shanhe Agricultural Experimental Zone of Nenjiang County, Heihe City, Heilongjiang Province, China.Precipitation in Nenjiang County is predominantly concentrated from June to September, accounting for approximately 80% of the annual precipitation.Notably, the months of July and August experience the highest levels of precipitation.Conversely, the period from November to February exhibits cold climatic conditions, rendering it unsuitable for agricultural production activities.Temperature patterns, including average high and low temperatures, as well as total precipitation for Nenjiang County from March to October, are illustrated in Figure 3.This study was divided into a total of four types of regression experimental design: single-factor experimental design, four groups; quadratic orthogonal rotating combination design, two groups of replicates; D-optimal regression design, one group; and random regression, one group.The single-factor regression design experiments consisted of four groups as follows: nitrogen, phosphorus, and potassium three-factor regression design experiments numbered 98~107, 108~117, and 118~127 in that order, and sowing density experiments numbered 88~97; two groups of replicated experimental design numbered 16~51 and 52~87; D-optimal experiments numbered 1~15; and a blank control experiment numbered 128.

Collection of Basic Data
This experiment was conducted in 2018 and 2019 as a two-year pre-experiment aimed at obtaining the optimal range of values for the single-factor test and laying the foundation for determining the 0-level values of nitrogen, phosphorus, and potassium fertilizer ap-plication and planting density in the second-order orthogonal rotating composite design.In 2020, the experimental steps of harvesting, threshing, and indoor drying and weighing were used to obtain the basic data of the field trial, construct the regression equation model, carry out the quantitative analysis, and optimize the solution [42].In this study, models were tested, constructed, and optimized using Design Expert 13 software for two sets of replicated experimental data from quadratic orthogonal rotated combinatorial design and one set of experimental data from D-optimal regression.They fitted to one-factor four sets of data using Origin 2022 software.

Single-Factor Regression Design Implementation Data
Fundamental research data concerning soybean yield were obtained by conducting measurements within specific land parcels classified by unique numerical identifiers.Figure 4 depicts the correlation between soybean yield, planting density, and fertilization levels.By examining the relationship diagram presented in Figure 4, which illustrates the correlation between planting density, fertilizer application rates, and soybean yield, it becomes evident that an incremental increase in planting density and fertilizer application rates leads to a gradual enhancement in soybean yield, ultimately reaching the maximum potential for soybean growth.Simultaneously, when planting density surpasses 46.35 × 10 4 plants/ha, the application rates of nitrogen, phosphorus, and potassium exceed 69.81 kg/ha, 169 kg/ha, and 25.29 kg/ha, respectively.Beyond these thresholds, the yield gradually diminishes, becoming a limiting factor.The experimental design employed to assess the unifactorial impact of planting den- By examining the relationship diagram presented in Figure 4, which illustrates the correlation between planting density, fertilizer application rates, and soybean yield, it becomes evident that an incremental increase in planting density and fertilizer application rates leads to a gradual enhancement in soybean yield, ultimately reaching the maximum potential for soybean growth.Simultaneously, when planting density surpasses 46.35 × 10 4 plants/ha, the application rates of nitrogen, phosphorus, and potassium exceed 69.81 kg/ha, 169 kg/ha, and 25.29 kg/ha, respectively.Beyond these thresholds, the yield gradually diminishes, becoming a limiting factor.

Data Processing and Model Optimization in Regression Experiments on Planting Density Effects
The experimental design employed to assess the unifactorial impact of planting density on crop yield underwent rigorous statistical analysis.The results of the parameter tests about the influence of planting density on crop yield are detailed in Table 3.Based on the data presented in Table 3, the column depicting variation probabilities, it is evident that the overall variation probability of the model is less than 0.0001.This finding indicates the significance of the relationship between soybean yield and planting density factors.Conducting examinations on other model parameters allows us to derive conclusions with a high degree of concordance between the regression equation and actual observations.A summary of the obtained model can be found in Table 4, while the model coefficients and corresponding t-test results are provided in Table 5.Whereas x 1 represents planting density (10 4 plants/ha), x 2 signifies nitrogen fertilizer application rate (kg/ha), x 3 denotes phosphorus fertilizer application rate (kg/ha), and x 4 signifies potassium fertilizer application rate (kg/ha).
Obtained regression equation: Based on the actual conditions at Shanhe Farm, the range of values for planting density was determined.A model was constructed using a regression equation [43] that captures the relationship between planting density and yield.Subsequently, the model was subjected to mathematical analysis and solution.y = −0.40x6 presents the parameter testing results for assessing the influence of nitrogen (N) fertilizer application on crop yield.These results were determined through a nitrogen efficiency regression trial design and subsequent calculations.The data presented in Table 6 demonstrates that the variation probability, as indicated in the corresponding column, is less than 0.0001.This finding suggests a substantial correlation between soybean yield and nitrogen fertilizer factors.Upon closer examination of the model, it becomes evident that the regression equation aligns favorably with realworld observations.
The obtained model summary is presented in Table 7, while the model coefficients and t-test results are available in Table 8.The regression equation for determining the relationship between nitrogen fertilizer application rates and crop yields: Considering the particular conditions at the Shanhe Farming Experimental Zone, we determined a range of nitrogen fertilizer application values.Afterward, we developed a regression equation using the collected yield data, which was then followed by computational modeling and solution.ŷ = −4.47x When all other influencing factors are maintained at zero levels, the optimal nitrogen fertilizer application rate is determined to be 69.81kg/ha, resulting in a maximum yield of 3157.39 kg/ha.

Data Processing and Model Optimization in Regression Experiments on Phosphorus (P) Efficiency
After performing calculations about the experimental design for phosphate (P) efficiency regression, we present the model parameters employed to examine the influence of phosphate fertilizer application rates on crop yield in Table 9.According to the data presented in Table 9, the column depicting variation probabilities reveals an overall variation probability of approximately 0.0002 for the model.This suggests a significant relationship between soybean yield and the phosphorus fertilizer factor.Further examination of other model parameters yields a favorable conclusion regarding the model's regression equation and its alignment with real-world conditions.
The model summary can be found in Table 10, while the model coefficients and t-test results are presented in Table 11.In this academic paper, we present the regression equation that illustrates the correlation between the application rates of phosphorus fertilizer and crop yields.ŷ = −0.02x 3 2 + 5.34x 3 + 2520.03(7) Based on the actual conditions of the Shanhe Farming Experimental Zone, the range of phosphorus fertilizer application rates was determined.A regression equation was constructed based on the obtained yield data, and subsequent computational analysis was performed to derive the model.ŷ = −0.02x 3 2 + 5.34x 3 + 2520.030 < x 3 < 262.50 (8) Solution: Yield ŷ = 2973.84Phosphate f ertiliser x 3 = 169 (9) When maintaining all other variables at their baseline levels, the optimal rate of phosphorus fertilizer application is established at 169 kg/ha.This point corresponds to the maximum observed yield, reaching a value of 2973.84 kg/ha.

Data Processing and Model Optimization in Regression Experiments on Potassium (K) Efficiency
Table 12 presents the results of the analysis of model parameters, which assess the influence of potassium (K) fertilizer application rates on crop yield.These assessments were conducted using a potassium efficiency response regression experimental design.As indicated by the data in Table 12, the variation probability column, the overall variation probability of the model is approximately 0.0001.This result suggests that there is a significant relationship between soybean yield and potassium fertilizer factors.Further examination of other aspects of the model yields the conclusion that the regression equation exhibits a good fit with real-world conditions.
The model-derived summaries can be found in Table 13, while the model coefficients and the results of t-tests are presented in Table 14.Based on the actual conditions at Shanhe Farming Experimental Zone, the range of potassium fertilizer application rates was determined.A regression equation for yield was established using the obtained data and subsequently solved through mathematical modeling.ŷ = −1.34x With all other factors at their baseline values, when the potassium fertilizer application rate is 25.29 kg/ha, at which point the maximum yield reaches 3125.48 kg/ha.The single-factor regression design only provides insights into the individual impacts of planting density, nitrogen fertilizer, phosphorus fertilizer, and potassium fertilizer on soybean yield.However, during the soybean cultivation process, these four factors exhibit interactions with each other.To investigate the combined effects of these interacting factors on soybean yield, we employed a quadratic orthogonal rotation regression design to compute the coded levels and spacing values for each factor.In contrast to single-factor regression experimental designs, the quadratic orthogonal rotation regression design introduces interactions between pairwise factors, with a focus on representing a comprehensive surface model within the coding range.

Second-Order Orthogonal
Under the coding scheme, where 'A' represents planting density, 'B' represents nitrogen fertilizer, 'C' represents phosphorus fertilizer, and 'D' represents potassium fertilizer; under the actual value scheme, where 'x 1 ' represents planting density, 'x 2 ' represents nitrogen fertilizer, 'x 3 ' represents phosphorus fertilizer, and 'x 4 ' represents potassium fertilizer, quadratic square terms were centered as follows: In this study, "n" represents the number of experiments, with "n" equal to 36.Specifically, 12 repetitions of experiments were conducted at the center point of the 0-level, resulting in the establishment of the experimental design matrix.Subsequently, based on the quadratic orthogonal rotation regression design matrix, an experimental plan implementation schedule was developed, and a partial sample table is presented in Table 15.The data from the quadratic orthogonal rotational regression tests were processed, and the model parameters and parameter test tables obtained are shown in Table 16.As shown in the F-value column of Table 16, the influence factors of the four factors on yield are as follows: planting density is 4.73, nitrogen fertilizer is 48.85, phosphorus fertilizer is 15.18, and potash fertilizer is 3.08.It can be seen that the magnitude of the factors contributing to yield is nitrogen fertilizer (N) > phosphorus fertilizer (P) > planting density > potassium fertilizer (K).
As indicated in Table 16, the variance probabilities of the model as a whole are less than 0.0001, signifying the significant regression relationships between the yield and the variables x 1 , x 2 , x 3 , and x 4 .Conducting an overall model, the F-test reveals that the Lack of Fit of model inadequacy is not significant.This suggests that the derived equation aligns well with the actual scenario and is deemed valid.
Under the coding scheme, where A represents planting density, B signifies nitrogen fertilizer, C represents phosphorus fertilizer, and D stands for potassium fertilizer, we derive the yield regression equation within the framework of the coding system: We established the model by removing non-significant terms (p > 0.05) from the yield regression equation.These excluded terms encompassed the first-order effect of potassium, as well as the interaction terms between planting density and nitrogen, planting density and phosphorus, and phosphorus and potassium, along with the second-order effect of potassium.Due to their orthogonality, it is evident that these exclusions do not influence the remaining coefficients.Subsequently, we conducted optimization calculations.The results of the variance analysis and F-test are presented in Table 17.Considering the presence of nitrogen elements in the diamine-based fertilizers employed, due consideration was given to this factor during the execution of the experiments.An inconsistency was observed between the numerical values of fertilizer application quantities in the experimental data table and the levels of the factors.As a result, the model was constructed using factor levels as independent variables to ensure alignment with the experimental design scheme.The ranges for the application rates of nitrogen, phosphorus, and potassium fertilizers, as well as planting density, were determined based on the actual conditions at Shanhe Farm.A coding scheme was established using the obtained yield regression equation and subsequently solved.Finally, actual values for fertilizer application rates and planting density were obtained through numerical conversion.
In the context of agricultural research, it has been observed that the nonlinear fitting model yields superior results.Specifically, when considering planting density, nitrogen fertilizer application, phosphorus fertilizer application, and potassium fertilizer application at values of 45.61 × 10 4 plants/ha, 104.42 kg/ha, 202.31 kg/ha, and 43.35 kg/ha, respectively, the maximum crop yield is achieved at 3703.76 kg/ha.This underscores the efficacy of nonlinear modeling techniques in optimizing crop yield under varying agricultural input conditions.

Data from the Second Set of Experimental Implementations in a Second-Order Orthogonal Rotational Composite Design
To mitigate the inadvertent errors introduced by experimental randomness, we implemented a strategy involving the conduction of duplicate trials across two distinct experimental groups.Table 18 presented herein comprises a subset of the data derived from the execution of the second set of trials.Upon processing the data from the secondary orthogonal rotational regression experiment, the model parameters and parameter validation table can be found in Table 19.As depicted in Table 19, the impact factors of the four factors on yield, denoted as F-values, are as follows: planting density at 9.45, nitrogen fertilizer at 83.57, phosphorus fertilizer at 45.73, and potassium fertilizer at 0.0111.Consequently, it can be inferred that the contribution order of these factors to yield is as follows: nitrogen fertilizer (N) > phosphorus fertilizer (P) > planting density > potassium fertilizer (K).
As shown in Table 19, the variation probability in the column for Model 19 is less than 0.0001, indicating an exceedingly significant regression relationship between yield and factors A, B, C, and D. Conducting an overall model F-test reveals that the variability probability of the model's lack of fit is not significant.This suggests that the obtained equation aligns well with the actual circumstances and is indeed applicable. The By the yield regression equation, a model was constructed by eliminating non-significant terms.Specifically, the interactive terms between planting density and nitrogen, planting density and phosphorus, and phosphorus and potassium were removed, along with the quadratic term for potassium.Given their orthogonality, it is evident that the removal of these terms does not impact the other coefficients.The variance analysis results and F-test for the second set of experiments are presented in Table 20.By the actual conditions observed at the Mountain River Farm, we determined the ranges of nitrogen (N), phosphorus (P), and potassium (K) fertilizer application rates, as well as planting densities.Subsequently, we established an encoding-based modeling framework utilizing regression equations derived from the obtained crop yields.These equations were subsequently solved to derive practical values for both fertilizer application rates and planting densities through numerical conversion.
In the realm of agricultural research, it has been observed that nonlinear regression models exhibit superior performance in optimizing crop yield.Notably, when the planting density, nitrogen fertilizer application rate, phosphorus fertilizer application rate, and potassium fertilizer application rate are set at 45.37 × 10 4 plants/ha, 98.4 kg/ha, 218.96 kg/ha, and 47.62 kg/ha, respectively, the maximum achievable yield stands at 3816.67 kg/ha.This outcome underscores the significance of employing nonlinear fitting techniques in agricultural yield optimization studies.

Implementation Data for D-Optimal Regression Experimental Design
In this study, we conducted an experimental design to optimize soybean yield.Given the high requirement for the goodness of fit of the regression equation and the challenge of dealing with a large number of experimental treatments, we employed the D-optimization regression approach as delineated in prior works [44][45][46].This method not only ensures precise parameter estimation but also reduces the requisite number of experimental data points, thus resulting in a cost-effective approach from an experimental perspective.A subset of the experimental designs and corresponding data samples are presented in Table 21 for reference.

Data Processing and Model Optimization in D-Optimal Regression Experiments
Processing of D-optimal regression experimental data yields model parameters and parameter test results, as shown in Table 22.In light of the actual conditions observed at the Shanhe Farm, we derived the ranges for nitrogen (N), phosphorus (P), and potassium (K) fertilizer application rates as well as planting densities.Subsequently, a regression equation for yield was established, encoded as Model (20), based on the data presented in Table 22.This model was then subjected to computational analysis, culminating in the attainment of practical values for fertilizer application rates and planting densities through numerical conversion.
Because the design scheme falls under a saturated design and no replicates have been introduced, the remaining degrees of freedom are zero.Consequently, conventional methods for conducting significance tests on the equation are not applicable.Instead, the goodness-of-fit Chi-squared test method, comparing observed values to predicted values, was employed for verification.
Following the Chi-squared test between the predicted values and the observed values using Design Expert 13, it was observed that the predicted values align well with the observed values.This observation suggests that the regression equation is by real-world conditions.Thus, we obtain the following result: In the realm of agricultural research, it has been empirically established that nonlinear regression models yield superior fitting outcomes.Specifically, when considering the variables of planting density, nitrogen application, phosphorus dosage, and potassium supplementation at rates of 35 × 10 4 plants/ha, 118.38 kg/ha, 201.09 kg/ha, and 47.15 kg/ha, respectively, the resultant maximum yield attains 3482.76 kg/ha.This finding underscores the efficacy of nonlinear regression as an optimal approach for modeling and predicting crop yield responses to various agronomic inputs.

Discussion
In this study, a comprehensive analysis was conducted to assess the impact of four factors: nitrogen (N), phosphorus (P), potassium (K) fertilizer application rates, and planting density on soybean yield.Additionally, the study examined the effects of individual factors, pairwise combinations of factors, and multiple factors and their interactions on soybean yield.The order of influence of these four factors on soybean yield was as follows: nitrogen fertilizer (N) > phosphorus fertilizer (P) > planting density > potassium fertilizer (K).It was also observed that there was a significant interaction between fertilizer application rates and planting density.Single-factor experiments in this study indicated that soybean yield decreased with increasing planting density and nitrogen application.
When considering the interactions between factors, orthogonal rotation combination experiments were performed at a significance level of α = 0.05.In the first group of experiments, it was found that the factor contribution rate of the interaction between N and P was 37.55, which was greater than the factor contribution rate of the interaction between planting density and K (22.33) and the interaction between N and K (20.72).In the second group of experiments, the factor contribution rate of the interaction between planting density and K was 45.16, higher than the factor contribution rate of the interaction between N and P (27.15) and the interaction between N and K (21.92).Both sets of experiments indicated that the interaction between planting density and K, as well as the interaction between N and P, had a greater impact on soybean yield compared to the interaction between N and K. Furthermore, interactions involving planting density with N, planting density with P, and P with K did not have a significant impact on soybean yield at a significance level of α = 0.05.
In addition, this study concluded that the influence of nitrogen on soybean yield was greater than that of phosphorus.Salvagiotti F's experiments highlighted the critical importance of crop yield and nutrient uptake rates for achieving a more sustainable farming system [19].Zhou Jing et al. conducted field experiments on the growth of soybeans and wheat in the northeastern region of China over two consecutive crop seasons [47].The findings of these aforementioned studies align consistently with the outcomes of our experiment.
The limitations of conducting seeding and fertilization experiments only once or twice are apparent, particularly in regions like Heilongjiang Province, China, which experiences a single annual cropping season.It takes several years to refine a model for a specific crop, a characteristic inherent to field experiments.In further collaboration with the Shanhe Farm, we plan to undertake multi-year soybean seeding and fertilization experiments at Shanhe Farm's test fields.The accumulated data will be instrumental in enhancing the model derived from this research.As the data gradually improves, optimization objectives such as 'maximizing economic efficiency', 'facilitating management', or 'enhancing ecological conditions' can be incorporated into the model.Additionally, future research will focus on fertilization ratios under different regional and soil conditions.Furthermore, numerous factors can influence final crop yields, including uncontrollable elements such as light intensity, precipitation, and pest species.In-depth investigations in these areas represent promising directions for future research.Additionally, this study employed quadratic modeling for the data.When experimental conditions are constrained, it may be worthwhile to explore alternative modeling approaches in future investigations.

Conclusions
In the context of advocating for sustainable agricultural development in China, emphasizing low carbon emissions and the "Three Reductions" initiative, this study focuses on soybeans as its subject of investigation.We have designed an experimental scheme to examine the relationship between soybean nutrients, planting density, and yield.The experimental trials were conducted at the Shanhe Farm in Nenjiang County, Heihe City, Heilongjiang Province, China, involving seeding and fertilization procedures.Based on the results obtained from these seeding and fertilization experiments, we have constructed a soybean planting model, which defines the relationship between yield, planting density, and the application of three key fertilizers.Subsequently, we optimized this soybean planting model.By aiming to maximize yield, we have identified a cultivation strategy that offers relative superiority.This optimal combination entails a planting density of 45.37 × 10 4 plants/ha, nitrogen application at 98.4 kg/ha, phosphorus application at 218.96 kg/ha, and potassium application at 47.62 kg/ha.The significance of the four experimental factors on soybean yield is ranked as follows: nitrogen fertilizer > phosphorus fertilizer > planting density > potassium fertilizer.This research not only establishes a theoretical foundation for practical seeding and fertilization practices but also furnishes a scientific basis for determining the optimal planting density and fertilizer ratios in different geographical regions' agricultural cultivation systems.Furthermore, it serves as a crucial decision-making tool to guide soil testing and tailored fertilization strategies, ultimately leading to the maximization of crop yields.Rational seeding and fertilization practices, as outlined in this study, not only harness the full potential of crops but also curtail production costs while preserving the ecological balance of soil ecosystems.This approach aligns seamlessly with China's sustainable development strategy.

Figure 2 .
Figure 2. Comparison of two sampling methods.Figure 2. Comparison of two sampling methods.(a) non-block sampling, (b) grid-based systematic sampling.

Figure 2 .
Figure 2. Comparison of two sampling methods.Figure 2. Comparison of two sampling methods.(a) non-block sampling, (b) grid-based systematic sampling.

Figure 3 .Figure 3 .
Figure 3. Annual temperature and total precipitation data for Nenjiang County in 2020.

Figure 4 .
Figure 4. Illustrative graph depicting the relationship between planting density, fertilizer application, and soybean yield.

Figure 4 .
Figure 4. Illustrative graph depicting the relationship between planting density, fertilizer application, and soybean yield.
Rotational Composite Design 3.3.1.Data from the First Set of Experimental Implementations in a Second-Order Orthogonal Rotational Composite Design

Table 1 .
Experimental site parameter configuration table.

Table 2 .
Table of partial plot experimental design and implementation plan.

Table 3 .
Model parameter examination table for the relationship between planting density and yield.

Table 4 .
Model summary of the relationship between planting density and yield.

Table 5 .
Table of model coefficients and parameter tests.
Within the experimental parameters of this study, the optimal planting density for achieving the highest yield was determined to be 46.35× 10 4 plants/ha, resulting in a peak yield of 2902.85 kg/ha.

Table 6 .
Table of the model parameters of the relationship between nitrogen (N) factor and yield.

Table 7 .
Model summary of the relationship between nitrogen (N) factor and yield.

Table 8 .
Model coefficient and parameter test table of the relationship between nitrogen (N) factor and yield.

Table 9 .
Table of the model parameters of the relationship between phosphorus (P) factor and yield.

Table 10 .
Model summary of the relationship between phosphorus (P) factor and yield.

Table 11 .
Model coefficient and parameter test table of the relationship between phosphorus (P) factor and yield.

Table 12 .
Table of the model parameters of the relationship between potassium (K) factor and yield.

Table 13 .
Model summary of the relationship between potassium (K) factor and yield.

Table 14 .
Model coefficient and parameter test table of the relationship between potassium (K) factor and yield.

Table 15 .
Partial actualizing data of design test plan.

Table 16 .
Quadratic orthogonal rotation regression design of the first set of the experimental model parameter test table.

Table 17 .
Table of quadratic orthogonal rotation design of the first set of test model parameter tests.

Table 18 .
Partial actualizing data of design test plan.

Table 19 .
Quadratic orthogonal rotation regression design of the second set of experimental model parameter test tables.

Table 20 .
Table of quadratic orthogonal rotation regressive design replication to yield model parameter test.

Table 21 .
Actualizing data of design test plan.

Table 22 .
D-optimal regression design model parameter test table.