CN113204909B - 基于地基观测光度信号的卫星几何特征与姿态估计方法 - Google Patents

基于地基观测光度信号的卫星几何特征与姿态估计方法 Download PDF

Info

Publication number
CN113204909B
CN113204909B CN202110624730.9A CN202110624730A CN113204909B CN 113204909 B CN113204909 B CN 113204909B CN 202110624730 A CN202110624730 A CN 202110624730A CN 113204909 B CN113204909 B CN 113204909B
Authority
CN
China
Prior art keywords
satellite
attitude
observation
model
photometric
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110624730.9A
Other languages
English (en)
Other versions
CN113204909A (zh
Inventor
智喜洋
鲍广震
巩晋南
张寅�
胡建明
施天俊
张鹏飞
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin Institute of Technology
Original Assignee
Harbin Institute of Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin Institute of Technology filed Critical Harbin Institute of Technology
Priority to CN202110624730.9A priority Critical patent/CN113204909B/zh
Publication of CN113204909A publication Critical patent/CN113204909A/zh
Application granted granted Critical
Publication of CN113204909B publication Critical patent/CN113204909B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Navigation (AREA)

Abstract

本发明公开了一种基于地基观测光度信号的卫星几何特征与姿态估计方法，包括：S1，建立地基观测条件下的卫星光度信号观测模型；S2，建立卫星的“几何‑姿态‑光度”数据库：S3，辨识卫星的几何模型和尺度；S4，建立被观测卫星的姿态运动学方程；S5，设置无损卡尔曼滤波器的初始参数；S6，将姿态运动学方程和卫星光度信号观测模型分别作为无损卡尔曼滤波算法的时间更新模型和观测更新模型，采用无损卡尔曼滤波算法对下一观测时刻卫星姿态参数进行更新估计；S7，将步骤S6估计的卫星姿态参数作为新的观测时刻卫星状态参数重复步骤S6，当卫星姿态参数估计值误差小于设定阈值或卫星超出观测范围时停止迭代，从而获得高精度的卫星姿态参数估计值。

Description

基于地基观测光度信号的卫星几何特征与姿态估计方法
技术领域
本发明涉及空间光学信息处理技术领域，特别是涉及一种基于地基观测光度信号的卫星几何特征与姿态估计方法。
背景技术
近年来，日益紧张的空间资源和国际关系对各国空间态势感知能力提出更高、更新的要求。空间目标关键特征的准确辨识是空间态势感知的主要任务之一，对目标类型识别、状态判断、威胁评估和战略决策等任务具有重要意义。地基光度信号只包含能量信息，是空间目标多种关键特征耦合作用的结果，可综合体现空间目标几何形状、尺度与姿态的时序变化，已经成为中、高轨卫星几何特征辨识和姿态估计的重要手段。
由于空间目标表面材料在各方向上反射率不同，地基观测系统属于典型非线性系统，基于地基光度信号的目标特征参数估计问题可转化为对非线性系统未知状态参数的估计问题。目前，针对非线性系统参数估计问题的分析方法较多，其中以无损卡尔曼滤波(UKF)为代表的非线性卡尔曼滤波算法，具有可定量、可实时估计的优点，因此具有广泛的应用前景。然而，目前多数针对光度曲线的研究仍有一些局限，一方面，这些研究多是建立在目标几何模型已知的基础上的，但实际上，目标的准确几何特征是难以获取的，导致接下来的姿态估计工作难以开展；另一方面，有约束的四元数通常被用来描述空间目标的姿态，而UKF的估计对象是无约束的参数，因此，UKF无法直接应用于空间目标姿态参数的估计。
在几何模型数据库已知的情况下，可以通过轨迹搜索计算匹配度，该匹配度表征衡量目标属于某几何模型的可能性。同时，罗德里格斯参数作为表征姿态的另一种方式，具有无约束、小角度内稳定无奇点的优点，适用于表征姿态的采样距离。综上分析，考虑到地基观测系统优势以及无损卡尔曼滤波算法(UKF)能够对非线性动态系统状态参数进行准确估计，提供一种结合地基光学观测系统以及无损卡尔曼滤波算法的卫星几何模型辨识和姿态估计方法具有重要意义。
发明内容
本发明的目的是提供一种基于地基观测光度信号的卫星几何特征与姿态估计方法，以地基观测的卫星时序光度信号作为输入，利用卫星几何数据库、卫星姿态运动学方程对卫星的几何特征进行辨识，在此基础上进行空间目标姿态的状态更新，并利用卫星地基观测模型对卫星地基观测光度信号进行观测更新，采用无损卡尔曼滤波算法对卫星地基观测系统进行非线性滤波，进而实现卫星姿态的高精度估计。
为实现上述目的，本发明提供了如下方案：
一种基于地基观测光度信号的卫星几何特征与姿态估计方法，该方法包括以下步骤：
S1，从卫星光度信号形成机理出发，建立地基观测条件下的卫星光度信号观测模型，并利用该卫星光度信号观测模型对卫星光度信号进行观测更新；
S2，建立卫星的“几何-姿态-光度”数据库：
S3，基于观测到的卫星光度信号计算卫星自转周期，并计算卫星光度信号的采样序列与实际观测序列的匹配度，根据匹配度辨识卫星的几何模型和尺度；
S4，在确定卫星的几何模型和尺度的基础上，建立被观测卫星的姿态运动学方程；
S5，根据被观测卫星的姿态运动学方程与地基观测系统噪声先验信息设置无损卡尔曼滤波器的初始参数；
S6，获取当前观测时刻的卫星光度信号，将被观测卫星的姿态运动学方程和卫星光度信号观测模型分别作为无损卡尔曼滤波算法的时间更新模型和观测更新模型，采用无损卡尔曼滤波算法对下一观测时刻卫星姿态参数进行更新估计；
S7，将步骤S6估计的卫星姿态参数作为新的观测时刻卫星状态参数重复步骤S6，当卫星姿态参数估计值误差小于设定阈值或卫星超出观测范围时停止迭代，获得卫星姿态参数估计值。
进一步的，所述步骤S1中，从卫星光度信号形成机理出发，建立地基观测条件下的卫星光度信号观测模型，并利用该卫星光度信号观测模型对卫星光度信号进行观测更新，具体包括以下步骤：
S101，构建被观测卫星的几何模型，将构建的卫星几何模型划分为有限个微小面元，并对利用三个基向量
Figure BDA0003101738790000031
对空间目标划分的每块面元进行描述，其中，单位向量
Figure BDA0003101738790000032
指向面元外法线方向，单位向量
Figure BDA0003101738790000033
和
Figure BDA0003101738790000034
在面元平面内与单位向量
Figure BDA0003101738790000035
成右手坐标系；
S102，基于BRDF构建卫星表面反射模型，对卫星表面对太阳辐射的反射进行准确计算，卫星表面对太阳辐射的反射分为镜面反射ρspec和漫反射ρdiff两部分，卫星表面面元总的反射系数ρtotal为：
ρtotal(i)＝ρspec(i)+ρdiff(i) (1)
其中，卫星表面面元的漫反射系数ρdiff表达式为：
Figure BDA0003101738790000036
卫星表面面元的镜面反射系数ρspec表达式为：
Figure BDA0003101738790000037
式中，菲涅尔反射率
Figure BDA0003101738790000038
Figure BDA0003101738790000039
为ECI坐标系下卫星面元中心指向太阳方向的单位向量，
Figure BDA00031017387900000310
为ECI坐标系下卫星面元中心指向观测方向的单位向量，
Figure BDA00031017387900000311
为向量
Figure BDA00031017387900000312
与向量
Figure BDA00031017387900000313
的平分方向单位向量，常数nuv描述目标材料表面粗糙度其值越大表示材料表面越光滑越接近于完全镜面反射，Rspec为面元材料镜面反射系数，Rdiff为面元材料漫反射系数；
S103，对卫星表面面元是否被其他面元遮挡进行判别，并计算未被遮挡的卫星表面面元在地基观测系统产生的照度，地基观测站观测到卫星的星等值表达式为：
Figure BDA0003101738790000041
式中，-26.7是太阳的星等值，Esun,vis是可见光波段内太阳在卫星表面产生的辐射照度，Epupil(i)为卫星第i块面元在地基观测站光学系统入瞳处产生的辐射照度；
由目标姿态求解目标面元的基向量
Figure BDA0003101738790000042
因此定义由目标姿态参数计算目标星等的模型为观测模型，即
Figure BDA0003101738790000043
其中，ρtotal,i与太阳向量
Figure BDA0003101738790000044
观测向量
Figure BDA0003101738790000045
以及目标姿态q有关，q的定义见公式(7)。
进一步的，所述步骤S2，建立卫星的“几何-姿态-光度”数据库中，所述数据库中包括卫星几何模型以及每种卫星几何模型在各姿态下的光度，其中卫星几何模型为有限元模型，对应姿态以固定步长遍历所有参数，对应的光度值采用式(4-1)计算，这一步骤具体包括：
S201，根据斐波那契网格，确定数据库中欧拉旋转轴的均匀离散采样点，斐波那契网格的离散序列为
Figure BDA0003101738790000046
式中，N是采样点总数，黄金分割比
Figure BDA0003101738790000047
S202，在欧拉旋转轴确定的基础上确定欧拉旋转角的采样间距，采样步长为
Figure BDA0003101738790000048
其中，M是一个旋转轴上的采样点数目，在实际应用中取M＝8可以兼顾运算速度和匹配精度；
S203，根据四元数的定义，将欧拉旋转角和欧拉旋转轴转换为四元数，建立卫星的姿态数据库；
S204，将姿态数据库中的姿态作为输入参数，输入到步骤S1中建立的卫星光度信号观测模型，利用式(1)至(4-1)求解空间目标的光度信号，通过组合不同卫星几何模型、不同姿态角及不同条件下求解获得的光度信号建立几何-姿态-光度数据库。
进一步的，所述步骤S203中，根据四元数的定义，将欧拉旋转角和欧拉旋转轴转换为四元数，建立卫星的姿态数据库，具体包括：
四元数的定义为：
q＝[q0μ]T (7)
欧拉角参数到四元数的转换公式为：
Figure BDA0003101738790000051
式中，欧拉旋转轴
Figure BDA0003101738790000052
v是欧拉旋转角，q0、μ是组成四元数的中间变量，其中q0是标量，μ是向量。
进一步的，所述步骤S3，基于观测到的卫星光度信号计算卫星自转周期，并计算卫星光度信号的采样序列与实际观测序列的匹配度，根据匹配度辨识卫星的几何模型和尺度，具体包括：
S301，记卫星光度信号的观测序列为
Figure BDA0003101738790000053
“几何-姿态-光度”数据库中的第m个旋转轴上不同旋转角度的姿态四元数集合为Qm，对应的光度信号观测为Mm，基于数字傅里叶变换分析观测序列的周期性，计算其主频并记为f，则卫星自转的角速度为：
w＝2πf (9)
S302，将Qm中N个等间距的光度值在实际观测的曲线上滑动，在对应位置取商并计算标准差，该标准差作为曲线在当前位置的幅值，匹配度函数定义为:
Figure BDA0003101738790000061
其中，pn(t)是中间变量，定义如式(10)所示，
Figure BDA0003101738790000062
表示pn(t)的均值；
Figure BDA0003101738790000063
表示t时刻实际观测的数值，mn表示第i个初始姿态、第m个旋转轴、第n个旋转角的光度值，该光度值由步骤S2中的“几何-姿态-光度”数据库获取；
S303，根据匹配度函数计算匹配度，与实际观测值匹配度最高的几何模型输出为几何辨识结果，几何模型的面积为：
Figure BDA0003101738790000064
其中，A0是“几何-姿态-光度”数据库中几何模型的总面积。
进一步的，所述步骤S4，在确定卫星的几何模型和尺度的基础上，建立被观测卫星的姿态运动学方程，具体包括：
卫星的姿态由四元数进行描述，对于空间目标的连续姿态四元数的运动学模型为：
Figure BDA0003101738790000065
Figure BDA0003101738790000066
其中，q是姿态的四元数，ω是卫星自转的角速度，[p×]表示一种运算规则，定义为：
Figure BDA0003101738790000067
进一步的，所述步骤S5，根据被观测卫星的姿态运动学方程与地基观测系统噪声先验信息设置无损卡尔曼滤波器的初始参数，具体包括：
根据建立的被观测卫星的姿态运动学方程和地基观测系统的噪声先验信息分别设置卫星系统状态参量噪声矩阵Q和观测噪声矩阵R，其中Q表征系统对姿态运动学方程的信任程度，信任度越高，Q的取值越低，收敛速度越快，但是鲁棒性越低，R表征系统对地基观测模型的信任程度，信任度越高，实际观测对估计的状态参数修正力度越大，收敛速度越快，同时鲁棒性越低；
设定卫星全局姿态参数初始估计值x＝[q0,q1,q2,q3,q4,ωx,ωy,ωz]，卫星局部姿态参数初始估计值δx＝[p1,p2,p3,ωx,ωy,ωz]；q0,q1,q2,q3,q4分别为四元数的四个数值，ωx,ωy,ωz是卫星旋转角速度的三个分量；
其中Rod＝[p1,p2,p3]是罗德里格斯参数，与四元数之间的转换关系为：
Figure BDA0003101738790000071
在此基础上利用单位矩阵初始化卫星初始状态估计误差协方差矩阵P。
进一步的，所述步骤S6，获取当前观测时刻的卫星光度信号，将被观测卫星的姿态运动学方程和卫星光度信号观测模型分别作为无损卡尔曼滤波算法的时间更新模型和观测更新模型，采用无损卡尔曼滤波算法对下一观测时刻卫星姿态参数进行更新估计；
S601，对卫星在当前观测时刻的Na维局部姿态估计参数向量
Figure BDA0003101738790000072
进行采样，获得局部姿态参数的σ采样点和对应的权重分布，卫星状态参数向量的2Na+1维σ采样点序列为：
Figure BDA0003101738790000073
σ采样点序列对应的均值计算权重为：
Figure BDA0003101738790000081
σ采样点序列对应的协方差计算权重为：
Figure BDA0003101738790000082
式中，常数γ控制σ点分布的尺寸，且满足0≤γ≤1一般取值为0.5，常数参数λ＝3γ2-Na，
Figure BDA0003101738790000083
表示对矩阵·的开根后的第i列，λ为中间变量；
S602，将卫星在当前时刻的用罗德里格斯参数表征的σ采样点的局部姿态参数
Figure BDA0003101738790000084
转换为四元数表征的全局姿态参数
Figure BDA0003101738790000085
S603，将四元数表征的全局姿态参数
Figure BDA0003101738790000086
代入运动学模型，计算卫星在下一观测时刻的姿态参数预测值
Figure BDA0003101738790000087
预测值均值
Figure BDA0003101738790000088
和预测方差矩阵
Figure BDA0003101738790000089
Figure BDA00031017387900000810
Figure BDA00031017387900000811
Figure BDA00031017387900000812
S604，利用卫星下一观测时刻姿态参数预测值
Figure BDA00031017387900000813
和卫星地基光度信号观测模型，计算卫星地基观测条件下光度信号预测值yk+1和预测值均值
Figure BDA00031017387900000814
Figure BDA00031017387900000815
Figure BDA00031017387900000816
式中，式中，h是观测函数，由式(4-2)计算得到；
S605，得到光度信号预测值后，根据σ采样点的协方差权重因子Wi cov计算预测值的自协方差矩阵
Figure BDA00031017387900000817
和预测值与卫星状态量的互协方差矩阵
Figure BDA00031017387900000818
Figure BDA0003101738790000091
Figure BDA0003101738790000092
S606，通过计算卡尔曼增益κk+1对卫星姿态状态参量
Figure BDA0003101738790000093
和状态参量估计协方差矩阵Pk+1进行更新计算：
Figure BDA0003101738790000094
Figure BDA0003101738790000095
Figure BDA0003101738790000096
根据本发明提供的具体实施例，本发明提供的基于地基观测光度信号的卫星几何特征与姿态估计方法，与现有技术相比具有如下技术效果：
(1)针对现有卫星姿态参数估计方法中的几何模型先验不确定和无损卡尔曼滤波不能直接应用于有约束的四元数的问题，加之考虑了地基光学系统的不受体积和质量制约、实现难度小且成本低等优势和难以对细节成像的局限性，提出了基于地基光度信号的卫星几何特征辨识和姿态估计方法，本发明利用少量信息完成卫星的几何特征和尺度辨识工作，在此基础上，将基于地基观测光度信号的卫星姿态参数估计转化为非线性系统未知状态参数的估计问题，进而采用无损卡尔曼滤波算法实现卫星姿态参数的高精度估计，对于空间目标监视体系建设具有一定的科学指导意义和应用价值；
(2)本发明从卫星光度信号形成的物理机理出发，通过建立卫星的几何模型，并基于双向反射分布函数(BRDF)建立卫星表面反射模型，对卫星表面反射的辐射能量进行准确计算，进而实现对地基观测条件下卫星光度信号的形成过程进行精准建模表征，该方法能够更加准确地描述地基观测条件下卫星光度信号形成的物理过程，从而提高地基观测条件下卫星光度信号的仿真准确度；
(3)本发明针对不同集合模型建立数据库，同时基于该数据库利用路径搜寻的方法计算该几何模型与目标卫星的匹配度，同时利用匹配函数可以求解卫星的尺度，该方法不依赖逐点匹配，抑制了几何突变、大气传输、成像系统等因素产生的观测噪声，相较于逐点匹配的方法具有更好的准确性和鲁棒性。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例基于地基观测光度信号的卫星几何特征与姿态估计方法的流程图；
图2为本发明实施例卫星表面面元反射几何示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明的目的是提供一种基于地基观测光度信号的卫星几何特征与姿态估计方法，以地基观测的卫星时序光度信号作为输入，利用卫星几何数据库、卫星姿态运动学方程对卫星的几何特征进行辨识，在此基础上进行空间目标姿态的状态更新，并利用卫星地基观测模型对卫星地基观测光度信号进行观测更新，采用无损卡尔曼滤波算法对卫星地基观测系统进行非线性滤波，进而实现卫星姿态的高精度估计。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
如图1所示，本发明实施例提供的基于地基观测光度信号的卫星几何特征与姿态估计方法，包括：
S1，从卫星光度信号形成机理出发，建立地基观测条件下的卫星光度信号观测模型，并利用该卫星光度信号观测模型对卫星光度信号进行观测更新；
S2，建立卫星的“几何-姿态-光度”数据库：
S3，基于观测到的卫星光度信号计算卫星自转周期，并计算卫星光度信号的采样序列与实际观测序列的匹配度，根据匹配度辨识卫星的几何模型和尺度；
S4，在确定卫星的几何模型和尺度的基础上，建立被观测卫星的姿态运动学方程；
S5，根据被观测卫星的姿态运动学方程与地基观测系统噪声先验信息设置无损卡尔曼滤波器的初始参数；
S6，获取当前观测时刻的卫星光度信号，将被观测卫星的姿态运动学方程和卫星光度信号观测模型分别作为无损卡尔曼滤波算法的时间更新模型和观测更新模型，采用无损卡尔曼滤波算法对下一观测时刻卫星姿态参数进行更新估计；
S7，将步骤S6估计的卫星姿态参数作为新的观测时刻卫星状态参数重复步骤S6，当卫星姿态参数估计值误差小于设定阈值或卫星超出观测范围时停止迭代，获得卫星姿态参数估计值。
其中，所述步骤S1中，从卫星光度信号形成机理出发，建立地基观测条件下的卫星光度信号观测模型，并利用该卫星光度信号观测模型对卫星光度信号进行观测更新，具体包括以下步骤：
S101，构建被观测卫星的几何模型，将构建的卫星几何模型划分为有限个微小面元，并对利用三个基向量
Figure BDA0003101738790000111
对空间目标划分的每块面元进行描述，如图2所示，其中，单位向量
Figure BDA0003101738790000112
指向面元外法线方向，单位向量
Figure BDA0003101738790000113
和
Figure BDA0003101738790000114
在面元平面内与单位向量
Figure BDA0003101738790000115
成右手坐标系；
S102，基于BRDF构建卫星表面反射模型，对卫星表面对太阳辐射的反射进行准确计算，卫星表面对太阳辐射的反射分为镜面反射ρspec和漫反射ρdiff两部分，卫星表面面元总的反射系数ρtotal为：
ρtotal(i)＝ρspec(i)+ρdiff(i) (1)
其中，卫星表面面元的漫反射系数ρdiff表达式为：
Figure BDA0003101738790000116
卫星表面面元的镜面反射系数ρspec表达式为：
Figure BDA0003101738790000121
式中，菲涅尔反射率
Figure BDA0003101738790000122
Figure BDA0003101738790000123
为ECI坐标系下卫星面元中心指向太阳方向的单位向量，
Figure BDA0003101738790000124
为ECI坐标系下卫星面元中心指向观测方向的单位向量，
Figure BDA0003101738790000125
为向量
Figure BDA0003101738790000126
与向量
Figure BDA0003101738790000127
的平分方向单位向量，常数nuv描述目标材料表面粗糙度其值越大表示材料表面越光滑越接近于完全镜面反射，Rspec为面元材料镜面反射系数，Rdiff为面元材料漫反射系数；
S103，对卫星表面面元是否被其他面元遮挡进行判别，并计算未被遮挡的卫星表面面元在地基观测系统产生的照度，地基观测站观测到卫星的星等值表达式为：
Figure BDA0003101738790000128
式中，-26.7是太阳的星等值，Esun,vis是可见光波段内太阳在卫星表面产生的辐射照度，Epupil(i)为卫星第i块面元在地基观测站光学系统入瞳处产生的辐射照度；
由目标姿态求解目标面元的基向量
Figure BDA0003101738790000129
因此定义由目标姿态参数计算目标星等的模型为观测模型，即
Figure BDA00031017387900001210
其中，ρtotal,i与太阳向量
Figure BDA00031017387900001211
观测向量
Figure BDA00031017387900001212
以及目标姿态q有关，q的定义见公式(7)。
所述步骤S2，所述步骤S2，建立卫星的“几何-姿态-光度”数据库中，所述数据库中包括卫星几何模型以及每种卫星几何模型在各姿态下的光度，其中卫星几何模型为有限元模型，对应姿态以固定步长遍历所有参数，对应的光度值采用式(4-1)计算，这一步骤具体包括：
S201，根据斐波那契网格，确定数据库中欧拉旋转轴的均匀离散采样点，斐波那契网格的离散序列为
Figure BDA0003101738790000131
式中，N是采样点总数，黄金分割比
Figure BDA0003101738790000132
S202，在欧拉旋转轴确定的基础上确定欧拉旋转角的采样间距，采样步长为
Figure BDA0003101738790000133
其中，M是一个旋转轴上的采样点数目，在实际应用中取M＝8可以兼顾运算速度和匹配精度；
S203，根据四元数的定义，将欧拉旋转角和欧拉旋转轴转换为四元数，建立卫星的姿态数据库；具体包括：
四元数的定义为：
q＝[q0μ]T (7)
欧拉角参数到四元数的转换公式为：
Figure BDA0003101738790000134
式中，欧拉旋转轴
Figure BDA0003101738790000135
v是欧拉旋转角，q0、μ是组成四元数的中间变量，其中q0是标量，μ是向量。
S204，将姿态数据库中的姿态作为输入参数，输入到步骤S1中建立的卫星光度信号观测模型，利用式(1)至(4-1)求解空间目标的光度信号，通过组合不同卫星几何模型、不同姿态角及不同条件下求解获得的光度信号建立几何-姿态-光度数据库。
所述步骤S3，基于观测到的卫星光度信号计算卫星自转周期，并计算卫星光度信号的采样序列与实际观测序列的匹配度，根据匹配度辨识卫星的几何模型和尺度，具体包括：
S301，记卫星光度信号的观测序列为
Figure BDA0003101738790000141
“几何-姿态-光度”数据库中的第m个旋转轴上不同旋转角度的姿态四元数集合为Qm，对应的光度信号观测为Mm，基于数字傅里叶变换分析观测序列的周期性，计算其主频并记为f，则卫星自转的角速度为：
w＝2πf (9)
由于空间目标的自转是周期性的，周期匹配法不需要考虑目标的具体起始姿态，仅须考虑旋转轴即可，因此周期匹配是针对每个欧拉轴进行的，针对欧拉轴的采样点一个周期内的旋转轨迹进行等间隔采样：
Figure BDA0003101738790000142
其中，N是采样点数，θj表示第j个旋转角度；
S302，针对第i个初始姿态，第j个旋转轴的观测与实际观测匹配度由匹配度函数曲线进行描述，将N个光度值在实际观测的曲线上滑动，在对应位置取商并计算标准差，该标准差作为曲线在当前位置的幅值，匹配度函数定义为:
Figure BDA0003101738790000143
其中，pn(t)是中间变量，定义如式(10)所示，
Figure BDA0003101738790000144
表示pn(t)的均值；
Figure BDA0003101738790000145
表示t时刻实际观测的数值，mn表示第i个初始姿态、第m个旋转轴、第n个旋转角的光度值，该光度值由步骤S2中的“几何-姿态-光度”数据库获取；
S303，根据匹配度函数计算匹配度，与实际观测值匹配度最高的几何模型输出为几何辨识结果，几何模型的面积为：
Figure BDA0003101738790000146
其中，A0是“几何-姿态-光度”数据库中几何模型的总面积，t＝argmin(f)。t是令式(10)中f取值最小时的值，表示观测到的光度曲线对应的几何模型表面积是“几何-姿态-光度”数据库中几何模型的总面积的
Figure BDA0003101738790000151
倍
所述步骤S4，在确定卫星的几何模型和尺度的基础上，建立被观测卫星的姿态运动学方程，具体包括：
卫星的姿态由四元数进行描述，对于空间目标的连续姿态四元数的运动学模型为：
Figure BDA0003101738790000152
Figure BDA0003101738790000153
其中，q是姿态的四元数，ω是卫星自转的角速度，[p×]表示一种运算规则，定义为：
Figure BDA0003101738790000154
所述步骤S5，根据被观测卫星的姿态运动学方程与地基观测系统噪声先验信息设置无损卡尔曼滤波器的初始参数，具体包括：
根据建立的被观测卫星的姿态运动学方程和地基观测系统的噪声先验信息分别设置卫星系统状态参量噪声矩阵Q和观测噪声矩阵R，其中Q表征系统对姿态运动学方程的信任程度，信任度越高，Q的取值越低，收敛速度越快，但是鲁棒性越低，R表征系统对地基观测模型的信任程度，信任度越高，实际观测对估计的状态参数修正力度越大，收敛速度越快，同时鲁棒性越低；
设定卫星全局姿态参数初始估计值x＝[q0,q1,q2,q3,q4,ωx,ωy,ωz]，卫星局部姿态参数初始估计值δx＝[p1,p2,p3,ωx,ωy,ωz]；q0,q1,q2,q3,q4分别为四元数的四个数值，ωx,ωy,ωz是卫星旋转角速度的三个分量；
其中Rod＝[p1,p2,p3]是罗德里格斯参数，与四元数之间的转换关系为：
Figure BDA0003101738790000161
在此基础上利用单位矩阵初始化卫星初始状态估计误差协方差矩阵P。
所述步骤S6，获取当前观测时刻的卫星光度信号，将被观测卫星的姿态运动学方程和卫星光度信号观测模型分别作为无损卡尔曼滤波算法的时间更新模型和观测更新模型，采用无损卡尔曼滤波算法对下一观测时刻卫星姿态参数进行更新估计；
S601，对卫星在当前观测时刻的Na维局部姿态估计参数向量
Figure BDA0003101738790000162
进行采样，获得局部姿态参数的σ采样点和对应的权重分布，卫星状态参数向量的2Na+1维σ采样点序列为：
Figure BDA0003101738790000163
σ采样点序列对应的均值计算权重为：
Figure BDA0003101738790000164
σ采样点序列对应的协方差计算权重为：
Figure BDA0003101738790000165
式中，常数γ控制σ点分布的尺寸，且满足0≤γ≤1一般取值为0.5，常数参数λ＝3γ2-Na，
Figure BDA0003101738790000166
表示对矩阵·的开根后的第i列，λ为中间变量；
S602，将卫星在当前时刻的用罗德里格斯参数表征的σ采样点的局部姿态参数
Figure BDA0003101738790000171
转换为四元数表征的全局姿态参数
Figure BDA0003101738790000172
S603，将四元数表征的全局姿态参数
Figure BDA0003101738790000173
代入运动学模型，计算卫星在下一观测时刻的姿态参数预测值
Figure BDA0003101738790000174
预测值均值
Figure BDA0003101738790000175
和预测方差矩阵
Figure BDA0003101738790000176
Figure BDA0003101738790000177
Figure BDA0003101738790000178
Figure BDA0003101738790000179
S604，利用卫星下一观测时刻姿态参数预测值
Figure BDA00031017387900001710
和卫星地基光度信号观测模型，计算卫星地基观测条件下光度信号预测值yk+1和预测值均值
Figure BDA00031017387900001711
Figure BDA00031017387900001712
Figure BDA00031017387900001713
式中，式中，h是观测函数，由式(4-2)计算得到；
S605，得到光度信号预测值后，根据σ采样点的协方差权重因子Wi cov计算预测值的自协方差矩阵
Figure BDA00031017387900001714
和预测值与卫星状态量的互协方差矩阵
Figure BDA00031017387900001715
Figure BDA00031017387900001716
Figure BDA00031017387900001717
S606，通过计算卡尔曼增益κk+1对卫星姿态状态参量
Figure BDA00031017387900001718
和状态参量估计协方差矩阵Pk+1进行更新计算：
Figure BDA00031017387900001719
Figure BDA00031017387900001720
Figure BDA00031017387900001721
本发明提供的基于地基观测光度信号的卫星几何特征与姿态估计方法，利用数据库的几何模型进行轨迹搜寻判断其几何模型，在确定几何模型的类型和尺度的基础上，利用卫星的运动学方程对卫星的姿态参数进行状态更新，从卫星光度信号形成机理出发，建立地基观测条件下的卫星光度信号观测模型，并利用观测模型对卫星光度信号进行观测更新，采用无损卡尔曼滤波算法对卫星姿态参数进行估计，该方法适用于基于地基观测光度信号的卫星几何特征估计和姿态参数的高精度估计。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (6)

1.一种基于地基观测光度信号的卫星几何特征与姿态估计方法，其特征在于，包括以下步骤：
S1，从卫星光度信号形成机理出发，建立地基观测条件下的卫星光度信号观测模型，并利用该卫星光度信号观测模型对卫星光度信号进行观测更新；
S2，建立卫星的“几何-姿态-光度”数据库：
S3，基于观测到的卫星光度信号计算卫星自转周期，并计算卫星光度信号的采样序列与实际观测序列的匹配度，根据匹配度辨识卫星的几何模型和尺度；
S4，在确定卫星的几何模型和尺度的基础上，建立被观测卫星的姿态运动学方程；具体包括：
卫星的姿态由四元数进行描述，对于空间目标的连续姿态四元数的运动学模型为：
Figure FDA0003610151670000011
Figure FDA0003610151670000012
其中，q是姿态的四元数，ω是卫星自转的角速度，[p×]表示一种运算规则，定义为：
Figure FDA0003610151670000013
S5，根据被观测卫星的姿态运动学方程与地基观测系统噪声先验信息设置无损卡尔曼滤波器的初始参数；具体包括：
根据建立的被观测卫星的姿态运动学方程和地基观测系统的噪声先验信息分别设置卫星系统状态参量噪声矩阵Q和观测噪声矩阵R，其中Q表征系统对姿态运动学方程的信任程度，信任度越高，Q的取值越低，收敛速度越快，但是鲁棒性越低，R表征系统对地基观测模型的信任程度，信任度越高，实际观测对估计的状态参数修正力度越大，收敛速度越快，同时鲁棒性越低；
设定卫星全局姿态参数初始估计值x＝[q0,q1,q2,q3,q4,ωx,ωy,ωz]，卫星局部姿态参数初始估计值δx＝[p1,p2,p3,ωx,ωy,ωz]；q0,q1,q2,q3,q4分别为四元数的四个数值，ωx,ωy,ωz是卫星旋转角速度的三个分量；
其中Rod＝[p1,p2,p3]是罗德里格斯参数，与四元数之间的转换关系为：
Figure FDA0003610151670000021
在此基础上利用单位矩阵初始化卫星初始状态估计误差协方差矩阵P；
S6，获取当前观测时刻的卫星光度信号，将被观测卫星的姿态运动学方程和卫星光度信号观测模型分别作为无损卡尔曼滤波算法的时间更新模型和观测更新模型，采用无损卡尔曼滤波算法对下一观测时刻卫星姿态参数进行更新估计；
S7，将步骤S6估计的卫星姿态参数作为新的观测时刻卫星状态参数重复步骤S6，当卫星姿态参数估计值误差小于设定阈值或卫星超出观测范围时停止迭代，获得卫星姿态参数估计值。
2.根据权利要求1所述的基于地基观测光度信号的卫星几何特征与姿态估计方法，其特征在于，所述步骤S1中，从卫星光度信号形成机理出发，建立地基观测条件下的卫星光度信号观测模型，并利用该卫星光度信号观测模型对卫星光度信号进行观测更新，具体包括以下步骤：
S101，构建被观测卫星的几何模型，将构建的卫星几何模型划分为有限个微小面元，并对利用三个基向量
Figure FDA0003610151670000022
对空间目标划分的每块面元进行描述，其中，单位向量
Figure FDA0003610151670000023
指向面元外法线方向，单位向量
Figure FDA0003610151670000024
和
Figure FDA0003610151670000025
在面元平面内与单位向量
Figure FDA0003610151670000026
成右手坐标系；
S102，基于BRDF构建卫星表面反射模型，对卫星表面对太阳辐射的反射进行准确计算，卫星表面对太阳辐射的反射分为镜面反射ρspec和漫反射ρdiff两部分，卫星表面面元总的反射系数ρtotal为：
ρtotal(i)＝ρspec(i)+ρdiff(i) (1)
其中，卫星表面面元的漫反射系数ρdiff表达式为：
Figure FDA0003610151670000031
卫星表面面元的镜面反射系数ρspec表达式为：
Figure FDA0003610151670000032
式中，菲涅尔反射率
Figure FDA0003610151670000033
Figure FDA0003610151670000034
为ECI坐标系下卫星面元中心指向太阳方向的单位向量，
Figure FDA0003610151670000035
为ECI坐标系下卫星面元中心指向观测方向的单位向量，
Figure FDA0003610151670000036
为向量
Figure FDA0003610151670000037
与向量
Figure FDA0003610151670000038
的平分方向单位向量，常数nuv描述目标材料表面粗糙度其值越大表示材料表面越光滑越接近于完全镜面反射，Rspec为面元材料镜面反射系数，Rdiff为面元材料漫反射系数；
S103，对卫星表面面元是否被其他面元遮挡进行判别，并计算未被遮挡的卫星表面面元在地基观测系统产生的照度，地基观测站观测到卫星的星等值表达式为：
Figure FDA0003610151670000039
式中，-26.7是太阳的星等值，Esun,vis是可见光波段内太阳在卫星表面产生的辐射照度，Epupil(i)为卫星第i块面元在地基观测站光学系统入瞳处产生的辐射照度；
由目标姿态求解目标面元的基向量
Figure FDA00036101516700000310
因此定义由目标姿态参数计算目标星等的模型为观测模型，即
Figure FDA00036101516700000311
其中，ρtotal,i与太阳向量
Figure FDA00036101516700000312
观测向量
Figure FDA00036101516700000313
以及目标姿态q有关。
3.根据权利要求2所述的基于地基观测光度信号的卫星几何特征与姿态估计方法，其特征在于，所述步骤S2，建立卫星的“几何-姿态-光度”数据库中，所述数据库中包括卫星几何模型以及每种卫星几何模型在各姿态下的光度，其中卫星几何模型为有限元模型，对应姿态以固定步长遍历所有参数，对应的光度值采用式(4-1)计算，这一步骤具体包括：
S201，根据斐波那契网格，确定数据库中欧拉旋转轴的均匀离散采样点，斐波那契网格的离散序列为
Figure FDA0003610151670000041
式中，N是采样点总数，黄金分割比
Figure FDA0003610151670000042
S202，在欧拉旋转轴确定的基础上确定欧拉旋转角的采样间距，采样步长为
Figure FDA0003610151670000043
其中，M是一个旋转轴上的采样点数目，取M＝8；
S203，根据四元数的定义，将欧拉旋转角和欧拉旋转轴转换为四元数，建立卫星的姿态数据库；
S204，将姿态数据库中的姿态作为输入参数，输入到步骤S1中建立的卫星光度信号观测模型，利用式(1)至(4-1)求解空间目标的光度信号，通过组合不同卫星几何模型、不同姿态角及不同条件下求解获得的光度信号建立几何-姿态-光度数据库。
4.根据权利要求3所述的基于地基观测光度信号的卫星几何特征与姿态估计方法，其特征在于，所述步骤S203中，根据四元数的定义，将欧拉旋转角和欧拉旋转轴转换为四元数，建立卫星的姿态数据库，具体包括：
四元数的定义为：
q＝[q0 μ]T (7)
欧拉角参数到四元数的转换公式为：
Figure FDA0003610151670000051
式中，欧拉旋转轴
Figure FDA0003610151670000052
v是欧拉旋转角，q0、μ是组成四元数的中间变量，其中q0是标量，μ是向量。
5.根据权利要求4所述的基于地基观测光度信号的卫星几何特征与姿态估计方法，其特征在于，所述步骤S3，基于观测到的卫星光度信号计算卫星自转周期，并计算卫星光度信号的采样序列与实际观测序列的匹配度，根据匹配度辨识卫星的几何模型和尺度，具体包括：
S301，记卫星光度信号的观测序列为
Figure FDA0003610151670000053
“几何-姿态-光度”数据库中的第m个旋转轴上不同旋转角度的姿态四元数集合为Qm，对应的光度信号观测为Mm，基于数字傅里叶变换分析观测序列的周期性，计算其主频并记为f，则卫星自转的角速度为：
w＝2πf (9)
S302，将Qm中N个等间距的光度值在实际观测的曲线上滑动，在对应位置取商并计算标准差，该标准差作为曲线在当前位置的幅值，匹配度函数定义为:
Figure FDA0003610151670000054
其中，pn(t)是中间变量，定义如式(10)所示，
Figure FDA0003610151670000055
表示pn(t)的均值；
Figure FDA0003610151670000056
表示t时刻实际观测的数值，mn表示第i个初始姿态、第m个旋转轴、第n个旋转角的光度值，该光度值由步骤S2中的“几何-姿态-光度”数据库获取；
S303，根据匹配度函数计算匹配度，与实际观测值匹配度最高的几何模型输出为几何辨识结果，几何模型的面积为：
Figure FDA0003610151670000061
其中，A0是“几何-姿态-光度”数据库中几何模型的总面积。
6.根据权利要求1所述的基于地基观测光度信号的卫星几何特征与姿态估计方法，其特征在于，所述步骤S6，获取当前观测时刻的卫星光度信号，将被观测卫星的姿态运动学方程和卫星光度信号观测模型分别作为无损卡尔曼滤波算法的时间更新模型和观测更新模型，采用无损卡尔曼滤波算法对下一观测时刻卫星姿态参数进行更新估计；
S601，对卫星在当前观测时刻的Na维局部姿态估计参数向量
Figure FDA0003610151670000062
进行采样，获得局部姿态参数的σ采样点和对应的权重分布，卫星状态参数向量的2Na+1维σ采样点序列为：
Figure FDA0003610151670000063
σ采样点序列对应的均值计算权重为：
Figure FDA0003610151670000064
σ采样点序列对应的协方差计算权重为：
Figure FDA0003610151670000065
式中，常数γ控制σ点分布的尺寸，且满足0≤γ≤1一般取值为0.5，常数参数λ＝3γ2-Na，
Figure FDA0003610151670000066
表示对矩阵·的开根后的第i列，λ为中间变量；
S602，将卫星在当前时刻的用罗德里格斯参数表征的σ采样点的局部姿态参数
Figure FDA0003610151670000067
转换为四元数表征的全局姿态参数
Figure FDA0003610151670000068
S603，将四元数表征的全局姿态参数
Figure FDA0003610151670000071
代入运动学模型，计算卫星在下一观测时刻的姿态参数预测值
Figure FDA0003610151670000072
预测值均值
Figure FDA0003610151670000073
和预测方差矩阵
Figure FDA0003610151670000074
Figure FDA0003610151670000075
Figure FDA0003610151670000076
Figure FDA0003610151670000077
S604，利用卫星下一观测时刻姿态参数预测值
Figure FDA0003610151670000078
和卫星地基光度信号观测模型，计算卫星地基观测条件下光度信号预测值yk+1和预测值均值
Figure FDA0003610151670000079
Figure FDA00036101516700000710
Figure FDA00036101516700000711
式中，h是观测函数，由式(4-2)计算得到；
S605，得到光度信号预测值后，根据σ采样点的协方差权重因子Wi cov计算预测值的自协方差矩阵
Figure FDA00036101516700000712
和预测值与卫星状态量的互协方差矩阵
Figure FDA00036101516700000713
Figure FDA00036101516700000714
Figure FDA00036101516700000715
S606，通过计算卡尔曼增益κk+1对卫星姿态状态参量
Figure FDA00036101516700000716
和状态参量估计协方差矩阵Pk+1进行更新计算：
Figure FDA00036101516700000717
Figure FDA00036101516700000718
Figure FDA00036101516700000719
CN202110624730.9A 2021-06-04 2021-06-04 基于地基观测光度信号的卫星几何特征与姿态估计方法 Active CN113204909B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110624730.9A CN113204909B (zh) 2021-06-04 2021-06-04 基于地基观测光度信号的卫星几何特征与姿态估计方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110624730.9A CN113204909B (zh) 2021-06-04 2021-06-04 基于地基观测光度信号的卫星几何特征与姿态估计方法

Publications (2)

Publication Number Publication Date
CN113204909A CN113204909A (zh) 2021-08-03
CN113204909B true CN113204909B (zh) 2022-07-19

Family

ID=77024429

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110624730.9A Active CN113204909B (zh) 2021-06-04 2021-06-04 基于地基观测光度信号的卫星几何特征与姿态估计方法

Country Status (1)

Country Link
CN (1) CN113204909B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113758373B (zh) * 2021-09-07 2023-01-10 重庆天箭惯性科技股份有限公司 一种提高弹载接收机定位、测速精度的方法、装置及设备
CN114111806B (zh) * 2022-01-21 2022-04-26 中国人民解放军32035部队 基于光度频谱特征的空间目标姿态稳定性估计方法及装置

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106595673A (zh) * 2016-12-12 2017-04-26 东南大学 面对地球静止轨道目标操作的空间多机器人自主导航方法
CN110109470A (zh) * 2019-04-09 2019-08-09 西安电子科技大学 基于无迹卡尔曼滤波的联合定姿方法、卫星姿态控制系统
CN110285815A (zh) * 2019-05-28 2019-09-27 山东航天电子技术研究所 一种可在轨全程应用的微纳卫星多源信息姿态确定方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
FR2975180B1 (fr) * 2011-05-09 2019-05-31 Airbus Defence And Space Sas Dispositif et procede de determination d'attitude d'un satellite et satellite embarquant un tel dispositif.

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106595673A (zh) * 2016-12-12 2017-04-26 东南大学 面对地球静止轨道目标操作的空间多机器人自主导航方法
CN110109470A (zh) * 2019-04-09 2019-08-09 西安电子科技大学 基于无迹卡尔曼滤波的联合定姿方法、卫星姿态控制系统
CN110285815A (zh) * 2019-05-28 2019-09-27 山东航天电子技术研究所 一种可在轨全程应用的微纳卫星多源信息姿态确定方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
魏小峰等.空间目标三维姿态估计方法.《武汉大学学报(信息科学版)》.2015,(第01期),全文. *

Also Published As

Publication number Publication date
CN113204909A (zh) 2021-08-03

Similar Documents

Publication Publication Date Title
CN111985093B (zh) 一种带噪声估计器的自适应无迹卡尔曼滤波状态估计方法
CN113204909B (zh) 基于地基观测光度信号的卫星几何特征与姿态估计方法
Subramani et al. Stochastic time-optimal path-planning in uncertain, strong, and dynamic flows
CN106707768B (zh) 一种基于模型校准的喷漆轨迹规划方法
CN109612472B (zh) 一种深空探测器自主导航系统构建方法及装置
CN111798491B (zh) 一种基于Elman神经网络的机动目标跟踪方法
CN111680870B (zh) 目标运动轨迹质量综合评估方法
US11294063B2 (en) System and method for fast wind flow measurement by LiDAR in a complex terrain
CN105678076A (zh) 点云测量数据质量评估优化的方法及装置
CN105738915A (zh) 三维雷达测量方法及装置
Yang et al. Bridge dynamic displacement monitoring using adaptive data fusion of GNSS and accelerometer measurements
JP2023517147A (ja) 複雑地形においてＬｉＤＡＲで風の流れの乱流を測定するためのシステムおよび方法
CN110231619B (zh) 基于恩克法的雷达交接时刻预报方法及装置
Proletarsky et al. Method for improving accuracy of INS using scalar parametric identification
CN114415157A (zh) 一种基于水声传感器网络的水下目标多模型跟踪方法
CN112346032A (zh) 基于一致性扩展卡尔曼滤波的单红外传感器目标定轨方法
CN108469604B (zh) Tws雷达空时联合数据模拟方法及系统
CN112926237B (zh) 一种基于光度信号的空间目标关键特征辨识方法
EP0580140B1 (en) Covert ranging method and system
Yun Sequential Monte Carlo filtering with Gaussian mixture models for highly nonlinear systems
Qin et al. Non-line-of-sight multipath classification method for BDS using convolutional sparse autoencoder with LSTM
CN117744742A (zh) 面向微小型无人机平台的深度学习模型压缩方法及系统
CN112346033A (zh) 一种针对量测数据带零偏的单红外传感器目标定轨方法
Tang et al. Underwater target tracking algorithm based on nonlinear RTS smoothing SVD-UKF with adaptive coefficients
CN114861122A (zh) 一种半分析半数值方式的轨道计算方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant