MATLAB算法判别热带气旋

Question

将六条要求合并到一起，做成一个条件嵌套的判别式，能够一层一层往下，符合第一个要求，就判断下一个。最好附带一个台风实例（任意都行），验证该台风是否符合要求，能否通过判别。

% 1.计算相对涡度
max_vorticity_850 = max(vorticity_850, [], 'all') % 计算相对涡度
critical_vorticity = 8*10^-5;
if max_vorticity_850 < critical_vorticity
is_typhoon = false;
end
 
% 2.计算850hPa的最大风速
max_wind_speed_850 = max(wind_speed_850, [], 'all'); % 计算最大风速
max_wind_speed_300 = max(wind_speed_300, [], 'all'); % 计算最大风速
critical_wind_speed = 13;
if max_wind_speed_850 < critical_wind_speed
is_typhoon = false;
end

% 3.计算暖芯
[~, max_vorticity_idx] = max(vorticity, [], 'all', 'linear'); % 最大相对涡度的索引
[lat_idx, lon_idx] = ind2sub(size(vorticity), max_vorticity_idx); % 最大相对涡度的经纬度索引
lat_box = lat(lat_idx-5:lat_idx+4, lon_idx-5:lon_idx+4); % 选取10°x 10°的经纬度格网
lon_box = lon(lat_idx-5:lat_idx+4, lon_idx-5:lon_idx+4);
t300_box = t300(lat_idx-5:lat_idx+4, lon_idx-5:lon_idx+4); % 选取相应区域内的300,500和700 hPa的温度
t500_box = t500(lat_idx-5:lat_idx+4, lon_idx-5:lon_idx+4);
t700_box = t700(lat_idx-5:lat_idx+4, lon_idx-5:lon_idx+4);
mean_t300_box = mean(t300_box, 'all'); % 计算均值
mean_t500_box = mean(t500_box, 'all');
mean_t700_box = mean(t700_box, 'all');
max_t_box = max([t300_box(:); t500_box(:); t700_box(:)]); % 计算最大温度
temp_deviation = [max_t_box - mean_t300_box, max_t_box - mean_t500_box, max_t_box - mean_t700_box]; % 计算温度偏差
sum_temp_deviation = sum(temp_deviation);
critical_deviation = 0.8;
if sum_temp_deviation < critical_deviation
is_typhoon = false;
end
 
% 4.比较850hPa和300hPa的最大风速
max_wind_speed_300 = max(wind_speed_300, [], 'all'); % 计算最大风速
if max_wind_speed_850 < max_wind_speed_300
is_typhoon = false;
end

% 5.在北印度洋应用半径阈值
if in_nio(lon_box, lat_box)
radius_threshold = 200e3; % 200 km 半径阈值
max_wind_speed_location = max(wind_speed_850, [], 'all', 'linear'); % 最大风速的经纬度索引
[storm_lat_idx, storm_lon_idx] = ind2sub(size(wind_speed_850), max_wind_speed_location); % 最大风速的经纬度
[storm_lat, storm_lon] = deal(lat(storm_lat_idx), lon(storm_lon_idx)); % 最大风速的经纬度值
% 计算最大平均风速的半径
[R, ~, ~] = distance(storm_lat, storm_lon, lat_box, lon_box, referenceSphere('Earth', 'meters'));
radius = mean(R(wind_speed_850 > 0.8*max_wind_speed_850));
% 如果最大平均风速的半径大于半径阈值，则认为不是台风
if radius > radius_threshold
is_typhoon = false;
end
end
 
% 6.计算台风的持续时间
time_threshold = 36*3600; % 36小时
typhoon_duration = time(end) - time(1); % 台风持续时间
if typhoon_duration < time_threshold
is_typhoon = false;
end

Answer 1

clc
close all
clera all
% 1.计算相对涡度
max_vorticity_850 = max(vorticity_850, [], 'all'); % 计算相对涡度
critical_vorticity = 8*10^-5;
if max_vorticity_850 < critical_vorticity
    is_typhoon = false;
else
    % 2.计算850hPa的最大风速
    max_wind_speed_850 = max(wind_speed_850, [], 'all'); % 计算最大风速
    max_wind_speed_300 = max(wind_speed_300, [], 'all'); % 计算最大风速
    critical_wind_speed = 13;
    if max_wind_speed_850 < critical_wind_speed
        is_typhoon = false;
    else
        
        % 3.计算暖芯
        [~, max_vorticity_idx] = max(vorticity, [], 'all', 'linear'); % 最大相对涡度的索引
        [lat_idx, lon_idx] = ind2sub(size(vorticity), max_vorticity_idx); % 最大相对涡度的经纬度索引
        lat_box = lat(lat_idx-5:lat_idx+4, lon_idx-5:lon_idx+4); % 选取10°x 10°的经纬度格网
        lon_box = lon(lat_idx-5:lat_idx+4, lon_idx-5:lon_idx+4);
        t300_box = t300(lat_idx-5:lat_idx+4, lon_idx-5:lon_idx+4); % 选取相应区域内的300,500和700 hPa的温度
        t500_box = t500(lat_idx-5:lat_idx+4, lon_idx-5:lon_idx+4);
        t700_box = t700(lat_idx-5:lat_idx+4, lon_idx-5:lon_idx+4);
        mean_t300_box = mean(t300_box, 'all'); % 计算均值
        mean_t500_box = mean(t500_box, 'all');
        mean_t700_box = mean(t700_box, 'all');
        max_t_box = max([t300_box(:); t500_box(:); t700_box(:)]); % 计算最大温度
        temp_deviation = [max_t_box - mean_t300_box, max_t_box - mean_t500_box, max_t_box - mean_t700_box]; % 计算温度偏差
        sum_temp_deviation = sum(temp_deviation);
        critical_deviation = 0.8;
        if sum_temp_deviation < critical_deviation
            is_typhoon = false;
        else
            
            % 4.比较850hPa和300hPa的最大风速
            max_wind_speed_300 = max(wind_speed_300, [], 'all'); % 计算最大风速
            if max_wind_speed_850 < max_wind_speed_300
                is_typhoon = false;
            else
                
                % 5.在北印度洋应用半径阈值
                if in_nio(lon_box, lat_box)
                    radius_threshold = 200e3; % 200 km 半径阈值
                    max_wind_speed_location = max(wind_speed_850, [], 'all', 'linear'); % 最大风速的经纬度索引
                    [storm_lat_idx, storm_lon_idx] = ind2sub(size(wind_speed_850), max_wind_speed_location); % 最大风速的经纬度
                    [storm_lat, storm_lon] = deal(lat(storm_lat_idx), lon(storm_lon_idx)); % 最大风速的经纬度值
                    % 计算最大平均风速的半径
                    [R, ~, ~] = distance(storm_lat, storm_lon, lat_box, lon_box, referenceSphere('Earth', 'meters'));
                    radius = mean(R(wind_speed_850 > 0.8*max_wind_speed_850));
                    % 如果最大平均风速的半径大于半径阈值，则认为不是台风
                    if radius > radius_threshold
                        is_typhoon = false;
                    else
                        % 6.计算台风的持续时间
                        time_threshold = 36*3600; % 36小时
                        typhoon_duration = time(end) - time(1); % 台风持续时间
                        if typhoon_duration < time_threshold
                            is_typhoon = false;
                        end
                    end
                else
                    % 6.计算台风的持续时间
                    time_threshold = 36*3600; % 36小时
                    typhoon_duration = time(end) - time(1); % 台风持续时间
                    if typhoon_duration < time_threshold
                        is_typhoon = false;
                    end
                end
            end
        end
    end
end