如何降低插值轨迹数据集计算二次极化损耗项的偏差？

核心问题 #

我编写了Python程序计算移动发射机与天线的通信接收信号强度，发射机轨迹数据非均匀间隔，需支持用户选择采样率进行插值。但发现采样率变化时，信号强度公式中极化损耗ploss的平方项数值差异极大——采样率30与60时，该平方项最大值偏差可达10倍，推测是插值的微小偏差被平方放大导致。我缺乏数据科学/scipy相关经验，不知道如何平滑数据，也不清楚搜索关键词。

现有核心代码 #

位置插值函数 #

def interp_position(times, sample_rate, position_vector, start_time=0, end_time= None, coord_system = "spherical",new_times=None):
    #times with parameterized sample rate
    if end_time is None:
        end_time = times[-1]
    if new_times is None:
        new_times = np.arange(times[0], times[-1],1/sample_rate)

    position_vector[:, 1] = np.radians(position_vector[:, 1])  
    position_vector[:, 2] = np.radians(position_vector[:, 2])

    # Find the indices corresponding to the time interval
    unique_indices = np.unique(times, return_index=True)[1]
    times = times[unique_indices]
    position_vector = position_vector[unique_indices]

    interval_indices = np.where((new_times >= start_time) & (new_times <= end_time))[0]


    # Interpolate position functions for the specified time interval
    pos_interp_funcs = [CubicSpline(times, position_vector[:,i], extrapolate=True) for i in range(position_vector.shape[1])]

    pos_interp = np.column_stack([func(new_times[interval_indices]) for func in pos_interp_funcs])

    # Get cartesian versions for the specified time interval
    pos_interp_cart = np.apply_along_axis(s_c_vec_conversion, 1, np.copy(pos_interp))

    
    new_times = new_times[interval_indices]
    if coord_system == "spherical":
        return pos_interp
    elif coord_system == "cartesian":
        return pos_interp_cart
    else:
        raise ValueError("Invalid coordinate system specified. Must be 'spherical' or 'cartesian'.")

极化损耗计算函数 #

def clip_norm_dots(vec1, vec2):
     return np.clip(np.einsum('ij,ij->i',vec1, vec2),0,0.9999)
def unit_vector(vector):
    return vector / np.linalg.norm(vector)
def get_polarization_loss(receiver, source, separation):

    receiver_hat = unit_vector(receiver)
    source_hat = unit_vector(source)
    separation_hat = unit_vector(separation)
    dot_rs = clip_norm_dots(receiver_hat,source_hat)
    dot_rsep = clip_norm_dots(receiver_hat,separation_hat)
    dot_ssep = clip_norm_dots(source_hat,separation_hat)
    denominator = np.sqrt(1-dot_rsep**2)*np.sqrt(1-dot_ssep**2)
    denominator = np.maximum(denominator,1e-6)
    ploss = (dot_rs - (dot_rsep*dot_ssep))/denominator
    return ploss

最小复现示例参数 #

receiver = np.array([1,0,0])
source = np.array([1,1,0])

辅助函数与调用逻辑 #

坐标转换与时间插值

def s_c_vec_conversion(spherical_vec):
    x = spherical_vec[0] * np.sin(spherical_vec[1]) * np.cos(spherical_vec[2])
    y = spherical_vec[0] * np.sin(spherical_vec[1]) * np.sin(spherical_vec[2])
    z = spherical_vec[0] * np.cos(spherical_vec[1])
    return np.array([x, y, z])

def c_s_vec_conversion(cartesian_vec):
    r = np.sqrt(cartesian_vec[0]**2+cartesian_vec[1]**2+cartesian_vec[2]**2)
    theta = np.arccos(cartesian_vec[2]/r)
    phi = np.arctan(cartesian_vec[1]/cartesian_vec[0])
    return np.array([r,theta,phi])

def interp_time(times, sample_rate, start_time=0, end_time=None):
    if end_time is None:
        end_time = times[-1]
    new_times = np.arange(times[0], times[-1],1/sample_rate)
    interval_indices = np.where((new_times >= start_time) & (new_times <= end_time))[0]
    return new_times[interval_indices]

轨迹位置处理

# Translate a point in spherical coordinates relative to a reference point
def translate_point_spherical(r1, theta1, phi1, r2, theta2, phi2):
    # Convert reference point and satellite point to Cartesian
    p1_cartesian = spherical_to_cartesian(r1, theta1, phi1)
    p2_cartesian = spherical_to_cartesian(r2, theta2, phi2)
    # Calculate the translated Cartesian coordinates
    translated_cartesian = p2_cartesian - p1_cartesian
    # Convert back to spherical coordinates
    return cartesian_to_spherical(*translated_cartesian)

def get_spherical_position(receiver_lat, receiver_lon, traj_arrays):
    pos_vec = np.zeros((len(traj_arrays["Time"]),3))
    times = traj_arrays["Time"]
# Process trajectory data
    for i, time in enumerate(traj_arrays["Time"]):
        latitude = traj_arrays["Latgd"][i]
        longitude = traj_arrays["Long"][i]
        altitude = traj_arrays["Altkm"][i] * 1000  # Convert altitude to meters
        time = traj_arrays["Time"][i]
        
        # Poker Flat as origin
        r, theta, phi = translate_point_spherical(
            R_EARTH, np.pi/2 - np.radians(receiver_lat), np.radians(receiver_lon),
            R_EARTH + altitude, np.pi/2 - np.radians(latitude), np.radians(longitude)
        )
        pos_vec[i][0] = r
        pos_vec[i][1] = (np.degrees(np.arccos(altitude/r)))
        phi = np.degrees(phi)
        if phi < 0:
            phi += 180
        
        pos_vec[i][2] = (phi)
    return pos_vec

接收功率计算

def calc_received_power(rocket_pos, gains_rx, gains_tx, ploss):
     result_watts = (txPwr * np.multiply(np.multiply(gains_tx,gains_rx),ploss**2))/path_loss(rocket_pos)
     result_watts[result_watts<=0]=1e-100
     result_watts = result_watts.astype(np.float64)
     result_dBm = 10*np.log10(result_watts)+30
     return result_dBm

函数调用示例 #

traj_arrays = ny.read_traj_data("Traj_Right.txt")
times = ny.get_times(traj_arrays)
positions = np.zeros((len(times),len(Receiver),3))
for i in range(len(Receiver)):
    positions[:,i] = ny.get_spherical_position(coords[i][0], coords[i][1], traj_arrays)
losses = np.zeros((len(times_interp),len(Receiver),2))
for i in range(len(Receiver)):
    losses[:,i,0] = ut.get_polarization_loss(receivers_ew,transmitters_aligned,positions_interp_cartesian[:,i])
    losses[:,i,1] = ut.get_polarization_loss(receivers_ns,transmitters_aligned,positions_interp_cartesian[:,i])

优化建议 #

• 切换插值坐标系：当前对球面坐标分量分别插值再转笛卡尔坐标，角度分量的非线性插值易引入偏差。建议先将原始轨迹转为笛卡尔坐标，对x/y/z分量做插值，再按需转回球面坐标，能大幅减少位置突变。
• 更换插值方法：用scipy.interpolate.make_interp_spline生成B样条插值替代三次样条，B样条的平滑性更好，能降低插值点的波动。
• 平滑插值结果：对插值后的位置向量或ploss序列做滑动窗口平滑，比如使用scipy.signal.savgol_filter，参数可根据轨迹变化幅度调整，避免过度平滑丢失真实趋势。
• 基准验证：用极高采样率（如1000）的插值结果作为基准，对比不同采样率下的ploss平方偏差，筛选最优插值方案。

内容的提问来源于stack exchange，提问作者Sean Wallace