CN117233799A - 基于虚拟基准站的煤矿采空区地表形变监测方法 - Google Patents
基于虚拟基准站的煤矿采空区地表形变监测方法 Download PDFInfo
- Publication number
- CN117233799A CN117233799A CN202311479876.4A CN202311479876A CN117233799A CN 117233799 A CN117233799 A CN 117233799A CN 202311479876 A CN202311479876 A CN 202311479876A CN 117233799 A CN117233799 A CN 117233799A
- Authority
- CN
- China
- Prior art keywords
- station
- delay
- difference
- monitoring
- satellite
- 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.)
- Granted
Links
- 238000012544 monitoring process Methods 0.000 title claims abstract description 123
- 238000000034 method Methods 0.000 title claims abstract description 56
- 239000003245 coal Substances 0.000 title claims abstract description 23
- 238000012937 correction Methods 0.000 claims abstract description 56
- 239000005436 troposphere Substances 0.000 claims abstract description 48
- 208000028257 Joubert syndrome with oculorenal defect Diseases 0.000 claims abstract description 28
- 238000004458 analytical method Methods 0.000 claims abstract description 22
- 238000001914 filtration Methods 0.000 claims abstract description 19
- 239000011159 matrix material Substances 0.000 claims description 29
- 239000005433 ionosphere Substances 0.000 claims description 26
- 238000004364 calculation method Methods 0.000 claims description 19
- 239000013598 vector Substances 0.000 claims description 16
- 230000008569 process Effects 0.000 claims description 10
- 230000009466 transformation Effects 0.000 claims description 7
- 230000009977 dual effect Effects 0.000 claims description 6
- 238000013441 quality evaluation Methods 0.000 claims description 6
- 238000012916 structural analysis Methods 0.000 claims description 6
- 238000007667 floating Methods 0.000 claims description 5
- 238000013461 design Methods 0.000 claims description 4
- 230000000694 effects Effects 0.000 claims description 4
- 230000036541 health Effects 0.000 claims description 3
- 230000007613 environmental effect Effects 0.000 claims description 2
- 238000013507 mapping Methods 0.000 claims description 2
- 238000012216 screening Methods 0.000 claims description 2
- 230000026676 system process Effects 0.000 claims description 2
- 230000007704 transition Effects 0.000 claims description 2
- 238000004891 communication Methods 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 3
- 230000000903 blocking effect Effects 0.000 description 2
- 238000003908 quality control method Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 238000012806 monitoring device Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 230000001360 synchronised effect Effects 0.000 description 1
Classifications
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
一种基于虚拟基准站的煤矿采空区地表形变监测方法。包括：获取CORS站点的GNSS数据、气象数据和坐标数据；对监测区域进行格网化，并计算对流层的干、湿延迟量；利用CORS站点数据组建双差观测方程，并求解CORS站点的整周模糊度；通过服务器，构建区域大气增强模型并播发格网改正数；监测站点使用格网改正数和概略坐标内插各种延迟量和误差，服务器则生成虚拟基准站观测值；监测站点利用这些观测值组建短基线双差观测方程，并应用卡尔曼滤波模型进行高精度准实时滤波；根据监测站点坐标分析形变，并根据分析结果进行安全预警。解决了现有技术不能满足数以千计乃至数以万计的并发监测请求的问题。
Description
技术领域
本申请涉及卫星定位技术领域，具体涉及一种基于虚拟基准站的煤矿采空区地表形变监测方法。
背景技术
为了掌握煤矿开采活动导致地表移动变形影响规律，需要对其开展监测，现有的地表变形监测技术主要采用传统的网络RTK进行监测，这种监测技术存在的问题是：用户并发规模只有几十到几百并发请求，不能满足数以千计乃至数以万计的并发监测请求。并且煤矿采空区一般面积较大，监测设备较多，且采空区多位于山区，通讯较为困难，对于足数以千计乃至数以万计的并发监测更是难上加难。
发明内容
为解决上述技术问题，本申请提供一种基于虚拟基准站的煤矿采空区地表形变监测方法。
第一方面，本申请实施例提供一种基于虚拟基准站的煤矿采空区地表形变监测方法，所述基于虚拟基准站的煤矿采空区地表形变监测方法包括：
步骤S01、获取CORS站的GNSS实时观测数据、气象观测数据和先验坐标数据，对获取的GNSS实时观测数据进行质量评估，过滤出不满足质量标准的GNSS实时观测数据，质量评估的指标包括高度角、周跳、多路径效应、信噪比和卫星健康状态；
步骤S02、利用CORS站的气象观测数据和对流层干延迟模型计算该站处的对流层干延迟量；
步骤S03、将过滤后的GNSS实时观测值组建双差观测方程，其中，利用先验坐标数据先验已知的特性，将对流层湿延迟、电离层延迟残余量、综合误差作为未知参数组建双差观测方程；
步骤S04、对每历元每颗卫星进行法方程叠加，构建整体法方程，求解CORS站的整周模糊度；
步骤S05、求解出每个CORS站处对流层延迟量、电离层延迟量和综合误差；
步骤S06、通过服务器构建区域大气增强模型，对服务区域进行格网化改造，并通过服务器播发区域大气模型格网改正数；
步骤S07、监测站根据监测站概略坐标和接收到的格网改正数，内插出监测站处的对流层延迟量、电离层延迟量和综合误差；
步骤S08、根据监测站概略坐标和步骤S07的内插结果，生成虚拟基准站观测值；
步骤S09、监测站将虚拟基准站观测值和监测站实际观测值组建短基线双差观测方程；
步骤S10、采用卡尔曼滤波模型对短基线双差观测方程进行卡尔曼滤波，得到监测站坐标时间序列；
步骤S11、对监测站坐标时间序列进行建模分析，提取形变量、形变方向、形变速率和形变原因分析，并分析变形的水平移动、沉降量、水平、倾斜和曲率，得到监测站变形分析结果；
步骤S12、根据监测站变形分析结果进行安全预警判断，若达到预警值则进行预警，若未达到预警值则输出变形分析结果。
本申请实施例提供的技术方案带来的有益效果包括：
通过将GNSS实时观测数据、气象观测数据和先验坐标数据的接收，以及对流层延迟量、电离层延迟量和综合误差的解算，都放在服务器处进行，通过格网化的办法内插计算格网的改正数，仅向监测站发送少量数据的改正数，然后监测站基于自身概略坐标和服务器播发的改正数，通过内插获得监测站处的对流层延迟量、电离层延迟量和综合误差，并生成虚拟基准站观测值，然后使用卡尔曼滤波获取监测站点坐标时间序列。以上办法使得大量数据的收发均是通过服务器完成，例如GNSS实时观测数据和气象观测数据的收发，仅需要在监测站处收发少量改正数和定位坐标，大大的减轻了单个监测站的数据收发量，从而使得同样的带宽下，可容纳数以千计乃至数以万计的并发监测请求；
解算服务通过服务器进行，解算复杂度与用户规模无关，播发服务与解算服务分离，解决用户请求与解算耦合问题；
通过卡尔曼滤波获取监测站点坐标时间序列，利用时间序列可进行形变量、形变方向和形变原因分析；
无需投资建设基准站，节约了大量的基准站建设投资，避免了基准站本身受采空区沉降影响导致的基准不稳定问题；
实现了对采空区地表形变得到实时毫米级监测，平面精度1~3mm，高程精度3~5mm，监测结果精度更稳定，更可靠；
避免了使用基准站监测地表形变缺少统一的起算基准，一旦基准站变动，前后两期数据不连续的问题。
附图说明
图1为本申请基于虚拟基准站的煤矿采空区地表形变监测方法一实施例的流程示意图。
具体实施方式
为了使本技术领域的人员更好地理解本申请方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。
为使本申请的目的、技术方案和优点更加清楚，下面将结合附图对本申请实施方式作进一步地详细描述。
第一方面，本申请实施例提供一种基于虚拟基准站的煤矿采空区地表形变监测方法。
一实施例中，参照图1，图1为本申请基于虚拟基准站的煤矿采空区地表形变监测方法一实施例的流程示意图。如图1所示，基于虚拟基准站的煤矿采空区地表形变监测方法包括：
步骤S01、获取CORS站的GNSS实时观测数据、气象观测数据和先验坐标数据，对获取的GNSS实时观测数据进行质量评估，过滤出不满足质量标准的GNSS实时观测数据，质量评估的指标包括高度角、周跳、多路径效应、信噪比和卫星健康状态；
步骤S02、利用CORS站的气象观测数据和对流层干延迟模型计算该站处的对流层干延迟量；
步骤S03、将过滤后的GNSS实时观测值组建双差观测方程，其中，利用先验坐标数据先验已知的特性，将对流层湿延迟、电离层延迟残余量、综合误差作为未知参数组建双差观测方程；
步骤S04、对每历元每颗卫星进行法方程叠加，构建整体法方程，求解CORS站的整周模糊度；
步骤S05、求解出每个CORS站处对流层延迟量、电离层延迟量和综合误差；
步骤S06、通过服务器构建区域大气增强模型，对服务区域进行格网化改造，并通过服务器播发区域大气模型格网改正数；
步骤S07、监测站根据监测站概略坐标和接收到的格网改正数，内插出监测站处的对流层延迟量、电离层延迟量和综合误差；
步骤S08、根据监测站概略坐标和步骤S07的内插结果，生成虚拟基准站观测值；
步骤S09、监测站将虚拟基准站观测值和监测站实际观测值组建短基线双差观测方程；
步骤S10、采用卡尔曼滤波模型对短基线双差观测方程进行卡尔曼滤波，得到监测站坐标时间序列；
步骤S11、对监测站坐标时间序列进行建模分析，提取形变量、形变方向、形变速率和形变原因分析，并分析变形的水平移动、沉降量、水平、倾斜和曲率，得到监测站变形分析结果；
步骤S12、根据监测站变形分析结果进行安全预警判断，若达到预警值则进行预警，若未达到预警值则输出变形分析结果。
监测站与CORS站点和服务器间的通信可采用非阻塞IO技术，避免了传统TCP通信IO阻塞问题，解决大规模数据并发接入。并且采用异步、事件驱动的并发网络框架，实现Ntrip、TCP协议北斗监测站同步请求。
进一步地，一实施例中，步骤S03的双差观测方程为：
其中，b表示基准站，r表示监测站，i表示参考卫星，j表示非参考卫星，表示站间星间差分算子，/>表示以周为单位的载波观测值，/>表示对应频率的波长，P表示以m为单位的伪距观测值，/>表示卫星与接收机之间的几何距离，T表示对流层延迟误差，I表示电离层延迟误差，N表示对应频率的整周模糊度，/>和/>分别表示载波观测值和伪距观测值的噪声，表示基准站b与监测站r之间，参考卫星i与非参考卫星 j 之间的载波观测值，/>表示基准站b与监测站r之间，参考卫星i与非参考卫星 j 之间的几何距离，/>表示基准站b与监测站r之间，参考卫星i与非参考卫星 j 之间的对流层延迟误差，/>表示基准站b与监测站r之间，参考卫星i与非参考卫星 j 之间的电离层延迟误差，/>表示基准站b与监测站r之间，参考卫星i与非参考卫星 j 之间的整周模糊度，/>表示基准站b与监测站r之间，参考卫星i与非参考卫星 j 之间的载波观测值的噪声，/>表示基准站b与监测站r之间，参考卫星i与非参考卫星 j 之间的伪距观测值，/>表示基准站b与监测站r之间，参考卫星i与非参考卫星 j 之间的伪距观测值的噪声。
进一步地，一实施例中，步骤S04中的整体法方程为：
其中，k表示一个历元；n表示观测时段内的历元总数；表示历元k的设计矩阵Bk的转置；/>表示观测历元的权阵。
进一步地，一实施例中，求解CORS站的整周模糊度包括：
基于预设选星规则分别从BDS系统和GPS系统中筛选出参考星；
对BDS系统中的卫星和GPS系统中的卫星按照波长从长到短的顺序分别固定选择观测值组合的双差模糊度，对于BDS系统中的卫星，先固定波长较长的观测值组合的双差模糊度，再固定波长较短的观测值组合的双差模糊度，依次固定超宽巷、宽巷和基本频点观测值3种不同波长的观测值；
对于GPS系统中的卫星，先固定宽巷双差模糊度，然后将固定的宽巷双差模糊度代入到无电离层组合，并利用卡尔曼滤波估计L1载波相位的整周模糊度和相对对流层天顶延迟，得到双差模糊度浮点解和方差阵，最后利用改进的LAMBDA方法实时搜索双差整周模糊度；
按照基线从短到长的顺序分别固定选择基线；
根据固定结果进行法方程更新，以得到模糊度参数实数解；
根据模糊度固定规则将模糊度参数实数解解固定模糊度整数解；
其中，所述根据固定结果进行法方程更新，以得到模糊度参数实数解包括以下步骤：
求解出所述双差模糊度的实数解；
将所述实数解中的模糊度参数固定为整数，并进行检验，得到整周模糊度；
将所述整周模糊度作为已知值代入法方程进行更新；
按上述过程进行迭代，直至得到所有整周模糊度参数；
所述预设选星规则为：
GPS和BDS每个系统各选一颗卫星作为参考星，在BDS系统中，首选GEO卫星作为参考星，若由于环境遮挡任何一颗GEO卫星均不为参考星，则依次选择IGSO、MEO卫星为参考星；在GPS卫星中，选择高度角较高的卫星作为参考星。
其中，基站间整周模糊度解算：
其一、通过载波观测值和P码伪距组成Melbourne-Wübbena组合（M-W组合），计算得到宽巷整周模糊度。其二、将载波相位和伪距观测值进行无电离层组合，解算得到无电离层模糊度的实数解。其三、恢复L1、L2双差模糊度浮点解，并固定。
通过Melbourne-Wübbena组合进行宽巷整周模糊度的解算，其计算公式为：
式中，为Melbourne-Wübbena组合观测值，/>为宽巷整周模糊度，/>为L1和L2波段的频率，c为真空环境下的光速。
无电离层延迟组合模糊度是L1和L2模糊度的线性组合，其计算公式如下：
利用固定的宽巷模糊度和无电离层组合模糊度的浮点解及协方差，恢复L1频点的双差模糊度浮点解及其协方差：
根据L1的双差模糊度是具有整周特性的，利用经典的LAMBDA方法对模糊度进行降相关处理，从而大大减少计算量。LAMBDA算法大致流程为：a)对模糊度方差阵使用整数Gauss变换，构建整数可逆矩阵Z，使用Z矩阵对模糊度进行Z变换，对模糊度相应的协方差阵进行降相关Z变换；b)根据降相关的模糊度协方差阵，对Z变换后的模糊度进行搜索；c)使用Z-1还原相应的模糊度。
需要进一步进行模糊度的质量控制。当前大部分的模糊度质量控制都是基于假设检验，给定显著性水平，检验最小的和次小的/>是否有显著的不同，如果显著相同，则需要继续滤波并在下一历元再对模糊度进行搜索和固定，直到找到最小的/>和次小的/>显著不同为止；如果显著不同，则判断模糊度N固定正确。
一般而言我们利用最小的和次小的/>两者的比值来判断模糊度是否固定正确。
该方法认为两组模糊度估值是相互独立且服从正态分布的，N和N’为这两组模糊度的期望值，和/>为其方差值。其中拒绝域为/>为F分布的上置信水平为/>的分位数，n为参与模糊度滤波的历元数（Frei，1990）。如果统计量ratio落入拒绝域，则认为/>以置信水平/>显著大于/>，所以判断该组模糊度N固定正确。
在L1和宽巷双差模糊度都成功固定的情况下，L2的双差模糊度也可以根据L1、L2和宽巷模糊度之间的线性关系轻松的计算得到：
。
进一步地，一实施例中，步骤S05中对流层延迟量的计算方法为：
式中，为双差对流层延迟干分量，通过经验模型计算得到，/>为天顶对流层湿延迟，/>为测站m和卫星i的高度角，/>为映射函数，/>表示双差对流层湿延迟分量；
所述步骤S05中电离层延迟量的计算方法为：
其中，和/>分别为L1频段和L2频段观测值的频率，/>和/>分别为L1频段和L2频段观测值波长，/>和/>分别为L1频段和L2频段双差载波观测值，/>和/>分别为L1频段和L2频段双差整周模糊度；
所述步骤S05中综合误差的计算方法为：
其中，表示卫星与接收机之间的几何距离，/>表示载波观测值，/>表示整周模糊度，/>表示载波的波长，/>表示电离层延迟量，/>表示对流层延迟量。
进一步地，一实施例中，步骤S07具体包括：
对内插区域内监测站的电离层延迟量的计算方法为：
在中心电离层高度处进行空间结构性分析和空间变异性分析；
将预设高度作为电离层误差改正基准面；
在所述电离层误差改正基准面进行监测站非差电离层延迟误差内插计算，得到监测站处的电离层改正数，根据监测站处的电离层改正数计算得到电离层延迟量；
所述非差电离层延迟误差内插计算包括：
在线性无偏、最优估计条件下，得到电离层总电子含量插值的克里金方程组，所述克里金方程组为：
式中，μ表示拉格朗日乘数因子，dij表示基准站与监测站间的距离，di0表示插值点与观测值间的距离，λj表示加权系数；
对克里金方程组进行求解，得到加权系数；
基于所述加权系数计算得到电离层总电子含量估计值；
通过所述电离层总电子含量估计值计算得到电离层改正数；
对内插区域内监测站的对流层延迟量的计算方法为：
将对流层相关误差分为干延迟部分和湿延迟部分，干延迟部分通过GPT2w或ITG模型进行改正，湿延迟部分根据各基准站解算的对流层湿延迟结果采用基于高度改正的Kriging内插算法，计算出监测站处的对流层改正数，根据监测站处的对流层改正数计算得到对流层延迟量；
所述基于高度改正的Kriging内插算法具体包括：
将各基准站的高程取平均值得到内插区域的对流层湿延迟改正的高程基准面；
根据对流层湿延迟高差改正模型计算出，各基准站到对流层误差改正基准面的对流层湿延迟改正数；
在高程基准面上对各个基准站对流层湿延迟进行空间结构性分析和变异性分析；
监测站根据监测站概略坐标，根据对流层湿延迟的空间结构性分析和变异性分析结果，采用Kriging算法内插出该处对流层湿延迟改正基准面的对流层湿延迟改正数；
根据监测站概略坐标高程，计算出监测站到对流层误差改正基准面的高差；
根据对流层湿延迟高差改正模型计算出监测站概略坐标高程处的对流层湿延迟改正数；
对内插区域内监测站的综合误差的计算方法为：
通过反距离权重法进行内插得到综合误差。
进一步地，一实施例中，步骤S08具体包括：
服务器基于监测站发来的概略坐标，构建虚拟基准站，通过内插算法获得虚拟基准站相对于CORS站的电离层、对流层和几何距离的相对差值，利用虚拟基准站和CORS站的坐标和卫星星历，将CORS站的观测值归算至虚拟基准站，然后再顾及对流层延迟误差和电离层延迟误差的站间相对差值，获得虚拟基准站的观测值；
计算CORS站的对流层延迟误差双差后的差值，计算方法如下：
计算CORS站处的电离层延迟误差双差后的差值，计算方法如下：
其中，和/>分别为L1频段和L2频段观测值的频率，/>和/>分别为L1频段和L2频段双差载波观测值，/>和/>分别为L1频段和L2频段双差整周模糊度，/>表示卫星与接收机之间的几何距离；
通过空间插值算法计算出虚拟基准站处的对流层延迟误差双差后的差值和电离层延迟误差双差后的差值/>；
计算虚拟基准站的观测值。
进一步地，一实施例中，所述虚拟基准站观测值包括虚拟基准站载波观测值和虚拟基准站伪距观测值，所述虚拟基准站载波观测值的计算方法如下：
其中，表示虚拟基准站载波观测值双差误差改正数，/>包括频点1和频点2分别对应的载波观测值双差误差改正数/>和/>，计算方法如下：
所述虚拟基准站伪距观测值的计算方法如下：
其中，表示虚拟基准站伪距双差误差改正数，/>包括频点1和频点2分别对应的虚拟基准站伪距双差误差改正数/>和/>，计算方法如下：
其中，下标v表示虚拟基准站。
进一步地，一实施例中，步骤S10具体包括：
预测状态向量及方差-协方差阵，其中，预测状态向量：
预测方差-协方差阵：
其中，表示/>时刻的坐标向量的估计值，/>表示/>时刻方差-协方差阵的估计值，/>表示/>时刻坐标向量预测值，/>表示/>时刻方差-协方差阵的估计值，/>为/>时刻至/>时刻的系统状态转移矩阵，/>是/>时刻系统过程噪声的方差矩阵；
更新状态向量及方差-协方差阵，包括：
计算滤波增益：
更新状态向量：
更新方差-协方差阵：
其中，表示/>时刻观测方程的系数矩阵，/>表示/>时刻观测噪声的方差矩阵，表示/>时刻的坐标观测值向量，I为单位矩阵，初始条件为/>和/>；
将卡尔曼滤波每一步迭代得到的坐标状态向量组成监测站坐标时间序列。
进一步地，一实施例中，计算两相邻点在竖直方向的相对移动与两相邻点间水平距离的比值确定倾斜变形；计算两相邻线段的倾斜差与两线段中点间的水平距离的比值确定曲率变形；计算两相邻点的水平移动差值与两点间水平距离的比值确定水平变形。
需要说明的是，上述本申请实施例序号仅仅为了描述，不代表实施例的优劣。
本申请的说明书和权利要求书及上述附图中的术语“包括”和“具有”以及它们任何变形，意图在于覆盖不排他的包含。例如包含了一系列步骤或单元的过程、方法、系统、产品或设备没有限定于已列出的步骤或单元，而是可选地还包括没有列出的步骤或单元，或可选地还包括对于这些过程、方法、产品或设备固有的其他步骤或单元。术语“第一”、“第二”和“第三”等描述，是用于区分不同的对象等，其不代表先后顺序，也不限定“第一”、“第二”和“第三”是不同的类型。
在本申请实施例的描述中，“示例性的”、“例如”或者“举例来说”等用于表示作例子、例证或说明。本申请实施例中被描述为“示例性的”、“例如”或者“举例来说”的任何实施例或设计方案不应被解释为比其它实施例或设计方案更优选或更具优势。确切而言，使用“示例性的”、“例如”或者“举例来说”等词旨在以具体方式呈现相关概念。
在本申请实施例的描述中，除非另有说明，“/”表示或的意思，例如，A/B可以表示A或B；文本中的“和/或”仅仅是一种描述关联对象的关联关系，表示可以存在三种关系，例如，A和/或B，可以表示：单独存在A，同时存在A和B，单独存在B这三种情况，另外，在本申请实施例的描述中，“多个”是指两个或多于两个。
在本申请实施例描述的一些流程中，包含了按照特定顺序出现的多个操作或步骤，但是应该理解，这些操作或步骤可以不按照其在本申请实施例中出现的顺序来执行或并行执行，操作的序号仅用于区分开各个不同的操作，序号本身不代表任何的执行顺序。另外，这些流程可以包括更多或更少的操作，并且这些操作或步骤可以按顺序执行或并行执行，并且这些操作或步骤可以进行组合。
通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例方法可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件，但很多情况下前者是更佳的实施方式。基于这样的理解，本申请的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品存储在如上所述的一个存储介质（如ROM/RAM、磁碟、光盘）中，包括若干指令用以使得一台终端设备执行本申请各个实施例所述的方法。
以上仅为本申请的优选实施例，并非因此限制本申请的专利范围，凡是利用本申请说明书及附图内容所作的等效结构或等效流程变换，或直接或间接运用在其他相关的技术领域，均同理包括在本申请的专利保护范围内。
Claims (10)
1.一种基于虚拟基准站的煤矿采空区地表形变监测方法，其特征在于，所述基于虚拟基准站的煤矿采空区地表形变监测方法包括：
步骤S01、获取CORS站的GNSS实时观测数据、气象观测数据和先验坐标数据，对获取的GNSS实时观测数据进行质量评估，过滤出不满足质量标准的GNSS实时观测数据，质量评估的指标包括高度角、周跳、多路径效应、信噪比和卫星健康状态；
步骤S02、利用CORS站的气象观测数据和对流层干延迟模型计算该站处的对流层干延迟量；
步骤S03、将过滤后的GNSS实时观测值组建双差观测方程，其中，利用先验坐标数据先验已知的特性，将对流层湿延迟、电离层延迟残余量、综合误差作为未知参数组建双差观测方程；
步骤S04、对每历元每颗卫星进行法方程叠加，构建整体法方程，求解CORS站的整周模糊度；
步骤S05、求解出每个CORS站处对流层延迟量、电离层延迟量和综合误差；
步骤S06、通过服务器构建区域大气增强模型，对服务区域进行格网化改造，并通过服务器播发区域大气模型格网改正数；
步骤S07、监测站根据监测站概略坐标和接收到的格网改正数，内插出监测站处的对流层延迟量、电离层延迟量和综合误差；
步骤S08、根据监测站概略坐标和步骤S07的内插结果，生成虚拟基准站观测值；
步骤S09、监测站将虚拟基准站观测值和监测站实际观测值组建短基线双差观测方程；
步骤S10、采用卡尔曼滤波模型对短基线双差观测方程进行卡尔曼滤波，得到监测站坐标时间序列；
步骤S11、对监测站坐标时间序列进行建模分析，提取形变量、形变方向、形变速率和形变原因分析，并分析变形的水平移动、沉降量、水平、倾斜和曲率，得到监测站变形分析结果；
步骤S12、根据监测站变形分析结果进行安全预警判断，若达到预警值则进行预警，若未达到预警值则输出变形分析结果。
2.如权利要求1所述的基于虚拟基准站的煤矿采空区地表形变监测方法，其特征在于，所述步骤S03的双差观测方程为：
其中，b表示基准站，r表示监测站，i表示参考卫星，j表示非参考卫星，表示站间星间差分算子，/>表示以周为单位的载波观测值，/>表示对应频率的波长，P表示以m为单位的伪距观测值，/>表示卫星与接收机之间的几何距离，T表示对流层延迟误差，I表示电离层延迟误差，N表示对应频率的整周模糊度，/>和/>分别表示载波观测值和伪距观测值的噪声，/>表示基准站b与监测站r之间，参考卫星i与非参考卫星 j 之间的载波观测值，/>表示基准站b与监测站r之间，参考卫星i与非参考卫星 j 之间的几何距离，/>表示基准站b与监测站r之间，参考卫星i与非参考卫星 j 之间的对流层延迟误差，/>表示基准站b与监测站r之间，参考卫星i与非参考卫星 j 之间的电离层延迟误差，/>表示基准站b与监测站r之间，参考卫星i与非参考卫星 j 之间的整周模糊度，/>表示基准站b与监测站r之间，参考卫星i与非参考卫星 j 之间的载波观测值的噪声，/>表示基准站b与监测站r之间，参考卫星i与非参考卫星 j 之间的伪距观测值，/>表示基准站b与监测站r之间，参考卫星i与非参考卫星 j 之间的伪距观测值的噪声。
3.如权利要求1所述的基于虚拟基准站的煤矿采空区地表形变监测方法，其特征在于，步骤S04中的整体法方程为：
其中，k表示一个历元；n表示观测时段内的历元总数；表示历元k的设计矩阵Bk的转置；/>表示观测历元的权阵。
4.如权利要求3所述的基于虚拟基准站的煤矿采空区地表形变监测方法，其特征在于，求解CORS站的整周模糊度包括：
基于预设选星规则分别从BDS系统和GPS系统中筛选出参考星；
对BDS系统中的卫星和GPS系统中的卫星按照波长从长到短的顺序分别固定选择观测值组合的双差模糊度，对于BDS系统中的卫星，先固定波长较长的观测值组合的双差模糊度，再固定波长较短的观测值组合的双差模糊度，依次固定超宽巷、宽巷和基本频点观测值3种不同波长的观测值；
对于GPS系统中的卫星，先固定宽巷双差模糊度，然后将固定的宽巷双差模糊度代入到无电离层组合，并利用卡尔曼滤波估计L1载波相位的整周模糊度和相对对流层天顶延迟，得到双差模糊度浮点解和方差阵，最后利用改进的LAMBDA方法实时搜索双差整周模糊度；
按照基线从短到长的顺序分别固定选择基线；
根据固定结果进行法方程更新，以得到模糊度参数实数解；
根据模糊度固定规则将模糊度参数实数解解固定模糊度整数解；
其中，所述根据固定结果进行法方程更新，以得到模糊度参数实数解包括以下步骤：
求解出所述双差模糊度的实数解；
将所述实数解中的模糊度参数固定为整数，并进行检验，得到整周模糊度；
将所述整周模糊度作为已知值代入法方程进行更新；
按上述过程进行迭代，直至得到所有整周模糊度参数；
所述预设选星规则为：
GPS和BDS每个系统各选一颗卫星作为参考星，在BDS系统中，首选GEO卫星作为参考星，若由于环境遮挡任何一颗GEO卫星均不为参考星，则依次选择IGSO、MEO卫星为参考星；在GPS卫星中，选择高度角较高的卫星作为参考星。
5.如权利要求1所述的基于虚拟基准站的煤矿采空区地表形变监测方法，其特征在于，所述步骤S05中对流层延迟量的计算方法为：
式中，为双差对流层延迟干分量，通过经验模型计算得到，/>为天顶对流层湿延迟，/>为测站m和卫星i的高度角，/>为映射函数，/>表示双差对流层湿延迟分量；
所述步骤S05中电离层延迟量的计算方法为：
其中，和/>分别为L1频段和L2频段观测值的频率，/>和/>分别为L1频段和L2频段观测值波长，/>和/>分别为L1频段和L2频段双差载波观测值，/>和/>分别为L1频段和L2频段双差整周模糊度；
所述步骤S05中综合误差的计算方法为：
其中，表示卫星与接收机之间的几何距离，/>表示载波观测值，/>表示整周模糊度，/>表示载波的波长，/>表示电离层延迟量，/>表示对流层延迟量。
6.如权利要求1所述的基于虚拟基准站的煤矿采空区地表形变监测方法，其特征在于，所述步骤S07具体包括：
对内插区域内监测站的电离层延迟量的计算方法为：
在中心电离层高度处进行空间结构性分析和空间变异性分析；
将预设高度作为电离层误差改正基准面；
在所述电离层误差改正基准面进行监测站非差电离层延迟误差内插计算，得到监测站处的电离层改正数，根据监测站处的电离层改正数计算得到电离层延迟量；
所述非差电离层延迟误差内插计算包括：
在线性无偏、最优估计条件下，得到电离层总电子含量插值的克里金方程组，所述克里金方程组为：
式中，μ表示拉格朗日乘数因子，dij表示基准站与监测站间的距离，di0表示插值点与观测值间的距离，λj表示加权系数；
对克里金方程组进行求解，得到加权系数；
基于所述加权系数计算得到电离层总电子含量估计值；
通过所述电离层总电子含量估计值计算得到电离层改正数；
对内插区域内监测站的对流层延迟量的计算方法为：
将对流层相关误差分为干延迟部分和湿延迟部分，干延迟部分通过GPT2w或ITG模型进行改正，湿延迟部分根据各基准站解算的对流层湿延迟结果采用基于高度改正的Kriging内插算法，计算出监测站处的对流层改正数，根据监测站处的对流层改正数计算得到对流层延迟量；
所述基于高度改正的Kriging内插算法具体包括：
将各基准站的高程取平均值得到内插区域的对流层湿延迟改正的高程基准面；
根据对流层湿延迟高差改正模型计算出，各基准站到对流层误差改正基准面的对流层湿延迟改正数；
在高程基准面上对各个基准站对流层湿延迟进行空间结构性分析和变异性分析；
监测站根据监测站概略坐标，根据对流层湿延迟的空间结构性分析和变异性分析结果，采用Kriging算法内插出该处对流层湿延迟改正基准面的对流层湿延迟改正数；
根据监测站概略坐标高程，计算出监测站到对流层误差改正基准面的高差；
根据对流层湿延迟高差改正模型计算出监测站概略坐标高程处的对流层湿延迟改正数；
对内插区域内监测站的综合误差的计算方法为：
通过反距离权重法进行内插得到综合误差。
7.如权利要求1所述的基于虚拟基准站的煤矿采空区地表形变监测方法，其特征在于，所述步骤S08具体包括：
服务器基于监测站发来的概略坐标，构建虚拟基准站，通过内插算法获得虚拟基准站相对于CORS站的电离层、对流层和几何距离的相对差值，利用虚拟基准站和CORS站的坐标和卫星星历，将CORS站的观测值归算至虚拟基准站，然后再顾及对流层延迟误差和电离层延迟误差的站间相对差值，获得虚拟基准站的观测值；
计算CORS站的对流层延迟误差双差后的差值，计算方法如下：
计算CORS站处的电离层延迟误差双差后的差值，计算方法如下：
其中，和/>分别为L1频段和L2频段观测值的频率，/>和/>分别为L1频段和L2频段双差载波观测值，/>和/>分别为L1频段和L2频段双差整周模糊度，/>表示卫星与接收机之间的几何距离；
通过空间插值算法计算出虚拟基准站处的对流层延迟误差双差后的差值和电离层延迟误差双差后的差值/>；
计算虚拟基准站的观测值。
8.如权利要求7所述的基于虚拟基准站的煤矿采空区地表形变监测方法，其特征在于，所述虚拟基准站观测值包括虚拟基准站载波观测值和虚拟基准站伪距观测值，所述虚拟基准站载波观测值的计算方法如下：
其中，表示虚拟基准站载波观测值双差误差改正数，/>包括频点1和频点2分别对应的载波观测值双差误差改正数/> ，计算方法如下：
所述虚拟基准站伪距观测值的计算方法如下：
其中，表示虚拟基准站伪距双差误差改正数，/>包括频点1和频点2分别对应的虚拟基准站伪距双差误差改正数/>和/>，计算方法如下：
其中，下标v表示虚拟基准站。
9.如权利要求1所述的基于虚拟基准站的煤矿采空区地表形变监测方法，其特征在于，所述步骤S10具体包括：
预测状态向量及方差-协方差阵，其中，预测状态向量：
预测方差-协方差阵：
其中，表示/>时刻的坐标向量的估计值，/>表示/>时刻方差-协方差阵的估计值，/>表示/>时刻坐标向量预测值，/>表示/>时刻方差-协方差阵的估计值，/>为/>时刻至/>时刻的系统状态转移矩阵，/>是/>时刻系统过程噪声的方差矩阵；
更新状态向量及方差-协方差阵，包括：
计算滤波增益：
更新状态向量：
更新方差-协方差阵：
其中，表示/>时刻观测方程的系数矩阵，/>表示/>时刻观测噪声的方差矩阵，/>表示/>时刻的坐标观测值向量，I为单位矩阵，初始条件为/>和/>；
将卡尔曼滤波每一步迭代得到的坐标状态向量组成监测站坐标时间序列。
10.如权利要求1所述的基于虚拟基准站的煤矿采空区地表形变监测方法，其特征在于，计算两相邻点在竖直方向的相对移动与两相邻点间水平距离的比值确定倾斜变形；计算两相邻线段的倾斜差与两线段中点间的水平距离的比值确定曲率变形；计算两相邻点的水平移动差值与两点间水平距离的比值确定水平变形。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311479876.4A CN117233799B (zh) | 2023-11-08 | 2023-11-08 | 基于虚拟基准站的煤矿采空区地表形变监测方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202311479876.4A CN117233799B (zh) | 2023-11-08 | 2023-11-08 | 基于虚拟基准站的煤矿采空区地表形变监测方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN117233799A true CN117233799A (zh) | 2023-12-15 |
| CN117233799B CN117233799B (zh) | 2024-02-09 |
Family
ID=89089618
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202311479876.4A Active CN117233799B (zh) | 2023-11-08 | 2023-11-08 | 基于虚拟基准站的煤矿采空区地表形变监测方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN117233799B (zh) |
Citations (15)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2011100680A2 (en) * | 2010-02-14 | 2011-08-18 | Trimble Navigation Limited | Gnss signal processing with regional augmentation positioning |
| EP3130943A1 (en) * | 2015-08-14 | 2017-02-15 | Trimble Inc. | Navigation satellite system positioning involving the generation of tropospheric correction information |
| CN106646538A (zh) * | 2016-10-31 | 2017-05-10 | 东南大学 | 一种基于单差滤波的变形监测gnss信号多路径改正方法 |
| CN107797126A (zh) * | 2017-09-26 | 2018-03-13 | 东南大学 | 基于星型网络的bds/gps广播式网络rtk算法 |
| US20180210089A1 (en) * | 2015-07-21 | 2018-07-26 | Comnav Technology Ltd. | Un-differential correction distributed processing system and method based on receiver of reference station |
| WO2019174113A1 (zh) * | 2018-03-16 | 2019-09-19 | 东南大学 | 一种gps/bds紧组合载波差分定位方法 |
| CN112902825A (zh) * | 2021-04-13 | 2021-06-04 | 长安大学 | 一种适用于高精度形变监测的北斗/gnss网络rtk算法 |
| WO2021237804A1 (zh) * | 2020-05-29 | 2021-12-02 | 湖南联智科技股份有限公司 | 基于北斗高精度定位的基础设施结构变形监测方法 |
| CN114152185A (zh) * | 2022-02-07 | 2022-03-08 | 山东科技大学 | 一种gnss形变监测系统及其工作方法 |
| CN114859390A (zh) * | 2022-04-21 | 2022-08-05 | 武汉大学 | 一种高精度cors电离层改正的ftk解算方法 |
| CN115421172A (zh) * | 2022-11-04 | 2022-12-02 | 南京市计量监督检测院 | 一种基于实时与准实时结合的北斗变形监测方法 |
| CN115963522A (zh) * | 2022-11-29 | 2023-04-14 | 国网思极位置服务有限公司 | 一种结合基准站卫星数据的定位方法与终端 |
| CN116026226A (zh) * | 2022-10-17 | 2023-04-28 | 南京凌远时空科技有限公司 | 一种半遮挡环境下的水闸变形监测方法及系统 |
| CN116125514A (zh) * | 2023-02-09 | 2023-05-16 | 广州市城市规划勘测设计研究院 | 基于北斗ppp-rtk虚拟观测值地灾监测方法、装置、终端及介质 |
| CN116520378A (zh) * | 2023-07-03 | 2023-08-01 | 武汉大学 | 非差rtk误差改正数确定方法、装置、设备及存储介质 |
-
2023
- 2023-11-08 CN CN202311479876.4A patent/CN117233799B/zh active Active
Patent Citations (15)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| WO2011100680A2 (en) * | 2010-02-14 | 2011-08-18 | Trimble Navigation Limited | Gnss signal processing with regional augmentation positioning |
| US20180210089A1 (en) * | 2015-07-21 | 2018-07-26 | Comnav Technology Ltd. | Un-differential correction distributed processing system and method based on receiver of reference station |
| EP3130943A1 (en) * | 2015-08-14 | 2017-02-15 | Trimble Inc. | Navigation satellite system positioning involving the generation of tropospheric correction information |
| CN106646538A (zh) * | 2016-10-31 | 2017-05-10 | 东南大学 | 一种基于单差滤波的变形监测gnss信号多路径改正方法 |
| CN107797126A (zh) * | 2017-09-26 | 2018-03-13 | 东南大学 | 基于星型网络的bds/gps广播式网络rtk算法 |
| WO2019174113A1 (zh) * | 2018-03-16 | 2019-09-19 | 东南大学 | 一种gps/bds紧组合载波差分定位方法 |
| WO2021237804A1 (zh) * | 2020-05-29 | 2021-12-02 | 湖南联智科技股份有限公司 | 基于北斗高精度定位的基础设施结构变形监测方法 |
| CN112902825A (zh) * | 2021-04-13 | 2021-06-04 | 长安大学 | 一种适用于高精度形变监测的北斗/gnss网络rtk算法 |
| CN114152185A (zh) * | 2022-02-07 | 2022-03-08 | 山东科技大学 | 一种gnss形变监测系统及其工作方法 |
| CN114859390A (zh) * | 2022-04-21 | 2022-08-05 | 武汉大学 | 一种高精度cors电离层改正的ftk解算方法 |
| CN116026226A (zh) * | 2022-10-17 | 2023-04-28 | 南京凌远时空科技有限公司 | 一种半遮挡环境下的水闸变形监测方法及系统 |
| CN115421172A (zh) * | 2022-11-04 | 2022-12-02 | 南京市计量监督检测院 | 一种基于实时与准实时结合的北斗变形监测方法 |
| CN115963522A (zh) * | 2022-11-29 | 2023-04-14 | 国网思极位置服务有限公司 | 一种结合基准站卫星数据的定位方法与终端 |
| CN116125514A (zh) * | 2023-02-09 | 2023-05-16 | 广州市城市规划勘测设计研究院 | 基于北斗ppp-rtk虚拟观测值地灾监测方法、装置、终端及介质 |
| CN116520378A (zh) * | 2023-07-03 | 2023-08-01 | 武汉大学 | 非差rtk误差改正数确定方法、装置、设备及存储介质 |
Non-Patent Citations (3)
| Title |
|---|
| YANG WANG: "A Refinement Method of Real-Time Ionospheric Model for China", REMOTE SENSING, pages 1 - 19 * |
| 吴北平, 李征航: "GPS网络RTK线性组合法与内插法关系的讨论", 测绘信息与工程, no. 05, pages 27 - 28 * |
| 梁霄;杨玲;黄涛;王延兵;: "网络RTK基准站间的模糊度及空间相关误差解算", 测绘工程, no. 01, pages 24 - 28 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN117233799B (zh) | 2024-02-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN108205150B (zh) | 差分定位方法及系统 | |
| Paziewski et al. | Assessment of GPS+ Galileo and multi-frequency Galileo single-epoch precise positioning with network corrections | |
| CN114518586B (zh) | 一种基于球谐展开的gnss精密单点定位方法 | |
| EP2995975B1 (en) | Precise gnss positioning system with improved ambiguity estimation | |
| Li et al. | The GFZ real-time GNSS precise positioning service system and its adaption for COMPASS | |
| CN111983654B (zh) | 一种基于gnss的北极区域电离层相位闪烁因子构建方法 | |
| CN111290005B (zh) | 载波相位的差分定位方法、装置、电子设备及存储介质 | |
| US20090135057A1 (en) | Real-time fast decimeter-level GNSS positioning | |
| CN108549095B (zh) | 一种区域cors网非差并行增强方法及系统 | |
| Paziewski | Study on desirable ionospheric corrections accuracy for network-RTK positioning and its impact on time-to-fix and probability of successful single-epoch ambiguity resolution | |
| EP3163324B1 (en) | Positioning device, positioning method, and program | |
| Wielgosz et al. | Troposphere modeling for precise GPS rapid static positioning in mountainous areas | |
| CN111694030A (zh) | 一种基于格网虚拟观测值的bds局域差分方法及系统 | |
| CN107966722B (zh) | 一种gnss钟差解算方法 | |
| CN112462396B (zh) | 一种高采样率导航卫星钟差的实时并行确定方法 | |
| Zus et al. | A forward operator and its adjoint for GPS slant total delays | |
| CN111381264A (zh) | 网络rtk中长基线模糊度固定方法和平台 | |
| JP5759676B2 (ja) | 伝搬経路推定システム及び伝搬経路推定方法 | |
| Chen et al. | A geometry-free and ionosphere-free multipath mitigation method for BDS three-frequency ambiguity resolution | |
| Zhao et al. | Accuracy and reliability of tropospheric wet refractivity tomography with GPS, BDS, and GLONASS observations | |
| Landau et al. | Trimble’s RTK and DGPS solutions in comparison with precise point positioning | |
| CN112731490A (zh) | 一种rtk定位方法及装置 | |
| Wang et al. | Effect of biases in integrity monitoring for RTK positioning | |
| CN114859390A (zh) | 一种高精度cors电离层改正的ftk解算方法 | |
| Brack et al. | Operational multi-GNSS global ionosphere maps at GFZ derived from uncombined code and phase observations |
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 |