CN112965128B - 一种无奇异性顾及局域保障条件计算外部重力异常的方法 - Google Patents

一种无奇异性顾及局域保障条件计算外部重力异常的方法 Download PDF

Info

Publication number
CN112965128B
CN112965128B CN202110181838.5A CN202110181838A CN112965128B CN 112965128 B CN112965128 B CN 112965128B CN 202110181838 A CN202110181838 A CN 202110181838A CN 112965128 B CN112965128 B CN 112965128B
Authority
CN
China
Prior art keywords
gravity anomaly
external
gravity
calculating
point
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
CN202110181838.5A
Other languages
English (en)
Other versions
CN112965128A (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.)
92859 TROOPS PLA
Original Assignee
92859 TROOPS PLA
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 92859 TROOPS PLA filed Critical 92859 TROOPS PLA
Priority to CN202110181838.5A priority Critical patent/CN112965128B/zh
Publication of CN112965128A publication Critical patent/CN112965128A/zh
Application granted granted Critical
Publication of CN112965128B publication Critical patent/CN112965128B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V7/00Measuring gravitational fields or waves; Gravimetric prospecting or detecting
    • G01V7/02Details
    • G01V7/06Analysis or interpretation of gravimetric records
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10Complex mathematical operations
    • G06F17/11Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems

Abstract

本发明涉及一种无奇异性顾及局域保障条件计算外部重力异常的方法,针对计算外部重力异常的奇异性问题,通过采用移去‑恢复技术,从流动点重力异常中移去外部空间计算点在球面上的投影点处的重力异常;同时通过积分恒等变换,恢复外部空间计算点在球面上的投影点处的重力异常的全球影响,消除积分核函数奇异性的影响;顾及实测数据局域保障条件,引入全球重力场位模型,恢复参考外部重力异常;利用积分恒等式转换,补偿由全球积分过渡到局域积分引起的模型误差,得到计算点处高精度的外部重力异常。本发明能够精确计算外部重力异常,利用全球重力场位模型建立的模拟标准场对解算结果进行了数值验证,证明了本发明的实用性和先进性,具有较高的应用价值,可广泛用于大地物理测量领域。

Description

一种无奇异性顾及局域保障条件计算外部重力异常的方法
技术领域
本发明属于地理测量领域,尤其是一种无奇异性顾及局域保障条件计算外部重力异常的方法。
背景技术
外部重力异常是为精密确定运载火箭、人造卫星、宇宙飞船、导弹武器、航天飞机等航天器飞行轨迹提供重力异常场补偿的关键要件,也是为近地空间开展重力测量数据质量评估提供比对基准必不可少的基础信息。在计算外部重力异常的积分式中,当计算点趋近于观测数据格网点时积分核函数会发生奇异,导致利用积分式计算超低空外部重力异常时,会出现由核函数奇异引起的奇异性问题,无法保证计算结果的有效性。同时,在实际应用中受观测数据覆盖范围限制,无法做到全球覆盖,实际计算过程中外部重力异常的全球积分式需要进行适用观测数据保障条件的改化,以保证计算结果的可靠性。但是目前对于重力异常尚未出现精度较高的计算方法。
发明内容
本发明的目的在于克服现有技术的不足,提出一种无奇异性顾及局域保障条件计算外部重力异常的方法,能够有效的对重力异常进行高精度计算。
本发明解决其技术问题是采取以下技术方案实现的:
一种无奇异性顾及局域保障条件计算外部重力异常的方法,包括以下步骤:
步骤1、计算外部重力异常Δg;
步骤2、采用移去-恢复技术,从流动点重力异常中移去外部空间计算点在球面上的投影点处的重力异常;同时通过积分恒等变换,恢复外部空间计算点在球面上的投影点处的重力异常的全球影响,消除积分核函数奇异性的影响;
步骤3、顾及实测数据局域保障条件,引入全球重力场位模型,利用移去-恢复技术,恢复参考外部重力异常;
步骤4、利用积分恒等式转换,补偿由全球积分过渡到局域积分引起的模型误差,得到计算点处高精度的外部重力异常。
而且,所述步骤3包括以下步骤:
步骤3.1、顾及实测数据局域保障条件,引入全球重力场位模型,利用移去-恢复技术,从重力异常观测值中移去参考重力异常,得到残差重力异常;
步骤3.2、从积分式核函数中移去与参考场对应阶次的核函数球谐表达式,得到截断核函数;
步骤3.3、利用全球重力场位模型高阶信息进行远区效应补偿,以削弱远区截断误差的影响。
而且,所述步骤1的具体实现方法为:
Figure BDA0002942348010000021
Figure BDA0002942348010000022
其中,Δg为外部空间计算点
Figure BDA0002942348010000023
重力异常;Δgq为球面上流动点
Figure BDA0002942348010000024
处的已知观测重力异常;R为地球椭球平均半径;r为计算点地心向径;
Figure BDA0002942348010000025
为计算点的纬度和经度;
Figure BDA0002942348010000026
为流动点的纬度和经度;σ为单位球面;dσ为单位球面的面积元;ψ为计算点至流动点之间的球面角距;
Figure BDA0002942348010000027
l为计算点至流动点之间的空间距离;K(r,ψ)为积分核函数。
而且,所述步骤2的具体实现方法为:
将外部重力异常Δg转换为;
Figure BDA0002942348010000028
其中,ΔgRp为外部空间计算点p在球面上的投影点PR处的重力异常。
而且,所述步骤3.1的具体实现方法为:
将外部重力异常Δg转换为:
Figure BDA0002942348010000029
计算残差重力异常δΔgq和δΔgRp
Figure BDA00029423480100000210
计算参考重力异常Δgref
Figure BDA0002942348010000031
其中,δΔgq为Δgq的残差重力异常;δΔgRp为ΔgRp的残差重力异常;KWG(r,ψ)为截断核函数;
Figure BDA0002942348010000032
为远区效应计算值;Δgref为参考外部重力异常;Δgqref和ΔgRpref分别为与Δgq和ΔgRp对应的由N阶参考场位模型计算得到的参考重力异常;GM为地球引力常数;N为由位模型定义的参考场最高阶数;
Figure BDA0002942348010000033
为完全规格化缔合勒让德函数;
Figure BDA0002942348010000034
Figure BDA0002942348010000035
为完全规格化地球位系数。
而且,所述步骤3.2的具体实现方法为:
Figure BDA0002942348010000036
其中,Pn(cosψ)为n阶勒让德函数。
而且,所述步骤3.3的具体实现方法为:
计算远区效应计算值
Figure BDA0002942348010000037
Figure BDA0002942348010000038
Figure BDA0002942348010000039
Figure BDA00029423480100000310
Figure BDA00029423480100000311
其中,Qn(Δg)为外部重力异常积分核截断系数;Tn为地球扰动位n阶拉普拉斯面球谐函数;Rn,m0)为勒让德函数的远区积分和。
而且,所述步骤4的具体实现方法为:
将外部重力异常Δg转换为:
Figure BDA00029423480100000312
计算δΔgRp在积分远区(σ-σ0)对外部重力异常Δg的影响
Figure BDA00029423480100000313
Figure BDA0002942348010000041
Figure BDA0002942348010000042
Figure BDA0002942348010000043
本发明的优点和积极效果是:
本发明针对计算外部重力异常的奇异性问题,通过采用移去-恢复技术,从流动点重力异常中移去外部空间计算点在球面上的投影点处的重力异常;同时通过积分恒等变换,恢复外部空间计算点在球面上的投影点处的重力异常的全球影响,消除积分核函数奇异性的影响;顾及实测数据局域保障条件,引入全球重力场位模型,利用移去-恢复技术,恢复参考外部重力异常;利用积分恒等式转换,补偿由全球积分过渡到局域积分引起的模型误差,得到计算点处高精度的外部重力异常。本发明利用全球重力场位模型建立的模拟标准场对本发明的解算结果进行了数值验证,证明了本发明新方法的实用性和先进性,具有较高的应用价值,可广泛用于大地物理测量领域。
具体实施方式
以下结合实施例对本发明做进一步详述。
一种无奇异性顾及局域保障条件计算外部重力异常的方法,包括以下步骤:
步骤1、计算外部重力异常Δg:
Figure BDA0002942348010000044
Figure BDA0002942348010000045
其中,Δg为外部空间计算点
Figure BDA0002942348010000046
重力异常;Δgq为球面上流动点
Figure BDA0002942348010000047
处的已知观测重力异常;R为地球椭球平均半径;r为计算点地心向径;
Figure BDA0002942348010000048
为计算点的纬度和经度;
Figure BDA0002942348010000049
为流动点的纬度和经度;σ为单位球面;dσ为单位球面的面积元;ψ为计算点至流动点之间的球面角距;
Figure BDA00029423480100000410
l为计算点至流动点之间的空间距离;K(r,ψ)为积分核函数。
步骤2、在步骤1中,当计算点趋近于数据点时,即当r→R和ψ→0时,会出现分母项l→0,积分核函数K(r,ψ)发生奇异,这说明利用步骤1计算地球外部超低空重力异常时,会出现由核函数奇异性引起的不确定性问题,无法保证计算结果的可靠性和有效性。为了消除步骤1的奇异性,确保外部重力异常的整体连续性,采用移去-恢复技术,从流动点重力异常中移去外部空间计算点在球面上的投影点处的重力异常;同时通过积分恒等变换,恢复外部空间计算点在球面上的投影点处的重力异常的全球影响,消除积分核函数奇异性的影响:
将外部重力异常Δg转换为;
Figure BDA0002942348010000051
其中,ΔgRp为外部空间计算点p在球面上的投影点PR处的重力异常。
步骤3、计算外部重力异常要求全球积分,但在实际应用中受观测数据覆盖范围限制,无法做到全球覆盖,实际计算过程中外部重力异常的全球积分式需要进行适用观测数据保障条件的改化,以保证计算结果的可靠性。顾及实测数据局域保障条件,引入全球重力场位模型,利用移去-恢复技术,恢复参考外部重力异常。
本步骤包括以下步骤:
步骤3.1、顾及实测数据局域保障条件,引入全球重力场位模型,利用移去-恢复技术,从重力异常观测值中移去参考重力异常,得到残差重力异常;
将外部重力异常Δg转换为:
Figure BDA0002942348010000052
计算残差重力异常δΔgq和δΔgRp
Figure BDA0002942348010000053
计算参考重力异常Δgref
Figure BDA0002942348010000054
其中,δΔgq为Δgq的残差重力异常;δΔgRp为ΔgRp的残差重力异常;KWG(r,ψ)为截断核函数;
Figure BDA0002942348010000055
为远区效应计算值;Δgref为参考外部重力异常;Δgqref和ΔgRpref分别为与Δgq和ΔgRp对应的由N阶参考场位模型计算得到的参考重力异常;GM为地球引力常数;N为由位模型定义的参考场最高阶数;
Figure BDA0002942348010000056
为完全规格化缔合勒让德函数;
Figure BDA0002942348010000057
Figure BDA0002942348010000058
为完全规格化地球位系数。
步骤3.2、从积分式核函数中移去与参考场对应阶次的核函数球谐表达式,得到截断核函数,使之与残差重力异常的频谱匹配:
Figure BDA0002942348010000061
其中,Pn(cosψ)为n阶勒让德函数。
步骤3.3、利用全球重力场位模型高阶信息进行远区效应补偿,以削弱远区截断误差的影响
计算远区效应计算值
Figure BDA0002942348010000062
Figure BDA0002942348010000063
Figure BDA0002942348010000064
Figure BDA0002942348010000065
Figure BDA0002942348010000066
其中,Qn(Δg)为外部重力异常积分核截断系数;Tn为地球扰动位n阶拉普拉斯面球谐函数;Rn,m0)为勒让德函数的远区积分和。
步骤4、利用积分恒等式转换,补偿由全球积分过渡到局域积分引起的模型误差,得到计算点处高精度的外部重力异常:
在既消除积分奇异性又顾及实测数据局域保障条件的
Figure BDA0002942348010000067
中,由于远区效应
Figure BDA0002942348010000068
只代表积分项δΔgq在远区的补偿,并未顾及另一积分项δΔgRp在远区的影响。对此利用积分恒等式转换,补偿由全球积分过渡到局域积分引起的模型误差,得到计算点处高精度的外部重力异常。将外部重力异常Δg转换为:
Figure BDA0002942348010000069
计算δΔgRp在积分远区(σ-σ0)对外部重力异常Δg的影响
Figure BDA00029423480100000610
Figure BDA0002942348010000071
Figure BDA0002942348010000072
Figure BDA0002942348010000073
根据上述一种无奇异性顾及局域保障条件计算外部重力异常的方法,采用全球重力场位模型EGM2008构建数值计算检验的标准场,即模拟产生地球表面1′×1′网格重力异常观测量“真值”,同时产生外部设定高度的重力异常理论“真值”。为了体现检验结果的代表性,特意选取重力异常场变化比较剧烈的马里亚纳海沟作为试验区,具体覆盖范围为:
Figure BDA0002942348010000074
Figure BDA0002942348010000075
选取r=R+h,R=6371km。
采用标准场模型EGM2008分别计算标准场7个高度面上的1′×1′网格重力异常“真值”Δg′ti(i=1,2,…,7),每一个高度面对应360×360=129600个网格点数据,9个高度分别为:hi=0km,0.1km,0.3km,1km,3km,5km,10km,30km,50km。如表1所示为其中的5个高度面上的重力异常理论“真值”统计结果。
表1由EGM2008模型计算得到的重力异常计结果/mGal
高度/km 最大值 最小值 平均值 均方根值
0 138.85 -78.48 0.38 22.63
1 122.29 -69.92 0.35 20.68
3 97.73 -56.82 0.29 17.40
30 10.96 -7.30 0.03 2.43
50 2.68 -2.05 0.01 0.66
为了对比分析本发明的有效性,引入外部重力异常Δg传统同时计算9个高度面上的外部重力异常。将计算值分别与相对应的理论“真值”Δgtri作比较,可获得精度评估信息,具体比对结果如表2所示。这里积分半径统一取为ψ0=2°,为了减小积分边缘效应对评估结果的影响,表2只列出中心区2°×2°方块内的比对结果。
表2由本发明计算的外部重力异常与“真值”的比较(单位:mGal)
Figure BDA0002942348010000076
Figure BDA0002942348010000081
对比表2和表1结果可以看出,传统在超低空高度段的均方根量值甚至超过了重力异常自身大小,说明传统在超低空高度段是失效的;本发明的计算结果与比对基准“真值”的最大互差均不超过1.0mGal,均方根不超过0.3mGal,验证了本发明的必要性和有效性,表明新方法实用易行,具有较高的应用价值。
需要强调的是,本发明所述的实施例是说明性的,而不是限定性的,因此本发明包括并不限于具体实施方式中所述的实施例,凡是由本领域技术人员根据本发明的技术方案得出的其他实施方式,同样属于本发明保护的范围。

Claims (8)

1.一种无奇异性顾及局域保障条件计算外部重力异常的方法,其特征在于:包括以下步骤:
步骤1、计算外部重力异常Δg;
步骤2、采用移去-恢复技术,从流动点重力异常中移去外部空间计算点在球面上的投影点处的重力异常;同时通过积分恒等变换,恢复外部空间计算点在球面上的投影点处的重力异常的全球影响,消除积分核函数奇异性的影响;
步骤3、顾及实测数据局域保障条件,引入全球重力场位模型,利用移去-恢复技术,恢复参考外部重力异常;
步骤4、利用积分恒等式转换,补偿由全球积分过渡到局域积分引起的模型误差,得到计算点处高精度的外部重力异常。
2.根据权利要求1所述的一种无奇异性顾及局域保障条件计算外部重力异常的方法,其特征在于:所述步骤3包括以下步骤:
步骤3.1、顾及实测数据局域保障条件,引入全球重力场位模型,利用移去-恢复技术,从重力异常观测值中移去参考重力异常,得到残差重力异常;
步骤3.2、从积分式核函数中移去与参考场对应阶次的核函数球谐表达式,得到截断核函数;
步骤3.3、利用全球重力场位模型高阶信息进行远区效应补偿,以削弱远区截断误差的影响。
3.根据权利要求1或2所述的一种无奇异性顾及局域保障条件计算外部重力异常的方法,其特征在于:所述步骤1的具体实现方法为:
Figure FDA0003531729780000011
Figure FDA0003531729780000012
其中,Δg为外部空间计算点
Figure FDA0003531729780000013
重力异常;Δgq为球面上流动点
Figure FDA0003531729780000014
处的已知观测重力异常;R为地球椭球平均半径;r为计算点地心向径;
Figure FDA0003531729780000015
为计算点的纬度和经度;
Figure FDA0003531729780000016
为流动点的纬度和经度;σ为单位球面;dσ为单位球面的面积元;ψ为计算点至流动点之间的球面角距;
Figure FDA0003531729780000017
l为计算点至流动点之间的空间距离;K(r,ψ)为积分核函数。
4.根据权利要求3所述的一种无奇异性顾及局域保障条件计算外部重力异常的方法,其特征在于:所述步骤2的具体实现方法为:
将外部重力异常Δg转换为;
Figure FDA0003531729780000021
其中,ΔgRp为外部空间计算点p在球面上的投影点PR处的重力异常。
5.根据权利要求2所述的一种无奇异性顾及局域保障条件计算外部重力异常的方法,其特征在于:所述步骤3.1的具体实现方法为:
将外部重力异常Δg转换为:
Figure FDA0003531729780000022
计算残差重力异常δΔgq和δΔgRp
Figure FDA0003531729780000023
计算参考重力异常Δgref
Figure FDA0003531729780000024
其中,δΔgq为Δgq的残差重力异常;δΔgRp为ΔgRp的残差重力异常;KWG(r,ψ)为截断核函数;Δgrq(σ-σ0)为远区效应计算值;Δgref为参考外部重力异常;Δgqref和ΔgRpref分别为与Δgq和ΔgRp对应的由N阶参考场位模型计算得到的参考重力异常;GM为地球引力常数;N为由位模型定义的参考场最高阶数;
Figure FDA0003531729780000025
为完全规格化缔合勒让德函数;
Figure FDA0003531729780000026
Figure FDA0003531729780000027
为完全规格化地球位系数。
6.根据权利要求2所述的一种无奇异性顾及局域保障条件计算外部重力异常的方法,其特征在于:所述步骤3.2的具体实现方法为:
Figure FDA0003531729780000028
其中,Pn(cosψ)为n阶勒让德函数。
7.根据权利要求2所述的一种无奇异性顾及局域保障条件计算外部重力异常的方法,其特征在于:所述步骤3.3的具体实现方法为:
计算远区效应计算值
Figure FDA0003531729780000031
Figure FDA0003531729780000032
Figure FDA0003531729780000033
Figure FDA0003531729780000034
Figure FDA0003531729780000035
其中,Qn(Δg)为外部重力异常积分核截断系数;Tn为地球扰动位n阶拉普拉斯面球谐函数;Rn,m0)为勒让德函数的远区积分和。
8.根据权利要求1所述的一种无奇异性顾及局域保障条件计算外部重力异常的方法,其特征在于:所述步骤4的具体实现方法为:
将外部重力异常Δg转换为:
Figure FDA0003531729780000036
计算δΔgRp在积分远区(σ-σ0)对外部重力异常Δg的影响
Figure FDA0003531729780000037
Figure FDA0003531729780000038
Figure FDA0003531729780000039
Figure FDA00035317297800000310
CN202110181838.5A 2021-02-08 2021-02-08 一种无奇异性顾及局域保障条件计算外部重力异常的方法 Active CN112965128B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110181838.5A CN112965128B (zh) 2021-02-08 2021-02-08 一种无奇异性顾及局域保障条件计算外部重力异常的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110181838.5A CN112965128B (zh) 2021-02-08 2021-02-08 一种无奇异性顾及局域保障条件计算外部重力异常的方法

Publications (2)

Publication Number Publication Date
CN112965128A CN112965128A (zh) 2021-06-15
CN112965128B true CN112965128B (zh) 2022-05-24

Family

ID=76284697

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110181838.5A Active CN112965128B (zh) 2021-02-08 2021-02-08 一种无奇异性顾及局域保障条件计算外部重力异常的方法

Country Status (1)

Country Link
CN (1) CN112965128B (zh)

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103513294B (zh) * 2013-09-29 2016-05-18 清华大学 一种低低星星跟踪卫星重力场测量性能解析计算方法
US9778360B2 (en) * 2013-12-17 2017-10-03 Fugro N.V. Method and system for generating a geoid via three computation spaces and airborne-acquired gravity data
CN104834320A (zh) * 2015-03-23 2015-08-12 中国人民解放军信息工程大学 一种空间分层扰动引力场网格模型快速构建方法
CN108376187B (zh) * 2018-01-19 2021-09-10 中国人民解放军92859部队 一种海域流动点外部扰动引力垂向分量的无奇异计算方法

Also Published As

Publication number Publication date
CN112965128A (zh) 2021-06-15

Similar Documents

Publication Publication Date Title
US11313698B2 (en) Method for initial alignment of radar assisted airborne strapdown inertial navigation system
CN108332649B (zh) 一种滑坡形变综合预警方法及系统
CN108317993B (zh) 一种集成gnss和激光跟踪仪的垂线偏差测量装置及方法
CN112965124B (zh) 一种顾及局域保障条件计算外部重力异常垂直梯度的方法
CN110187369B (zh) 一种基于gnss卫星位置观测的垂线偏差测量和验证方法
CN109613583A (zh) 基于单星与地面站测向及联合测时差的无源目标定位方法
CN116182795B (zh) 普速铁路纵断面精密测量方法
Novák et al. Ellipsoidal representation of the topographical potential and its vertical gradient
CN108875174A (zh) 一种基于多段打靶法的不变拟周期轨道确定方法
CN109581523B (zh) 一种采用卫星跟踪卫星装置对加速度计校准的方法和系统
CN112965128B (zh) 一种无奇异性顾及局域保障条件计算外部重力异常的方法
Sośnica LAGEOS sensitivity to ocean tides
CN112965127B (zh) 一种基于重力异常计算外部扰动重力径向分量的方法
CN114353744A (zh) 一种基于srtm1和egm2008的gps高程拟合方法及拟合系统
CN112949049B (zh) 一种利用带限思想计算重力异常低阶径向导数的方法
Simonetti et al. Robust modeling of geodetic altitude from barometric altimeter and weather data
CN112818285B (zh) 一种计算外部扰动重力北向分量中央区效应的方法
CN112965125B (zh) 一种基于重力异常计算外部扰动重力东向分量的方法
CN113324539A (zh) 一种sins/srs/cns多源融合自主组合导航方法
CN114924270A (zh) 基于GNSS的InSAR形变监测基准建立方法及装置
CN112965123B (zh) 一种基于重力异常计算外部扰动重力北向分量的方法
CN112987118B (zh) 一种利用带限思想计算重力异常高阶径向导数的方法
CN112965126B (zh) 一种计算外部扰动重力东向分量中央区效应的方法
CN112836378B (zh) 基于Poisson理论计算外部重力异常垂直梯度中央区效应的方法
Tao et al. Research on Accurate Impact Point Prediction Methods of Heavy Rocket Debris

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