本发明涉及移动通信技术领域，具体涉及一种利用移动通信大数据对一人多号用户进行识别的方法。
背景技术：
鉴于统计人口信息的重要性，目前已衍生出很多种间接估计人口信息的方法。比如地区生活必需品消费量的变化来估算人口数量，如用电量法、用水量法、食盐销售量法等。这些方法的缺陷在于获取的数据精度不高，且有一定滞后性，成本也比较高，无法快速准确获取区域人口数据信息。另外基本的人口普查、人口抽查以及少数省市不定时的专题性人口调查，也是获取人口信息的方法。但总体来说，缺乏针对人口信息的快速动态监测的有效方法与方案，这已成为各级政府在进行科学决策和采取应急行动时的主要制约因素之一。同时伴随着我国人口流动日益频繁，城市人口分布及结构正面临重大变迁，而及时掌握人口信息的手段却仍然十分有限。
鉴于当前城市人口调控需求，以及随着移动通信规模的不断扩大以及技术的不断发展，大规模通信数据特别是轨迹数据的存储已实现，因此以移动通信大数据分析为研究手段对城市人口规模及流动情况进行估算已成为可能。统计部门需要运用大数据开展人口数据的动态监测，完善现有的人口监测体系，探索产业调控、功能布局与人口发展之间的关系，并实现定期追踪掌握疏解人口流向，及时预警区域人口变化趋势。相比于传统的研究方法，基于大数据统计及人口监测将具有更高的可信度和准确性。
目前电信运营商在面对日益激烈的市场竞争压力时，不断降低移动通信资费并不断推出新的促销方式，以推动移动用户数量的增长，因此一人多卡甚至一人多机、多运营商现象日益增多，在基于通信大数据进行人口统计和监测的过程中，对一人双(多)卡用户的识别是目前遇到的一类突出问题。而现有的方法和系统中没有有效的手段对多卡用户进行判断识别，会造成一定的统计误差。对这一问题的有效解决，可以大大增加通信大数据进行人口统计的准确性。
技术实现要素：
本发明的目的在于提供一种基于距离的一人多号(卡/机)用户识别方法，解决通信数据中的一人多卡情况影响大数据人口统计分析准确性的问题。
为实现上述目的，本发明采用以下技术方案：
一种基于距离的一人多号移动用户识别方法，其特征在于：
(1)数据采集：利用电信运营商信令数据，对每个imsi识别号所在的基站位置、进出基站时间点信息进行采集并存储；
(2)数据筛选：每次读取一个imsi对，设为aimsi和bimsi，作为待确定的一人多号对；在一定时间段内，按时间点先后排序，分别得到这两个imsi的各个时间点以及在各个时间点所对应的基站位置的数据；
(3)距离计算：根据每个具体基站所对应的经纬度，计算两个imsi在各时间点下的距离，计算方法如下：
(31)先找到aimsi的一个时间点t，再找到bimsi在t时间点前的最后一个时间点t1以及在t时间点后的第一个时间点t2；
(32)对应出aimsi在时间点t所在位置的经纬度，以及bimsi在时间点t1和t2时的位置经纬度；
(33)分别计算出aimsi在t时间点和bimsi在t1时间点时的距离，以及aimsi在t时间点和bimsi在t2时间点时的距离，计算公式如下：
dis1[t]＝r*acos(sinpi(y1/180)*sinpi(y2/180)+cospi(y1/180)*cospi(y2/180)*cospi((x1-x2)/180))；
公式中r表示地球半径；dis1[t]表示aimsi在t时刻与bimsi在该时刻前的最后一个时间点t1的距离；x1、x2分别表示两个loc-ci的经度，y1、y2分别表示两个loc-ci的纬度；
同理可以计算出aimsi在t时刻与bimsi在该时刻后的第一个时间点t2的距离dis2[t]；
同时，计算出两个时间差：
time1[t]＝t-t1；time2[t]＝t2-t；
(34)计算aimsi在t时间点下和bimsi在t1时间点的距离阈值disthr1[t]和最大距离dismax1[t]：
disthr1[t]＝speedthr*time1[t]；
dismax1[t]＝speedmax*time1[t]；
这里speedthr表示速度阈值，设为10公里/小时；speedmax表示最大速度，设为20公里/小时；
同理，计算出aimsi在t时间点下和bimsi在t2时间点的距离阈值disthr2[t]和最大距离dismax2[t]；
(4)按照上述步骤(3)的方法进行循环计算，分别计算出aimsi在各个时间点i下与bimsi最接近的两个时间点的距离dis1[i]和dis2[i]，以及两个距离阈值disthr1[i]和disthr2[i],两个最大距离值dismax1[i]和dismax2[i]；
(5)根据上述计算得到的距离值、距离阈值、和最大距离值，进行一人多号判断：
(51)如果在任意时间点i，距离dis1[i]和dis2[i]分别大于dismax1[i]和dismax2[i]，则直接判断该对imsi为非同一轨迹，即非一人多号；否则，进一步按以下步骤(52)进行置信度判断；
(52)如果距离dis1[i]小于或等于距离阈值disthr1[i]时，按以下公式计算似然值l1[i]：
l1[i]＝1-(dis1[i]/disthr1[i])*α；
如果距离dis1[i]大于距离阈值disthr1[i]，而小于或等于最大距离dismax1[i]时，则似然值l1[i]按以下公式计算：
l1[i]＝(1-α)-(dis1[i]/dismax1[i])*0.1；
上述公式中，α为经验值，默认为α＝0.6；
同理，计算出似然值l2[i]；
计算出一对imsi在各时间点上似然值的均值；
进一步，设定一个阈值β，β＝0.8；
如果一对imsi在各时间点上的似然值的均值小于阈值β，则认为该对imsi为非同一轨迹，认定为非一人多号；
如果一对imsi在各时间点上的似然值的均值大于等于阈值β，则认为该对imsi在所述的时间段内为同一轨迹；设定一个包含有多个连续的所述时间段的时间周期，如果该对imsi在该时间周期内有大于设定数量的时间段被判定为同一轨迹，则认定该对imsi为一人多号。
本发明通过位置间的距离来估算不同的imsi对之间的轨迹相似性，进而对一人多号现象进行判断，较好的实现了一人多号用户的筛选，对利用通讯大数据进行人口统计和监测的做法提供了有利支持。
具体实施方式
本发明的具体实现方式如下：
(1)数据采集：利用电信运营商信令数据，根据每个imsi识别号的轨迹数据，对每个imsi所在的基站位置、进出小区时间信息进行采集并存储。
本发明所采用的数据来自移动运营商的信令数据，包括：用户手机号码-imsi(国际移动用户识别码:internationalmobilesubscriberidentificationnumber)；位置区识别码lac：用于标识不同的位置区；基站编号ci：与位置区识别码(lac)结合，用于识别网络中覆盖的小区；imsi进入基站的时间、离开该基站的时间；
(2)数据筛选：每次读取一个imsi对，设为aimsi和bimsi，作为待确定的一人多号对；在一定时间段内，按时间点先后排序，分别得到这两个imsi的各个时间点以及在各个时间点所对应的基站位置的数据；所述的时间段，一般设定为为1天(24小时)，所述时间点来自于所使用的轨迹数据。
(3)距离计算：根据每个具体基站所对应的经纬度，计算两个imsi在各时间点下的距离，计算方法如下：
(31)先找到aimsi的一个时间点t，再找到bimsi在t时间点前的最后一个时间点t1以及在t时间点后的第一个时间点t2；
(32)对应出aimsi在时间点t所在位置的经纬度，以及bimsi在时间点t1和t2时的位置经纬度；
(33)分别计算出aimsi在t时间点和bimsi在t1时间点时的距离，以及aimsi在t时间点和bimsi在t2时间点时的距离，计算公式如下：
dis1[t]＝r*acos(sinpi(y1/180)*sinpi(y2/180)+cospi(y1/180)*cospi(y2/180)*cospi((x1-x2)/180))；
公式中r表示地球半径；dis1[t]表示aimsi在t时刻与bimsi在该时刻前的最后一个时间点t1的距离；x1、x2分别表示两个loc-ci的经度，y1、y2分别表示两个loc-ci的纬度；
同理可以计算出aimsi在t时刻与bimsi在该时刻后的第一个时间点t2的距离dis2[t]；
同时，计算出两个时间差：
time1[t]＝t-t1；time2[t]＝t2-t；
(34)计算aimsi在t时间点下和bimsi在t1时间点的距离阈值disthr1[t]和最大距离dismax1[t]：
disthr1[t]＝speedthr*time1[t]；
dismax1[t]＝speedmax*time1[t]；
这里speedthr表示速度阈值，设为10公里/小时；speedmax表示最大速度，设为20公里/小时；
同理，计算出aimsi在t时间点下和bimsi在t2时间点的距离阈值disthr2[t]和最大距离dismax2[t]；
(4)按照上述步骤(3)的方法进行循环计算，分别计算出aimsi在各个时间点i下与bimsi最接近的两个时间点的距离dis1[i]和dis2[i]，以及两个距离阈值disthr1[i]和disthr2[i],两个最大距离值dismax1[i]和dismax2[i]；
(5)根据上述计算得到的距离值、距离阈值、和最大距离值，进行一人多号判断：
(51)如果在任意时间点i，距离dis1[i]和dis2[i]分别大于dismax1[i]和dismax2[i]，则直接判断该对imsi为非同一轨迹，即非一人多号；否则，进一步按以下步骤(52)进行置信度判断；
(52)如果距离dis1[i]小于或等于距离阈值disthr1[i]时，按以下公式计算似然值l1[i]：
l1[i]＝1-(dis1[i]/disthr1[i])*α；
如果距离dis1[i]大于距离阈值disthr1[i]，而小于或等于最大距离dismax1[i]时，则似然值l1[i]按以下公式计算：
l1[i]＝(1-α)-(dis1[i]/dismax1[i])*0.1；
上述公式中，α为经验值，默认为α＝0.6；
同理，计算出似然值l2[i]；
计算出一对imsi在各时间点上似然值的均值；
进一步，设定一个阈值β，根据经验值，β取值为0.8；
如果一对imsi在各时间点上的似然值的均值小于阈值β，则认为该对imsi为非同一轨迹，认定为非一人多号；
如果一对imsi在各时间点上的似然值的均值大于等于阈值β，则认为该对imsi在所述的时间段内为同一轨迹；设定一个包含有多个连续的所述时间段的时间周期，如果该对imsi在该时间周期内有大于设定数量的时间段被判定为同一轨迹，则认定该对imsi为一人多号。
上述时间周期以1个月(30天)为例，所述时间段设定为1天(24小时)；如果存在一组imsi，在该月中超过一定天数(设定为10天以上)被认定为同一轨迹，则该组imsi可视为一人多号。