【农气项目】基于适宜度的产量预报

直接上干货（复制到开发工具即可运行的代码）

1. 适宜度模型及作物適宜度计算方法

2. 产量分离

3. 基于适宜度计算产量预报

1. 适宜度模型及作物適宜度计算方法

    // 三基点温度配置private final double tempMin;private final double tempOpt;private final double tempMax;// 降水/日照理想值private final double idealPrecipitation;private final double idealSunshineRatio;// 权重配置private final double tempWeight;private final double precWeight;private final double sunWeight;public SimpleCropSuitability(double tempMin, double tempOpt, double tempMax,double idealPrecipitation, double idealSunshineRatio,double tempWeight, double precWeight, double sunWeight) {// 参数校验if (tempMin >= tempOpt || tempOpt >= tempMax) {throw new IllegalArgumentException("温度参数必须满足: tMin < tOpt < tMax");}if (idealPrecipitation <= 0 || idealSunshineRatio <= 0) {throw new IllegalArgumentException("降水和日照理想值必须大于0");}this.tempMin = tempMin;this.tempOpt = tempOpt;this.tempMax = tempMax;this.idealPrecipitation = idealPrecipitation;this.idealSunshineRatio = idealSunshineRatio;this.tempWeight = tempWeight;this.precWeight = precWeight;this.sunWeight = sunWeight;}// 温度适宜度（三基点模型）private double calculateTempScore(double meanTemp) {if (meanTemp <= tempMin || meanTemp >= tempMax) {return 0;}return (meanTemp < tempOpt)? (meanTemp - tempMin) / (tempOpt - tempMin)  // 低温区间线性增长: 1 - (meanTemp - tempOpt) / (tempMax - tempOpt); // 高温区间线性衰减}// 降水适宜度（新线性模型）private double calculatePrecipitationScore(double precipitation) {return Math.min(precipitation / idealPrecipitation, 1.0);}// 日照适宜度（新线性模型）private double calculateSunshineScore(double sunshineRatio) {return Math.min(sunshineRatio / idealSunshineRatio, 1.0);}/*** 综合适宜度计算** @param meanTemp      旬平均温度(℃)* @param precipitation 旬累计降水(mm)* @param sunshineRatio 旬日照比率(实际日照/理论最大日照)* @return 0-1标准化适宜度*/public double evaluate(double meanTemp, double precipitation, double sunshineRatio) {double tScore = calculateTempScore(meanTemp);double pScore = calculatePrecipitationScore(precipitation);double sScore = calculateSunshineScore(sunshineRatio);// 加权求和（权重自动归一化）double totalWeight = tempWeight + precWeight + sunWeight;return (tScore * tempWeight + pScore * precWeight + sScore * sunWeight) / totalWeight;}public static void main(String[] args) {// 配置小麦抽穗期参数SimpleCropSuitability wheatModel = new SimpleCropSuitability(10, 22, 35,    // 温度三基点(℃)60,            // 旬理想降水(mm)0.65,          // 理想日照比率0.4, 0.4, 0.2  // 温度、降水、日照权重);// 测试不同场景testCase(wheatModel, 22, 60, 0.65);  // 理想条件testCase(wheatModel, 15, 30, 0.5);   // 各项不足testCase(wheatModel, 25, 80, 0.7);   // 降水日照超额}private static void testCase(SimpleCropSuitability model,double temp, double prec, double sun) {double suitability = model.evaluate(temp, prec, sun);System.out.printf("温度=%.1f℃ 降水=%.1fmm 日照=%.2f → 适宜度=%.2f\n",temp, prec, sun, suitability);}

这里的降水和日照都是只有适宜值，如果也类似气温一样有最高最低气温和适宜值，那么按照气温修改一下即可使用 

2. 产量分离

根据历年的产量计算出趋势产量和气象产量

    public enum SmoothingOption {THREE_YEAR(3), FIVE_YEAR(5);public final int windowSize;SmoothingOption(int size) { this.windowSize = size; }}public static class YieldData {private final int year;private final double yield;public YieldData(int year, double yield) {this.year = year;this.yield = yield;}public int getYear() { return year; }public double getYield() { return yield; }}public static class DecomposedYield {private final int year;private final Double socialYield;private final Double trendYield;private final Double weatherYield;private final Double relativeYield;public DecomposedYield(int year, Double social, Double trend,Double weather, Double relative) {this.year = year;this.socialYield = social;this.trendYield = trend;this.weatherYield = weather;this.relativeYield = relative;}@Overridepublic String toString() {return String.format("%d\t%s\t%s\t%s\t%s",year,socialYield != null ? String.format("%.2f", socialYield) : "",trendYield != null ? String.format("%.2f", trendYield) : "",weatherYield != null ? String.format("%.2f", weatherYield) : "",relativeYield != null ? String.format("%.1f%%", relativeYield) : "");}}public static void main(String[] args) {// 示例数据（2011-2020）List<YieldData> yields = Arrays.asList(new YieldData(2011, 5.2), new YieldData(2012, 5.5),new YieldData(2013, 5.1), new YieldData(2014, 5.8),new YieldData(2015, 6.0), new YieldData(2016, 5.7),new YieldData(2017, 6.2), new YieldData(2018, 6.5),new YieldData(2019, 6.3), new YieldData(2020, 6.4));// 预报年份设置Integer forecastYear = null;// 3年滑动窗口分析System.out.println("=== 3年滑动窗口分析 ===");process(yields, SmoothingOption.THREE_YEAR, forecastYear);// 5年滑动窗口分析System.out.println("\n=== 5年滑动窗口分析 ===");
//        process(yields, SmoothingOption.FIVE_YEAR, forecastYear);}private static void process(List<YieldData> yields, SmoothingOption option, Integer forecastYear) {// 1. 产量分解List<DecomposedYield> decomposition = decomposeYield(yields, option);// 2. 获取建模有效数据（排除窗口期不足的年份）List<YieldData> modelData = yields.subList(option.windowSize - 1, yields.size());// 3. 打印趋势产量公式printTrendFormula(modelData);// 4. 预报指定年份if (forecastYear != null){double forecast = calculateTrendForecast(modelData, forecastYear);System.out.printf("预报 %d 年趋势产量: %.2f\n", forecastYear, forecast);}// 5. 打印分解结果printDecompositionResults(decomposition);}private static List<DecomposedYield> decomposeYield(List<YieldData> yields, SmoothingOption option) {List<DecomposedYield> results = new ArrayList<>();int windowSize = option.windowSize;// 计算趋势产量（仅使用有效数据）List<YieldData> modelData = yields.subList(windowSize - 1, yields.size());Map<Integer, Double> trendYields = calculateTrendYield(modelData);for (int i = 0; i < yields.size(); i++) {YieldData current = yields.get(i);Double social = (i >= windowSize - 1) ?calculateMovingAverage(yields, i, windowSize) : null;if (social == null) {results.add(new DecomposedYield(current.year, current.yield, null, null, null));} else {Double trend = trendYields.get(current.year);Double weather = trend != null ? current.yield - social : null;Double relative = weather != null ? (weather / social) * 100 : null;results.add(new DecomposedYield(current.year, social, trend, weather, relative));}}return results;}private static double calculateMovingAverage(List<YieldData> yields, int endIndex, int windowSize) {return yields.subList(endIndex - windowSize + 1, endIndex + 1).stream().mapToDouble(YieldData::getYield).average().orElse(0);}private static Map<Integer, Double> calculateTrendYield(List<YieldData> yields) {double[] coeff = calculateRegressionCoefficients(yields);int baseYear = yields.get(0).getYear();return yields.stream().collect(Collectors.toMap(YieldData::getYear,yd -> coeff[0] + coeff[1] * (yd.getYear() - baseYear)));}private static double[] calculateRegressionCoefficients(List<YieldData> yields) {int n = yields.size();int baseYear = yields.get(0).getYear();double sumX = 0, sumY = 0, sumXY = 0, sumX2 = 0;for (YieldData yd : yields) {double x = yd.getYear() - baseYear;double y = yd.getYield();sumX += x;sumY += y;sumXY += x * y;sumX2 += x * x;}double b = (n * sumXY - sumX * sumY) / (n * sumX2 - sumX * sumX);double a = (sumY - b * sumX) / n;return new double[]{a, b};}private static double calculateTrendForecast(List<YieldData> yields, int year) {double[] coeff = calculateRegressionCoefficients(yields);return coeff[0] + coeff[1] * (year - yields.get(0).getYear());}private static void printTrendFormula(List<YieldData> yields) {double[] coeff = calculateRegressionCoefficients(yields);System.out.printf("趋势产量公式: y = %.4f + %.4f * (年份 - %d)\n",coeff[0], coeff[1], yields.get(0).getYear());}private static void printDecompositionResults(List<DecomposedYield> results) {System.out.println("\n产量分解结果:");System.out.println("年份\t社会产量\t趋势产量\t气象产量\t相对产量");results.forEach(System.out::println);}

一般都会进行滑动操作如果不需要删除相关代码即可

3. 基于适宜度计算产量预报

/*** 基于适宜度的产量预报，根据适宜度和产量，计算一元一次回归方程，并计算其相关系数和Significance F值* @param args*/public static void main(String[] args) {// 历史数据：适宜度(X)和产量(Y)double[] suitability = {0.714202, 1.2389059999999996, 0.9274519999999999, 1.175564, 1.150908, 0.9320140000000001, 1.207076, 1.012762, 1.321816};double[] yield = {272.288, 119.1564, 27.4286, 217.9705, 58.2261, 9.9326, -83.3608, -101.6542, -115.4477};// 创建并拟合回归模型SimpleRegression regression = new SimpleRegression();for (int i = 0; i < suitability.length; i++) {regression.addData(suitability[i], yield[i]);}// 获取回归系数double intercept = regression.getIntercept();double slope = regression.getSlope();// 计算F统计量double ssr = regression.getRegressionSumSquares();double sse = regression.getSumSquaredErrors();int n = suitability.length;int k = 1;double fStatistic = (ssr / k) / (sse / (n - k - 1));// 计算Significance F (p值)FDistribution fDist = new FDistribution(k, n - k - 1);double significanceF = 1 - fDist.cumulativeProbability(fStatistic);// 设置4位小数格式DecimalFormat df = new DecimalFormat("0.0000");// 处理极小的p值String formattedPValue;if (significanceF < 0.0001) {formattedPValue = "<0.0001";} else {formattedPValue = df.format(significanceF);}// 输出结果System.out.println("回归方程: Y = " + slope + " * X + " + intercept);System.out.println("相关系数(R): " + regression.getR());System.out.println("决定系数(R²): " + regression.getRSquare());System.out.println("F统计量: " + fStatistic);System.out.println("Significance F (p值): " + formattedPValue);// 模型显著性判断double alpha = 0.05;if (significanceF < alpha) {System.out.println("结论: 回归模型显著 (p < " + alpha + ")");} else {System.out.println("结论: 回归模型不显著 (p ≥ " + alpha + ")");}}

根据方法三我们可以拿到气象产量公式，结合方法二的趋势产量公式代入年份进行计算即可计算出基于适宜度的产量预报