更新:
日期 | 内容 |
---|---|
2022.09.13 | 代码中 Vg_Holland2 = r * abs(f) / 2; 中 r 的单位是 km ,需要乘以 1000 转化为 m 。 |
仅使用Matlab实现Holland圆对称风场,不考虑热带气旋的移动速度和前进方位角。
主要的参数设定:
参数 | 公式 |
---|---|
空气密度 | ρ = 1.15 k g ⋅ m − 3 \rho=1.15 kg·m^{-3} ρ=1.15kg⋅m−3 |
外围环境气压 | P n = 1010.0 h P a P_n = 1010.0 hPa Pn=1010.0hPa |
气压差 | Δ P = P n − P c \Delta P = P_n - P_c ΔP=Pn−Pc |
科氏力参数 | f = 2 ∗ 7.292 ∗ 0.00001 ∗ s i n ( L a t c ∗ π / 180.0 ) f = 2 * 7.292 * 0.00001 * sin(Lat_c * \pi / 180.0) f=2∗7.292∗0.00001∗sin(Latc∗π/180.0) |
最大风速半径 | R m = − 18.18 ∗ l o g ( Δ P ) + 112.2 R_m = -18.18 * log(\Delta P) + 112.2 Rm=−18.18∗log(ΔP)+112.2 |
Holland B系数 | B = 1.881 − 0.00557 ∗ R m − 0.01295 ∗ ∣ L a t c ∣ ) B = 1.881 - 0.00557 * R_m - 0.01295 * |Lat_c|) B=1.881−0.00557∗Rm−0.01295∗∣Latc∣) |
气压场 | P r = P c + Δ P ∗ e x p ( − ( R m / r ) B ) P_r = P_c + \Delta P * exp(-(R_m / r) ^ B) Pr=Pc+ΔP∗exp(−(Rm/r)B) |
梯度风场 | V g ( r ) = Δ P B ρ ( R m r ) B e x p ( − ( R m r ) B ) ) + ( r ⋅ f 2 ) 2 + r ⋅ ∣ f ∣ 2 V_g(r)=\sqrt{\Delta P \frac{B}{\rho}(\frac{R_m}{r})^Bexp(-(\frac{R_m}{r})^B))+(\frac{r·f}{2})^2} + \frac{r·|f|}{2} Vg(r)=ΔPρB(rRm)Bexp(−(rRm)B))+(2r⋅f)2+2r⋅∣f∣ |
近地表10m处的平均风速 | V 10 m ( r ) = K m ∗ V g ( r ) V_{10m}(r) = K_m * V_g(r) V10m(r)=Km∗Vg(r) |
K m K_m Km取值标准 | K m = { 0.81 , V g < 6 0.81 − 2.96 ∗ 1 0 − 3 ∗ V g , 6 ≤ V g < 19.5 0.77 − 4.31 ∗ 1 0 − 3 ∗ V g , 19.5 ≤ V g < 45 0.66 , V g ≥ 45 K_m=\begin{cases}0.81, V_g < 6 \\0.81 - 2.96 * 10 ^ {-3} * V_g, 6\leq V_g<19.5 \\0.77 - 4.31 * 10 ^ {-3} * V_g, 19.5\leq V_g<45\\0.66, V_g \ge 45\end{cases} Km=⎩ ⎨ ⎧0.81,Vg<60.81−2.96∗10−3∗Vg,6≤Vg<19.50.77−4.31∗10−3∗Vg,19.5≤Vg<450.66,Vg≥45 |
生成0.25°分辨率的风场:
附Matlab代码:
其中,basetif.tif是一个0.25°的影像
中心经纬度以及中心气压:
P
c
=
995
h
P
a
P_c = 995 hPa
Pc=995hPa
L
a
t
c
=
10.125
°
Lat_c = 10.125 °
Latc=10.125°
L
o
n
c
=
120.125
°
Lon_c = 120.125 °
Lonc=120.125°
%%
clc;
clear;
%%
fullpath = mfilename('fullpath');
[path, name] = fileparts(fullpath);
%%
%
basetif = strcat(path, '\basetif.tif');
[A_M, R_M] = geotiffread(basetif);
info = geotiffinfo(basetif);
A_I = zeros(size(A_M));
A_J = zeros(size(A_M));
for i = 1 : size(A_I, 1)
A_I(i, :) = i;
end
for j = 1 : size(A_J, 2)
A_J(:, j) = j;
end
[A_Lat, A_Lon] = pix2latlon(info.RefMatrix, A_I, A_J);
%%
P_c = 995; % hPa
Lat_c = 10.125; % °
Lon_c = 120.125; % °
%%
%
rou_air = 1.15; % kg m-3
%
P_n = 1010.0;
dertP = P_n - P_c;
Rm = -18.18 * log(dertP) + 112.2; % R2 = 0.84
B = 1.881 - 0.00557 * Rm - 0.01295 * abs(Lat_c);
f = 2 * 7.292 * 0.00001 * sin(Lat_c * pi / 180.0);
% great cricle arc
R_earth = 6371;
arcLAT1 = A_Lat * pi / 180.0;
arcLAT2 = Lat_c * pi / 180.0;
arcLON1 = A_Lon * pi / 180.0;
arcLON2 = Lon_c * pi / 180.0;
cos_sita = sin(arcLAT1) .* sin(arcLAT2) + ...
cos(arcLAT1) .* cos(arcLAT2) .* cos(arcLON1 - arcLON2);
sita = acos(cos_sita);
r = R_earth * sita + 0.001;
P_r = P_c + dertP * exp(-(Rm ./ r) .^ B);
Vg_Holland1 = 100 * dertP * B / rou_air * (Rm ./ r) .^ B .* exp(-(Rm ./ r) .^ B);
Vg_Holland2 = r * abs(f) / 2 * 1000;
Vg_Holland = (Vg_Holland1 + Vg_Holland2 .^ 2) .^ 0.5 - Vg_Holland2;
%%
Vg = Vg_Holland;
Km = zeros(size(Vg_Holland));
Km(Vg < 6) = 0.81;
Km(Vg >= 6 & Vg < 19.5) = 0.81 - 2.96 * 10 ^ (-3) * (Vg(Vg >= 6 & Vg < 19.5) - 6);
Km(Vg >= 19.5 & Vg < 45) = 0.77 - 4.31 * 10 ^ (-3) * (Vg(Vg >= 19.5 & Vg < 45) - 19.5);
Km(Vg >= 45) = 0.66;
V_10m = Km .* Vg;
outfile = strcat(path, '\V_10m.tif');
geotiffwrite(outfile, single(V_10m), R_M);