It is crucial to accurately reproduce the climatological features of the thermospheric neutral winds. We used the 2020–2022 data of the Michelson Interferometer for Global High-Resolution Thermospheric Imaging (MIGHTI) on the Ionospheric Connection Explorer (ICON) to develop an empirical model of thermospheric winds with the methods of the Non-Uniform Rational B-Spline and harmonic fittings (NURBS-Harmonic model). The NURBS-Harmonic model exhibits good adaptability for winds in every season with longitude, latitude, local time, and altitude. Over 70 % of the NURBS-Harmonic data are present with errors within 15 m/s, showing a better performance than the Horizontal Wind Model (HWM14). The percentage of HMW14 data errors within 15 m/s is about 40 %. The NURBS-Harmonic effectively captures the spatial structure of the wind fields and exhibits seasonal variations well. At 250 km, zonal winds in the low and middle latitudes at 12 LT exhibit a wavenumber-3 (WN3) structure in all the seasons, while they display a slight wavenumber-4 (WN4) longitude structure trend during June Solstice and September Equinox. Meridional winds at the same altitude show a distinct WN4 structure in every season. Different from the HWM14, which only includes migrating tides, the NURBS-Harmonic model also considers non-migrating tidal components. To validate the significance of non-migrating tides, we conducted a wave analysis on meridional winds at 250 km in September Equinox near the equator where non-migrating tides are relatively important. The analysis indicates that the amplitudes of the predominant non-migrating tidal components (SE2, DW2, and D0) are nearly equal and about 40 % of the dominant migrating tidal component DW1. This suggests that the non-migrating tidal components can contribute significantly to the overall wind field variability and should not be overlooked. In addition, the extensive height coverage from 91 km to 300 km of the MIGHTI data enables the model to provide precise altitude information. The zonal WN4 structure of zonal winds is most prominent in the altitude ranges of 100–250 km at low latitudes in the Northern Hemisphere during June Solstice and September Equinox. The model accurately reflects the phases variations of the winds at various altitudes. The phases of WN3 and WN4 at 110 km are opposite to those at 150 km, existing in all seasons. This may partially explain previous reports that the topside ionosphere equatorial vertical plasma drift is positively correlated with the zonal winds at 110 km and inversely correlated with those at 150 km.