KrigingInterpolation.java

package heyingzhe;

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

public class KrigingInterpolation {
    public static boolean areAllZeros(double[][] data) {
        for (int i = 0; i < data.length; i++) {
            for (int j = 0; j < data[i].length; j++) {
                if (data[i][j] != 0) {
                    return false;
                }
            }
        }
        return true;
    }
    /**
     * 克里金插值方法
     *
     * @param inputArray 原始二维数组，0表示未知值，其他值为已知值
     * @param sill 半变异函数参数：平坦值
     * @param range 半变异函数参数：影响范围
     * @param nugget 半变异函数参数：起点偏移
     * @return 插值后的二维数组
     */
    public static double[][] krigingInterpolation(double[][] inputArray, double sill, double range, double nugget) {

        // 如果没有有效站点数据值，就会返回空数组
        if (areAllZeros(inputArray)) {
            return inputArray;
        }

        int rows = inputArray.length;
        int cols = inputArray[0].length;

        // 自动提取已知点的坐标和值
        List<double[]> samplePointsList = new ArrayList<>();
        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < cols; j++) {
                if (inputArray[i][j] != 0) { // 非0点为已知点
                    samplePointsList.add(new double[]{i, j, inputArray[i][j]});
                }
            }
        }

        // 将List转换为二维数组
        double[][] samplePoints = samplePointsList.toArray(new double[0][]);

        // 深拷贝原始数组，防止修改原始数据
        double[][] result = Arrays.stream(inputArray).map(double[]::clone).toArray(double[][]::new);

        // 对每个未知点（值为0）进行克里金插值
        for (int i = 0; i < rows; i++) {
            for (int j = 0; j < cols; j++) {
                if (result[i][j] == 0) { // 仅对未知点进行插值
                    double[] weights = calculateWeights(i, j, samplePoints, sill, range, nugget);
                    double interpolatedValue = 0.0;

                    // 根据权重和已知点值计算插值
                    for (int k = 0; k < samplePoints.length; k++) {
                        interpolatedValue += weights[k] * samplePoints[k][2];
                    }

                    result[i][j] = interpolatedValue;
                }
            }
        }

        return result;
    }

    // 计算插值点的权重
    private static double[] calculateWeights(int x, int y, double[][] samplePoints,
                                             double sill, double range, double nugget) {
        int n = samplePoints.length;
        double[][] matrix = new double[n + 1][n + 1];
        double[] vector = new double[n + 1];

        // 构建克里金矩阵
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                double h = distance(samplePoints[i][0], samplePoints[i][1], samplePoints[j][0], samplePoints[j][1]);
                matrix[i][j] = variogram(h, sill, range, nugget);
            }
            matrix[i][n] = 1.0;
            matrix[n][i] = 1.0;

            double h = distance(samplePoints[i][0], samplePoints[i][1], x, y);
            vector[i] = variogram(h, sill, range, nugget);
        }

        matrix[n][n] = 0.0;
        vector[n] = 1.0;

        // 求解线性方程组得到权重
        return solveLinearSystem(matrix, vector);
    }

    // 半变异函数模型（例如球状模型）
    private static double variogram(double h, double sill, double range, double nugget) {
        if (h > range) {
            return sill + nugget;
        } else {
            return nugget + sill * (1.5 * (h / range) - 0.5 * Math.pow(h / range, 3));
        }
    }

    // 计算两点之间的欧几里得距离
    private static double distance(double x1, double y1, double x2, double y2) {
        return Math.sqrt(Math.pow(x2 - x1, 2) + Math.pow(y2 - y1, 2));
    }

    // 高斯消去法求解线性方程组
    private static double[] solveLinearSystem(double[][] matrix, double[] vector) {
        int n = vector.length;

        for (int i = 0; i < n; i++) {
            // 寻找主元
            int max = i;
            for (int j = i + 1; j < n; j++) {
                if (Math.abs(matrix[j][i]) > Math.abs(matrix[max][i])) {
                    max = j;
                }
            }

            // 交换行
            double[] temp = matrix[i];
            matrix[i] = matrix[max];
            matrix[max] = temp;

            double t = vector[i];
            vector[i] = vector[max];
            vector[max] = t;

            // 归一化主对角线
            for (int j = i + 1; j < n; j++) {
                double factor = matrix[j][i] / matrix[i][i];
                vector[j] -= factor * vector[i];
                for (int k = i; k < n; k++) {
                    matrix[j][k] -= factor * matrix[i][k];
                }
            }
        }

        // 回代
        double[] solution = new double[n];
        for (int i = n - 1; i >= 0; i--) {
            double sum = 0.0;
            for (int j = i + 1; j < n; j++) {
                sum += matrix[i][j] * solution[j];
            }
            solution[i] = (vector[i] - sum) / matrix[i][i];
        }

        return solution;
    }

    public static void main(String[] args) {
        // 示例：初始化一个2000x1000的二维数组，部分位置设置为已知值
        double[][] inputArray = new double[2000][1000];
        inputArray[100][200] = 5.0;
        inputArray[500][600] = 10.0;
        inputArray[1000][200] = 15.0;

        // 设置半变异函数参数
        double sill = 50.0;   // 数据变化的方差
        double range = 900.0; // 空间相关性消失的距离
        double nugget = 0.1;  // 小尺度变异

        // 调用插值方法
        double[][] interpolatedArray = krigingInterpolation(inputArray, sill, range, nugget);

        // 打印部分结果（仅打印前10行）
        for (int i = 0; i < 10; i++) {
            System.out.println(Arrays.toString(interpolatedArray[i]));
        }
    }
}