公式摘自 GB/T 33703—2017《自动气象站观测规范》。

公式

计算公式如下：

\(\overline{W_D} = \arctan\left(\frac{\overline{X}}{\overline{Y}}\right)\)

式中：

\(\overline{W_D}\)——观测时段内的平均风向；

\(\overline{X}\)——观测时段内单位矢量在\(x\)轴（西东方向）上的平均分量；

\(\overline{X} = \frac{1}{N} \times \sum_{i=1}^{N} \sin D_i\)

\(\overline{Y}\)——观测时段内单位矢量在\(y\)轴（南北方向）上的平均分量；

\(\overline{Y} = \frac{1}{N} \times \sum_{i=1}^{N} \cos D_i\)

\(D_i\)——观测时段内第\(i\)个风矢量的幅角（与\(y\)轴的夹角）；

\(N\)——观测时段内的样本总数。

根据\(\overline{X}\)、\(\overline{Y}\)的正负，按以下方法对\(\overline{W_D}\)进行修正：

\(\overline{X} > 0\)、\(\overline{Y} > 0\)，\(\overline{W_D}\)无需修正；

\(\overline{X} > 0\)、\(\overline{Y} < 0\)或\(\overline{X} < 0\)、\(\overline{Y} < 0\)，\(\overline{W_D}\)加180°。

\(\overline{X} < 0\)、\(\overline{Y} > 0\)，\(\overline{W_D}\)加360°。

参考代码

Python

import math

# 单位矢量平均法计算平均风向
def calculate_average_wind_direction(wind_directions):
    # 转为弧度制
    wind_directions_rad = [math.radians(d) for d in wind_directions]

    # 计算平均值
    X_avg = sum(math.sin(d) for d in wind_directions_rad) / len(wind_directions)
    Y_avg = sum(math.cos(d) for d in wind_directions_rad) / len(wind_directions)

    # 计算风向平均值的弧度
    average_wind_direction_rad = math.atan2(X_avg, Y_avg)

    # 转回角度制
    average_wind_direction_deg = math.degrees(average_wind_direction_rad)

    # 出界修正
    if average_wind_direction_deg < 0:
        average_wind_direction_deg += 360

    return average_wind_direction_deg

C++

#include <vector>
#include <cmath>
#include <numeric>

double calculate_average_wind_direction(const std::vector<double>& wind_directions) {
    // 转为弧度制
    std::vector<double> wind_directions_rad;
    wind_directions_rad.reserve(wind_directions.size());
    for (double d : wind_directions) {
        wind_directions_rad.push_back(d * M_PI / 180.0);
    }

    // 计算平均值
    double X_avg = std::accumulate(wind_directions_rad.begin(), wind_directions_rad.end(), 0.0, 
        [](double sum, double d) { return sum + std::sin(d); }) / wind_directions.size();

    double Y_avg = std::accumulate(wind_directions_rad.begin(), wind_directions_rad.end(), 0.0, 
        [](double sum, double d) { return sum + std::cos(d); }) / wind_directions.size();

    // 计算风向平均值的弧度
    double average_wind_direction_rad = std::atan2(X_avg, Y_avg);

    // 转回角度制
    double average_wind_direction_deg = average_wind_direction_rad * 180.0 / M_PI;

    // 出界修正
    if (average_wind_direction_deg < 0) {
        average_wind_direction_deg += 360.0;
    }

    return average_wind_direction_deg;
}

C#

using System;
using System.Collections.Generic;
using System.Linq;

public class WindDirectionCalculator
{
    public static double CalculateAverageWindDirection(List<double> windDirections)
    {
        // 转为弧度制
        var windDirectionsRad = windDirections.Select(d => d * Math.PI / 180.0).ToList();

        // 计算平均值
        double X_avg = windDirectionsRad.Select(d => Math.Sin(d)).Average();
        double Y_avg = windDirectionsRad.Select(d => Math.Cos(d)).Average();

        // 计算风向平均值的弧度
        double averageWindDirectionRad = Math.Atan2(X_avg, Y_avg);

        // 转回角度制
        double averageWindDirectionDeg = averageWindDirectionRad * 180.0 / Math.PI;

        // 出界修正
        if (averageWindDirectionDeg < 0)
        {
            averageWindDirectionDeg += 360.0;
        }

        return averageWindDirectionDeg;
    }
}

Java

import java.util.List;

public class WindDirectionCalculator {

    public static double calculateAverageWindDirection(List<Double> windDirections) {
        // 转为弧度制
        double[] windDirectionsRad = windDirections.stream()
                .mapToDouble(d -> Math.toRadians(d))
                .toArray();

        // 计算平均值
        double X_avg = 0.0;
        double Y_avg = 0.0;

        for (double d : windDirectionsRad) {
            X_avg += Math.sin(d);
            Y_avg += Math.cos(d);
        }

        X_avg /= windDirections.size();
        Y_avg /= windDirections.size();

        // 计算风向平均值的弧度
        double averageWindDirectionRad = Math.atan2(X_avg, Y_avg);

        // 转回角度制
        double averageWindDirectionDeg = Math.toDegrees(averageWindDirectionRad);

        // 出界修正
        if (averageWindDirectionDeg < 0) {
            averageWindDirectionDeg += 360.0;
        }

        return averageWindDirectionDeg;
    }
}