CN108897068B - 气象台站面密度度量方法 - Google Patents

气象台站面密度度量方法 Download PDF

Info

Publication number
CN108897068B
CN108897068B CN201810364309.7A CN201810364309A CN108897068B CN 108897068 B CN108897068 B CN 108897068B CN 201810364309 A CN201810364309 A CN 201810364309A CN 108897068 B CN108897068 B CN 108897068B
Authority
CN
China
Prior art keywords
station
measurement method
areal density
density measurement
density
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
CN201810364309.7A
Other languages
English (en)
Other versions
CN108897068A (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.)
Nanjing University
Original Assignee
Nanjing University
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 Nanjing University filed Critical Nanjing University
Priority to CN201810364309.7A priority Critical patent/CN108897068B/zh
Publication of CN108897068A publication Critical patent/CN108897068A/zh
Application granted granted Critical
Publication of CN108897068B publication Critical patent/CN108897068B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01W—METEOROLOGY
    • G01W1/00—Meteorology
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Environmental & Geological Engineering (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Environmental Sciences (AREA)
  • Algebra (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Ecology (AREA)
  • Evolutionary Biology (AREA)
  • Biodiversity & Conservation Biology (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Atmospheric Sciences (AREA)
  • Length Measuring Devices With Unspecified Measuring Means (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了气象台站面密度度量方法，该方法利用核函数近似对台站面密度进行平滑处理，并引入局地尺度等，最终获得自适应的气象台站面密度；该方法可以消除现有计算台站密度方法的主观性、不连续性及单一尺度性，可以满足台站密度的客观度量及进一步使用插值、EOF分析等统计方法的要求。

Description

气象台站面密度度量方法
技术领域
本发明涉及大气科学研究领域的气象台站面密度的度量，尤其涉及气象台站面密度度量的方法。
背景技术
气象台站的面密度是一个反映气象台站分布特点的基本量，有着广泛地应用，如许多插值方法需要用到台站的面密度(Stahl et al.2006)，如观测站网的代表性分析(Renand Ren 2012)，如台站资料的经验正交函数(EOF)分析(王盘兴等2011)等。现有研究中的台站面密度主要是用原始的面密度定义求得的，即用给定区域的台站数除以区域面积：
Figure GSB0000188459250000011
其中子区域
Figure GSB0000188459250000012
Ω为待研究的总区域。S(ω)表示子区域ω的面积，n(ω)表示子区域ω内的台站数。(注：本文所用的面密度一词是相对于线密度或体密度而言的，表示二维平面中的密度。)
实际使用时常选取某特定半径d的圆形区域进行统计，此时的子区域ω记为
Figure GSB0000188459250000013
且有
Figure GSB0000188459250000014
其中A为任意变量场，即可得
Figure GSB0000188459250000015
S(ω)＝πd2，故尺度d下的面密度为
Figure GSB0000188459250000016
其中
Figure GSB0000188459250000017
为总区域内的坐标，子区域ω取为以
Figure GSB0000188459250000018
为中心、半径为d的圆形区域。
然而，这种直接的定义有几点不足：1.需要人为划分区域或给定统计半径d，因而结果具有主观性，不利于客观化分析。2.面密度会呈现不合理的不连续跃变，影响面密度分布的判断及后续的分析。3.通常只能得到一个统一尺度d下的面密度分布特征，而实际的台站分布疏密差异极大，呈现多尺度特征，这会导致基于散点数据的统计方法的不准确，如散点数据的EOF分析(王盘兴等2011)。
发明内容
本发明所要解决的技术问题是针对现有技术中的不足，提供一种新的计算气象台站面密度的方法，该方法可以消除现有计算台站密度方法的主观性、不连续性及单一尺度性，可以满足台站面密度的客观度量及进一步使用插值、EOF分析等统计方法的要求。
为了解决上述技术问题，本发明气象台站面密度度量方法，包含对台站面密度进行平滑处理；具体处理方法是，首先，子区域ω内的台站数n(ω)用狄拉克函数表示为
Figure GSB0000188459250000021
其中Ω为待研究的总区域，N总的台站数，
Figure GSB0000188459250000022
是以
Figure GSB0000188459250000023
为圆心、半径为d的圆形区域，
Figure GSB0000188459250000024
为区域ω内的坐标，
Figure GSB0000188459250000025
为台站j的坐标，ds为由
Figure GSB0000188459250000026
决定的面积微元，
Figure GSB0000188459250000027
为狄拉克函数，满足
Figure GSB0000188459250000028
且
Figure GSB0000188459250000029
其次，将狄拉克函数
Figure GSB00001884592500000210
进一步用核函数
Figure GSB00001884592500000211
来近似，即
Figure GSB00001884592500000212
使用核函数，获得光滑台站面密度
Figure GSB00001884592500000213
所述
Figure GSB00001884592500000214
其中〈>表示对变量的平滑近似。
为了能够使计算更简化，所述
Figure GSB00001884592500000215
其中权重
Figure GSB00001884592500000216
为了在区域边界附近及台站密度大于网格的区域时仍有足够的精度，所述
Figure GSB00001884592500000217
其中权重
Figure GSB00001884592500000218
为积分变量。
以上得到的台站面密度解决了原始方法不连续的问题，但仍然需要事先选定一个统一的尺度d，下面进一步得到自适应的台站面密度，本文中也称为局地光滑台站面密度。(注1：本文中自适应指无需人为选定适合的一个或多个尺度d，而能够自动得到适宜于台站分布特点的随空间变化的尺度及台站面密度。注2：本文所用的局地一词表示随着空间点的变化而变化，即相对于所有空间点都使用一个全局的d而言的。)
上述技术方案中，使d随着位置而变化，也就是将d替换为局地尺度
Figure GSB00001884592500000219
所述
Figure GSB0000188459250000031
其中
Figure GSB0000188459250000032
满足关系
Figure GSB0000188459250000033
其中M为正整数，
Figure GSB0000188459250000034
为以
Figure GSB0000188459250000035
为中心、包含M个台站的圆的半径，κ≥1为平滑参数，κ取值越大获得的结果越平滑，j＝1，...，N为各台站标号，min()为取最小值的函数，count()为计数函数，即满足条件的参数的个数；整个式子表示的意义是，局地尺度为
Figure GSB0000188459250000036
包含M个台站的最小半径的κ倍。
进而获得局地光滑台站面密度
Figure GSB0000188459250000037
所述
Figure GSB0000188459250000038
其中〈>local表示核函数近似及引入局地尺度后的变量场。
为了减少计算量，所述
Figure GSB0000188459250000039
其中权重
Figure GSB00001884592500000310
为了在区域边界附近及台站密度大于网格的区域仍有足够的精度，
所述
Figure GSB00001884592500000311
其中权重
Figure GSB00001884592500000312
为积分变量。
上述技术方案中，所述局地尺度
Figure GSB00001884592500000313
中M＝2，κ＝1.3。核函数取为三次样条核函数(Monaghan&Lattanzio 1985)：
Figure GSB00001884592500000314
其中
Figure GSB00001884592500000315
本发明的技术方案，可实现以下的有益效果：
1.该方法能解决传统方法得到的台站面密度分布不连续的缺陷，得到物理上更为合理的连续分布的台站面密度场。以满足定性及定量应用的要求。
2.该方法能摆脱选取尺度的主观性，得到自适应的台站面密度分布。
3.该方法能解决传统方法只能得到单一尺度的台站面密度分布的缺陷，对于台站疏密差异很大的情形，能得到符合台站多尺度特点的面密度分布。
4.该方法能进一步应用于后续的气象统计，以提高后续统计的效果。如散点资料的EOF分析。
附图说明
图1为理想的台站分布，网格间距为0.5，黑点为台站；
图2为图1相应的局地面密度；
图3为IGRA全球探空台站分布图；
图4为IGRA全球探空台站局地面密度；
图5为计算局地光滑台站面密度的流程图。
具体实施方式
设Ω为待研究的总区域。总的台站数为N。区域
Figure GSB0000188459250000041
为子区域其面积为S(ω)，区域内有n(ω)个台站，则单位面积的台站数即台站面密度ρ(ω)的原始定义为：
Figure GSB0000188459250000042
实际使用时常选取某特定半径d的圆形区域进行统计，此时的子区域ω记为
Figure GSB0000188459250000043
且有
Figure GSB0000188459250000044
S(ω)＝πd2，故尺度d下的面密度为
Figure GSB0000188459250000045
其中
Figure GSB0000188459250000046
为总区域内的坐标，子区域ω取为以
Figure GSB0000188459250000047
为中心、半径为d的圆形区域。
为了在区域边界附近及台站密度大于网格的区域仍有足够的精度，，(2)式应写为
Figure GSB0000188459250000048
使用该式可直接衡量台站的疏密程度，下文中称该方法为原方法。
下面分成两步对原方法进行改进。
1.光滑台站面密度
将区域ω内的台站数n(ω)用狄拉克函数表示为
Figure GSB0000188459250000051
其中
Figure GSB0000188459250000052
为狄拉克函数，满足
Figure GSB0000188459250000053
且
Figure GSB0000188459250000054
ds为由
Figure GSB0000188459250000055
决定的面积微元。而狄拉克函数
Figure GSB0000188459250000056
可进一步用核函数
Figure GSB0000188459250000057
来近似，即
Figure GSB0000188459250000058
该过程称为核函数近似(Monaghan 2005)，其中核函数
Figure GSB0000188459250000059
需满足(1)归一性
Figure GSB00001884592500000510
(2)趋向于狄拉克函数
Figure GSB00001884592500000511
(3)正定性
Figure GSB00001884592500000512
(4)紧性：
Figure GSB00001884592500000513
c为某一常数。
使用核函数近似，获得光滑台站面密度
Figure GSB00001884592500000514
Figure GSB00001884592500000515
其中〈>表示对变量的平滑近似，关于核函数的描述详见于光滑粒子流体动力学((Monaghan 2005；Liu and Liu 2010))。使用核函数，能解决原方法求得的台站面密度不连续的缺陷。
为提高计算效率，也可采用下式计算
Figure GSB00001884592500000516
其中权重
Figure GSB00001884592500000517
为了在区域边界附近及台站密度大于网格的区域仍有足够的精度，(5)式应写为
Figure GSB00001884592500000518
其中权重
Figure GSB00001884592500000519
为积分变量。
至此，能够解决原方法的密度不连续跃变问题。
2.局地光滑台站面密度
上述求得的台站密度必须事先给定一个固定的尺度，这是主观化的，且会导致其余尺度信息的丢失：若尺度取得过大，无法分辨出较小尺度的台站密度变化；若尺度取得过小，又会丢失大尺度的台站分布信息。为解决这一问题，下面定义一种局地光滑台站面密度。
首先，令d随着位置而变化，定义局地尺度
Figure GSB0000188459250000061
Figure GSB0000188459250000062
其中
Figure GSB0000188459250000063
满足关系
Figure GSB0000188459250000064
其中M可取为正整数，
Figure GSB0000188459250000065
为以
Figure GSB0000188459250000066
为中心、包含M个台站的圆的半径，κ≥1为平滑参数，越大得到的结果越平滑，j＝1，...，N为各台站标号，min()为取最小值的函数，count()为计数函数，即满足条件的参数的个数；整个式子表示的意义是，局地尺度为
Figure GSB0000188459250000067
包含M个台站的最小半径的κ倍。现在，可进一步定义局地光滑台站面密度
Figure GSB0000188459250000068
其中〈>local表示核函数近似及引入局地尺度后的变量场，然而该式的计算量太大，故实际使用时，采用以下近似式：
Figure GSB0000188459250000069
其中权重
Figure GSB00001884592500000610
若考虑区域边界及网格离散化的影响，(9)式应写为
Figure GSB00001884592500000611
其中权重
Figure GSB00001884592500000612
为积分变量。
这样就得到了包含台站多尺度特点的一种自适应的光滑面密度分布，即进一步解决了尺度选取的主观性问题及不同尺度信息的丢失问题。
实际使用时，上述技术方案中，按以下方式选取参数及核函数效果较好：
所述局地尺度
Figure GSB0000188459250000071
中M＝2，κ＝1.3。核函数取为三次样条核函数(Monaghan&Lattanzio 1985)：
Figure GSB0000188459250000072
其中
Figure GSB0000188459250000073
以图1理想的台站分布及图3的IGRA全球探空台站分布为例来说明本发明具体的实现。
先计算局地尺度。
所述公式(7)中，取M＝2，κ＝1.3，即
Figure GSB0000188459250000074
其中
Figure GSB0000188459250000075
满足关系
Figure GSB0000188459250000076
本实施例中，子区域Ωi取为以
Figure GSB0000188459250000077
为圆心，半径为d的圆形区域。对于图1的理想台站分布，网格划分为0.5×0.5，使用常规的欧氏距离进行计算。对于图3的IGRA(IntegratedGlobal Radiosonde Archive)全球探空台站集，划分为0.5°×0.5°的经纬度网格，使用Haversine球面距离：
Figure GSB0000188459250000078
其中r＝6378km为地球半径，
Figure GSB0000188459250000079
及λ1、λ2分别为
Figure GSB00001884592500000710
的纬度及经度。球坐标中的面积微元
Figure GSB00001884592500000711
其中r为地球半径，λ为经度，
Figure GSB00001884592500000712
为纬度。
而后，选取三次样条核函数(11)式为核函数，可用(10)式求得自适应的台站面密度分别如图2、4所示。
其计算流程见图5。
下面是对局地面密度计算结果的说明。
图2为图1中理想台站分布的局地面密度，该局地面密度含有不同尺度的信息。上方从左至右，显示出了台站由疏变密的分布，这种由疏变密不是由于台站间距变化引起的，而是由于台站从一维的排布渐变到二维的排布而引起的。左下方区域体现出了整个区域的低密度及中心区域的高密度。右下方区域从左至右反应出了密度逐渐减小的特征，同时保留了图1右侧几个密度最大的中心的信息。可见，局地面密度能够解决台站分布的疏密差异巨大时密度的度量问题。
图4为图3的IGRA全球探空台站集的局地面密度分布。计算结果表明，亚洲东南部、马来群岛、西欧的台站面密度最大，基本超过10个/1000km2，其次是北美洲、非洲东南部等区域，北冰洋的台站面密度大于南极洲，台站最稀疏的区域是南太平洋、南印度洋和南大西洋，小于0.25个/1000km2。可见，全球IGRA探空台站面密度的差异量级达102以上。

Claims (9)

1.气象台站面密度度量方法，其特征在于：对台站面密度进行平滑处理；具体处理方法是，首先，子区域ω内的台站数n(ω)用狄拉克函数表示为
Figure FSB0000189239090000011
其中Ω为待研究的总区域，N为总的台站数，
Figure FSB0000189239090000012
是以
Figure FSB0000189239090000013
为圆心、半径为d的圆形区域，
Figure FSB0000189239090000014
为区域ω内的坐标，
Figure FSB0000189239090000015
为台站j的坐标，ds为由
Figure FSB0000189239090000016
决定的面积微元，
Figure FSB0000189239090000017
为狄拉克函数，满足
Figure FSB0000189239090000018
且
Figure FSB0000189239090000019
其次，将狄拉克函数
Figure FSB00001892390900000110
进一步用核函数
Figure FSB00001892390900000111
来近似，即
Figure FSB00001892390900000112
使用核函数，获得光滑台站面密度
Figure FSB00001892390900000113
所述
Figure FSB00001892390900000114
其中<>表示对变量的平滑近似。
2.如权利要求1的气象台站面密度度量方法，其特征在于：所述
Figure FSB00001892390900000115
其中权重
Figure FSB00001892390900000116
3.如权利要求1的气象台站面密度度量方法，其特征在于：所述
Figure FSB00001892390900000117
其中权重
Figure FSB00001892390900000118
为积分变量。
4.如权利要求1的气象台站面密度度量方法，其特征在于：将d替换为局地尺度
Figure FSB00001892390900000119
所述
Figure FSB00001892390900000120
其中
Figure FSB00001892390900000121
满足关系
Figure FSB00001892390900000122
M为正整数，
Figure FSB00001892390900000123
为以
Figure FSB00001892390900000124
为中心、包含M个台站的圆的半径，κ≥1为平滑参数，j＝1，...，N为各台站标号，min()为取最小值的函数，count()为计数函数，进而获得局地光滑台站面密度
Figure FSB00001892390900000125
所述
Figure FSB0000189239090000021
其中<>local表示核函数近似及引入局地尺度后的变量场。
5.如权利要求4的气象台站面密度度量方法，其特征在于：所述
Figure FSB0000189239090000022
其中权重
Figure FSB0000189239090000023
6.如权利要求4的气象台站面密度度量方法，其特征在于：所述
Figure FSB0000189239090000024
其中权重
Figure FSB0000189239090000025
为积分变量。
7.如权利要求4的气象台站面密度度量方法，其特征在于：所述局地尺度
Figure FSB0000189239090000026
κ＝1.3。
8.如权利要求1至6中任意一项的气象台站面密度度量方法，其特征在于：所述核函数取为
Figure FSB0000189239090000027
其中
Figure FSB0000189239090000028
9.如权利要求7的气象台站面密度度量方法，其特征在于：所述核函数取为
Figure FSB0000189239090000029
其中
Figure FSB00001892390900000210
CN201810364309.7A 2018-04-20 2018-04-20 气象台站面密度度量方法 Active CN108897068B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810364309.7A CN108897068B (zh) 2018-04-20 2018-04-20 气象台站面密度度量方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810364309.7A CN108897068B (zh) 2018-04-20 2018-04-20 气象台站面密度度量方法

Publications (2)

Publication Number Publication Date
CN108897068A CN108897068A (zh) 2018-11-27
CN108897068B true CN108897068B (zh) 2020-11-13

Family

ID=64342407

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810364309.7A Active CN108897068B (zh) 2018-04-20 2018-04-20 气象台站面密度度量方法

Country Status (1)

Country Link
CN (1) CN108897068B (zh)

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004501358A (ja) * 2000-05-11 2004-01-15 ベクトン・ディキンソン・アンド・カンパニー 最適な境界を有する平滑化された多角形を使用して散布図中のクラスタを識別するシステム
JP5670832B2 (ja) * 2011-05-24 2015-02-18 富士通株式会社 シミュレーション方法、シミュレーション装置及びシミュレーションプログラム
CN102254090B (zh) * 2011-06-08 2013-07-03 南京信息工程大学 采用核密度估计的闪电密度分布及ng值估计方法

Also Published As

Publication number Publication date
CN108897068A (zh) 2018-11-27

Similar Documents

Publication Publication Date Title
Zhao et al. A method for simplifying ship trajectory based on improved Douglas–Peucker algorithm
CN108491757B (zh) 基于多尺度特征学习的光学遥感图像目标检测方法
WO2022016884A1 (zh) 一种基于K-means聚类算法的海面风速方法
CN107730527B (zh) 一种基于遥感卫星影像的高原地区冰湖提取方法
CN101599120B (zh) 一种遥感影像建筑物识别方法
CN107016677A (zh) 一种基于fcn和cnn的云图分割方法
CN106778749B (zh) 基于聚集度和Delaunay三角重构的巡回作业区域边界提取方法
CN102938161B (zh) 一种基于Mean Shift的三维形状自动分割方法
CN110596008B (zh) 基于地块的中国洪积平原农业区土壤养分数字制图方法
CN107301649B (zh) 一种基于超像素的区域合并sar图像海岸线检测算法
WO2018168165A1 (ja) 気象予測装置、気象予測方法、およびプログラム
CN103606154A (zh) 基于jseg和谱聚类的多尺度海面溢油sar图像分割方法
CN114067118B (zh) 一种航空摄影测量数据的处理方法
Zou et al. A method of radar echo extrapolation based on TREC and Barnes filter
Ren et al. Detection of SST fronts from a high-resolution model and its preliminary results in the south China sea
CN117541930A (zh) 流域尺度高时空分辨率河流网络遥感信息提取方法及系统
CN107274361B (zh) Landsat TM遥感影像数据除云方法及系统
CN115223000A (zh) 面向耕地资源监测的深度学习样本制作方法
CN108897068B (zh) 气象台站面密度度量方法
CN113011520B (zh) 基于热力图与图论聚类的雷达数据过滤方法
CN109191484B (zh) 一种从机载激光雷达点云中快速提取平面片的方法
Cabral et al. Space–time trends and dependence of precipitation extremes in North‐Western Germany
CN107632328A (zh) 一种近海海洋气象信息分析处理方法
CN117171855A (zh) 一种基于Delaunay三角剖分的丘陵区流场模型建模方法
CN110175638B (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