基于Emanuel(1986)热机理论，我们提出了一种低聚方法以海平面变压等实测数据反演热带气旋(TC)风场结构。本贴旨介绍历史研究背景及这个方法的实现思路，对于针眼或实测密度稀缺等不确定的情况，建议谨慎使用方法给出的风速结论。
Emanuel(1986)理论的详细介绍请参考我们在新风吧的帖子：【闲聊吹水/理论分析】关于热带气旋垂直...
我们暂定的分析目标桦1镇楼
萌新第一次在风吧开分析贴，望85莫删

膜

我们先来介绍并回顾一下传统风压对应模型。
在热带气旋地面观测数据中，由于下垫面摩擦、稀疏的时间或空间观测密度，导致地面观测较难录得热带气旋的最大平均风速。然而，气压观测几乎排除了任何外界干扰，可以极其准确地作为飓风中的物理量。因此长期以来，热带气旋的风速-气压对应关系一直是备受关注的目标。
在过去数十年里，人们提出了多种风-压关系，譬如美国国家飓风中心（NHC）和中太平洋飓风中心(CPHC)所使用的Dvorak大西洋对应方法。在2009年Courtney、Knaff和Zehr的研究中，确立了一种标准的风-压关系，目前已广泛适用于全球各大热带气旋预警中心（如JTWC），这就是我们经常使用的KZC风压对应。这里给出其表达式：
其中V_srm1是热带气旋处于准静止状态下10米1分钟最大平均风速，以下是其与经过风暴移动速度（C）调整后的10米1分钟最大平均风速（V_max1）之间的关系：
S是位于400~600km半径处的实际切向风速（V500）与500km半径处的气候态平均切向风速（V500_c）之比，而V500又可以表示为烈风半径（R34）的函数：
然而，KZC方法只是一种统计学关系。对于RMW非常小的热带气旋而言，比如飓风Wilma（2005），从理论和统计关系两个角度均呈现对于这类系统的Δp低估现象。他们在论文中表示，如果将RMW作为回归拟合参数，可能会取得更好的效果。

因此，本文提出了一种方法，用来实现海平面气压观测数据反演近海表1分钟最大平均风速，该方法基于Emanuel于1986年提出的海气耦合模型（以下简称E86）。在这一方法中，采用一个热带气旋过境时的海平面变压数据，并使用一些可以简单获取的物理量，便可以计算中心附近的最大平均风速。这一方法对现有的中心海平面气压-最大10米风速的统计学对应关系（以下简称风-压关系）进行了改进。
此方法可以应用至大部分热带气旋观测数据中，并且相比于传统风-压关系，此方法不受最大风速半径（RMW）大小的影响。此外，该算法还可以应用于热带气旋结构分析，如眼墙置换（ERC）过程中的风场结构。
另外需要注意的是，本文提出的方法对定位的准确度要求相对较高，如果有雷达辅助的话可能会得到更精确的结果。
在我们写的低聚程序以及吧内大佬的研究成果中，可以使用雷达测风对中心气压进行估算（详见我们的帖子）。我们可以通过雷达扫描的径向速度进行CAPPI处理，并用径向速度除以径向与风向的夹角的余弦值，得到TC给定高度（或气压）面上的切向风速，将其视为梯度风后代入梯度风平衡方程
静力平衡方程
及状态方程
去计算海平面气压，这就是JMA相关论文中使用的DR法。由于TC具有暖心结构，如果不对分析高度上的暖心距平进行处理，则会造成分析层面的降压高估。反过来，我们提出一种猜想，是否可以对海平面的变压实测通过类似“反向DR”的方式反演海平面风速？然而这并不容易解决，边界层以下由于下垫面摩擦显然是次梯度风，想要实现这一点必须先通过海平面的气压分布反推700mb的高度分布，用此计算FL梯度风，转化为10sFL后乘以0.9折算系数计算海平面实际1min平均风。我们对此思路进行了尝试，首先手动输入FL气温，对700mb至海平面气压进行气温线性插值，并使用静力平衡+理想气体定律计算700mb位势高度，使用位势高度的梯度风平衡式计算FL风速，并使用前面提到的方法得到海平面风速。然而效果并不理想，手动输入的暖心距平容易使得FL梯度风出现不严谨且不合理的波动。因此我们尝试提出一种相对而言更严谨的方式去实现这一目的。

历史上，美国国家海洋和大气管理局（NOAA）进行过数十次强烈飓风的飞机侦察。对于强飓风而言，飞行层面通常为700mb或750mb。Stern和Nolan（2011，以下简称SN11）的研究指出，部分热带气旋在此层面呈现出非平衡气流。WRF模拟显示，对于拥有较小的RMW的热带气旋，飞行层面倾向于出现一股次梯度气流，对应4-5km高空的一股超梯度气流。SN11指出，正是因为这种垂直结构，导致部分飓风的最大切向风的垂直衰减速率明显偏小（包括折算系数偏大）。然而，尽管如此，大部分热带气旋的飞行层面通常符合梯度风平衡。
飞行层面风速的确定是关键的。Franklin（2003）基于飞机投落式GPS探空仪的数据，得到了热带气旋的一条归一化的平均垂直风速廓线。他们的研究指出，WS10m/WS700mb通常为平均的0.91，这一折算因子已经广泛应用于热带气旋近海表风速的评定。
在热带气旋实际观测数据与ERA5数据集中，可以获取以下变量：海平面气温T_s、近似流出温度T_0、热带气旋外围海平面相对湿度RH_as、背景气压p_a、参考气压p_0=1000mb与热带气旋中心的距离及对应的海平面气压数列
。首先，我们需要得到一条连贯的降压曲线，所以我们放弃了线性插值，对数据(p_i, r_i)s进行三次样条插值，得到p_s作为r的函数，其性质可导可微。

在E86理论中的所有推导基于且仅基于斜压、热成风平衡、气块沿等角动量面上升时呈广义中性浮力状态，在实际观测数据中证明了这些假设的合理性；E86认为热带气旋的低压是由海平面气温T_s和相当位温梯度导致的，反映在相对湿度的径向增长上，且根据边界层物理量的径向分布建模可以得到TC风场结构，（关于该理论详见我们在新台风吧内的帖子）。其中热带气旋中的相当位温作为一个关键量，由于湿绝热过程其在等角动量面上均匀分布。为了方便起见，我们首先尝试计算海平面背景气块混合比q_as。Tetens方程给出了饱和水汽压es(T)的积分近似表达式，并可以用此导出海平面背景饱和混合比q*_as：
根据
可以轻松得到海平面实际气块混合比的值。
再由相当位温的定义
可以求出背景相当位温。
接下来，我们假设热带气旋近海表同样满足斜压湿中性和热成风平衡两个假设，E86中的式(18)应可以重写为：
其中pi_s是由Exner函数(6)定义的。通过解上述方程，可以得到θ_es(r).
显然上式是一个仅关于r作为θ_e的函数，其中的变量和常量均可以根据已知和先前的推导得出，通过求解该方程可以解出整个热带气旋近海表的相当位温θ_es(r)。
值得注意的是，这里我们可以把T_out近似为T_0进行处理。然而在实际热带气旋中的流出温度并不完全为定值。比如可以通过下式进行修改：
其中的
, 近似值为2.0。
然而其变异性对海平面的影响相当小，因此我们在这里并没有采用相关修正。

同样根据相当位温的定义，由Exner函数的定义和对数定理，有：
那么
我们可以通过上式去计算海平面混合比q_s作为r的函数，其适用于整个热带气旋。上面的推导先计算相当位温分布，后计算比湿分布，实际上和Emanuel的处理方式一致。
下一步，我们进行计算海平面的梯度风场。湿空气的状态方程可以写为：
其中T_v是虚温，α为比容，其定义为1/ρ。
等高梯度风方程如下式所示：
我们根据上式，最终得出整个热带气旋中的海平面梯度风场。
在实际热带气旋边界层底部，由于粘性以及下垫面摩擦，其实际风速明显低于梯度风。然而尽管如此，由于我们所求的飞行层面风速遵循梯度风平衡原理，解出海平面梯度风的过程显然是必要的。另外，虽然上述讨论可能从理论上并不完全适用于眼区，然而眼区的下沉绝热增温很可能可以基本抵消相对湿度的高估。因此，我们认为使用上述模型解出的梯度风可以基本适用于整个热带气旋，对于进行眼墙置换(ERC)过程的热带气旋，该方法仍可以较好的捕捉到次级眼墙产生的第二峰风速峰值。
我们目前构建出了一种基于海平面实际气压观测的梯度风模型。相比于经验模型，譬如Holland，它同样给出了梯度风的预测。然而其对于紧凑的热带气旋(如飓风Tracy)，其无法较好的反应核心区高达5mb/km的强烈气压梯度。Holland的变压模型认为
然而其对于ERC过程中的热带气旋容易失效。并且其没有考虑空气密度对梯度风计算结果的影响，因此可能造成一定程度上的梯度风低估。由于本模型基于实际的气压梯度观测，所以可以有效避免上述问题。
另外需要指出的是，该方法对TC外围的相当位温垂直分布并不一定完全准确（比较类似于平均值），因为TC外围存在水平或垂直湍流混合（即并不完全遵守湿绝热抬升）。然而用本方法算出的梯度风结构等应大致无问题（尤其是核心区）。

接下来我们尝试对飞行层面的梯度风进行处理。我们首先需要得到700mb高度的气温。
由于接下来的推导使用的是饱和混合比（而非实际混合比），假设海平面气温恒定，我们首先根据公式
计算TC内部的海平面饱和混合比，其中外围饱和混合比已经由前面的推导计算得出。
参考Clausius-Clapeyron式，给定气压层面的饱和混合比可以表示为：
又由位温的定义（θ=T/pi），我们分类讨论饱和大气（LCL以上）和非饱和大气中饱和混合比与实际混合比的大小关系，并代入相当位温的定义，有以下方程组：
这个方程组直接给出了指定气压层面M面的温度T。
然而尽管在Python代码编写时可以使用fsolve函数直接求解上述复杂非线性方程组，但是作为理论E86给出了这个方程组的迭代求解方法。我们有必要介绍一下该迭代法：
首先猜测某特定M面上某特定气压点pi对应的T值（或者也可以反过来），根据
Clausius-Clapeyron式求饱和混合比，並根据q_s与q*的大小关系判断是否饱和；
求位温
；
由pi 和 θ求得新的T值： T = θ * pi；
这个T值作为新的估计值，并重复迭代上述过程，几次迭代后T将收敛于迭代次数，便可以
得出精确的M面上T、θe、q*和pi值。
接下来我们计算指定位势高度M面的半径并用其计算梯度风。在E86中，作者给出了M面
半径随气温变化的方程，其表示为：
对于某个给定的位势高度z(p)，使用上述过程中迭代出的值，代入上式即得

最后一步，我们直接求出FL梯度风。对于已求得的海平面的梯度风，根据
求出M面的M值。对于FL，根据前述步骤中求得的
代入
解得
至此，我们得到了一套使用海平面变压数据计算FL梯度风的方法。以Ragasa(2025)的一巅浮标数据为例，尝试使用该模型进行具体的梯度风分析：
I. 输入参数:
观测点数: 10
海平面气温 T_s: 300.00 K (26.9°C)
流出温度 T_0: 200.00 K (-73.1°C)
外围相对湿度 RH_as: 0.800
纬度: 15.00°
风暴外半径 r_0: 500.0 km
目标气压: 700 hPa
背景气压 p_a: 1006.5 hPa
----------------------------------------------------------------------
II. 计算海平面气块混合比
科氏参数 f: 3.77e-05 s^-1
饱和水汽压 es(T_s): 35.35 hPa
外围饱和混合比 q_as*: 22.137 g/kg
外围实际混合比 q_as: 17.710 g/kg
外围相当位温 θ_ea: 346.81 K
热力学参数:
热力学效率 ε: 0.3333
中心点 (r=0) 计算示例:
气压 p: 893.6 hPa
相当位温 θ_es: 382.25 K
实际混合比 q_s: 25.344 g/kg
饱和混合比 q_s*: 24.773 g/kg
----------------------------------------------------------------------
III. 计算海平面梯度风并确定最大风速半径(RMW)
海平面梯度风分析:
最大风速 V_max: 79.65 m/s
最大风速半径 RMW: 13.53 km
中心气压: 893.6 hPa
M面起点 (RMW处):
半径: 13.53 km
风速: 79.65 m/s
气压: 932.8 hPa
Exner函数 π_s: 0.980315
实际混合比 q_s: 18.623 g/kg
饱和混合比 q_s*: 23.733 g/kg
相当位温 θ_e: 357.13 K
----------------------------------------------------------------------
IV. 计算M面上的绝对角动量
M面角动量 M: 1080766 m²/s
----------------------------------------------------------------------
V. 计算700 hPa高度的温度
目标气压: 700 hPa
目标Exner函数 π_target: 0.903141
海平面实际混合比 q_s: 18.623 g/kg
海平面饱和混合比 q_s*: 23.733 g/kg
700 hPa温度 T_700: 286.61 K (13.5°C)
700 hPa混合比: 13.605 g/kg (饱和)
----------------------------------------------------------------------
VI. 计算700 hPa的M面半径
700 hPa M面半径: 14.417 km
----------------------------------------------------------------------
VII. 计算700 hPa梯度风
700 hPa梯度风: 74.69 m/s
对应半径: 14.417 km
温度: 286.61 K (13.5°C)
混合比: 13.605 g/kg
与海平面(RMW)对比:
梯度风速: 74.69 m/s (700 hPa) --- 79.65 m/s (海平面)
半径: 14.417 km (700 hPa) --- 13.525 km (海平面)
实验结果的各项物理量均比较符合预期。尽管根据内眼墙外缘至风眼65mb的气压降幅参考其他实测案例气压梯度的估算问题不大，但还是由于Ragasa的使用IR/VIS中心定位不确定性相对较大（且实测过程中系统内眼正在不断升压减弱），该实验结果在具体定强过程中仅供参考。

目前我们已经求出了FL梯度风速。在飞机实测中，使用30sFL风速作为梯度风分布计算得出的中心气压更符合实际中心气压，因此使用上述方法计算出的700mb梯度风可以基本视为30sFL风速。我们现在需要得到10sFL与30sFL的比值关系，依此计算10sFL风速并使用
~0.9的折算系数得到10m一分钟最大持续风。通过读取HDOB报文，我们将简单对30sFL风速大于50kt的数据进行分析。
由于RMW内的水平风速梯度非常剧烈，导致30sFL风速并不一定是梯度风，且会造成10sFL/30sFL的值异常偏大，因此我们决定手动剔除R<RMW区域的数据。接下来基于最小二乘法拟合10sFL与30sFL的一元线性回归方程，最小二乘法公式为：
β1是回归系数（直线斜率），β0是截距。
拟合优度良好，R^2为0.99，显著性检验p值接近0。因此我们直接使用
作为10sFL与30sFL的关系回归方程。

我们知道，热带气旋的飞行层面风场通常高于Ekman层，其风场是梯度风矢量与移动速度矢量之和。在垂直于移动矢量的危险半圆处，风速在该半径的环上取到最大，其速度等于梯度风速与移速之和。因此将前面求出的700mb梯度风与移速直接相加，并使用10sFL与30sFL的回归方程计算TC中最大10sFL，可以由最大10sFL转换为近海平面1min持续风速：
最后我们使用程序进行100m间隔网格根据以上描述的方法计算对应的700mb半径、气温、梯度风，并使用上式反推海平面持续风。并且也可以用方程（11）计算外围的700hPa温度，依此计算等压面暖心距平（尽管外围并不完全符合湿绝热）：
q*_as和q_as由方程1-3的Tetens公式给出，θ_ea由方程(4)给出，π_as = (pa/p0)^(R/Cp)由Exner函数定义。
以桦1的浮标数据为例，我们计算出的Ragasa第一次巅峰的最大1分钟10米持续风达到147节(还是相对而言比较符合吧内的主流评价)，最大700hPa暖心距平为8.8K.

根据我们的算法可以成功捕捉到约60km处正在形成的次级眼墙，可以看出此时系统处于置换初期，外眼墙尚未对内眼造成显著压制。这一结果比较符合附近扫过的微波表现。

其实如果算近表风的话，摩擦，湍流交换系数等都需要考虑，甚至假设有个测站在10m高处采样了1分钟，其湍流风谱特性可能也会导致出现偏差，梯度风转实际风里面这个部分其实影响还是非常大的，而且还很麻烦（）这也间接导致目前10m高度风体系其实比较的不可靠）