[已解决]代码提速问题！麻烦下(过程已详细补充清楚)

香川群子

已知经纬度两点距离的公式
两点距离=6378137*2*ASIN(SQRT(SUMSQ(SIN((RADIANS(起点纬度) - RADIANS(终点纬度))/2))+ COS(RADIANS(起点纬度))*COS( RADIANS(终点纬度))*SUMSQ(SIN((RADIANS(起点经度)-RADIANS(终点经度))/2))))


或者正确的弧度公式是：
C = sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(lon2-lon1)

需要首先确认一下楼主距离计算公式的准确性。

香川群子

经验近似公式：
=SQRT(((B4-E4)*110965)^2+((C4-F4)*95839)^2)

这个近似公式计算速度会快很多。

直接用经纬度按弧度计算的公式：
=ACOS(SIN(C4)*SIN(F4)+COS(C4)*COS(F4)*COS(E4-B4))*PI()/180*6371004

地球半径和弧度转角度部分的计算=PI()/180*6371004 可固化为常数=111195
=ACOS(SIN(C4)*SIN(F4)+COS(C4)*COS(F4)*COS(E4-B4))*111195

这样的函数改进本身就可以大幅度提高计算速度了。

suxiong

香川群子 发表于 2014-1-8 14:15
你再上一个12000行的数据附件吧。

现在的附件只有50多行用不了多少时间没有意义。

12000.zip (561.77 KB, 下载次数: 4)

2014-1-8 17:31 上传

点击文件名下载附件

这是附件，麻烦你帮我改下可以吗?

suxiong

香川群子 发表于 2014-1-8 16:20
经验近似公式：
=SQRT(((B4-E4)*110965)^2+((C4-F4)*95839)^2)

=SQRT(((B4-E4)*110965)^2+((C4-F4)*95839)^2)

这个算出来也算不多r

香川群子

计算部分已经得到验证：

A列1000行 对B列 7839行、耗时4秒……推算A列12000行  对B列 7839行、耗时大约1分钟。

比你的5个小时要快 300倍!


…………
代码未最后完成，因为A列每一行对应B列距离最小值往往不止1个……

我不知道你需要如何输出这些结果。
是只要任意一个最小值呢？还是要列出全部最小值？或者列出最小的10个值？

香川群子

上附件。

提速300倍。


代码已经做了简单注释，楼主可以自己修改、完善了。

香川群子

在高精度计算中，因为浮点计算的问题，居然还会造成数据类型溢出错误。呵呵。

后来想办法解决了。

香川群子

1. Sub kagawa()
   
2.     Dim i&, j&, k&, l&, m&, k1&, k2&, j1#, w1#, t#, cnt&
   
3.     tms = Timer
   
4.    
   
5.     [e1].CurrentRegion.Sort [f2], 1, [g2], , 1, , , 1
   
6.     'B组数据先按经度、维度排序 以便排除经度差距大距离肯定较远的值，可减少80%计算量
   
7.    
   
8.     m = [e1].End(4).Row - 1 '得到B组数据个数m
   
9.     nr2 = [e2].Resize(m) 'B组序号
   
10.     jr2 = [f2].Resize(m) 'B组经度值
    
11.     wr2 = [g2].Resize(m) 'B组维度值
    
12.    
    
13.     arr = [a1].CurrentRegion '获取A组原始数据
    
14. '    For i = 2 To 1000 '测试时遍历检查A组1000行
    
15.     For i = 2 To UBound(arr) '遍历检查A组各行
    
16.         j1 = arr(i, 2): w1 = arr(i, 3) '经度1、维度1数据存入变量
    
17.         k1 = Application.Match(j1 - 0.02, jr2, 1) '检查比经度1小0.02的起始位置k1
    
18.         k2 = Application.Match(j1 + 0.02, jr2, 1) '检查比经度1大0.02的结束位置k2
    
19.         cnt = cnt + k2 - k1 + 1 '累计本次对比检查B组的次数k2-k2+1
    
20.         
    
21.         ReDim br(10, 1): l = 0 '定义存放10个最小值数据的临时数组br
    
22.         For j = k1 To k2 '遍历B组中k1 to k2 范围 （这样就已经减少了85%左右的计算量）
    
23.             
    
24.             t = f(j1, w1, jr2(j, 1), wr2(j, 1)) '用经验公式f计算两点坐标间距离
    
25.                                                 '这个方法可以提速 5-6倍
    
26.             For k = l - 1 To 0 Step -1
    
27.                 If t > br(k, 1) Then
    
28.                     Exit For
    
29.                 Else
    
30.                     br(k + 1, 1) = br(k, 1)
    
31.                     br(k + 1, 0) = br(k, 0)
    
32.                 End If
    
33.             Next
    
34.             br(k + 1, 1) = t: br(k + 1, 0) = j
    
35.             If l < 10 Then l = l + 1
    
36.             '以上部分为比较法排序筛选得到最小的10个值，以及它们的B组行位置
    
37.         Next
    
38.         For k = 0 To l - 1
    
39.             '仅仅对这10个最小值，用精确公式重新计算其距离值
    
40.             If br(k, 1) Then br(k, 1) = ff(j1, w1, jr2(br(k, 0), 1), wr2(br(k, 0), 1))
    
41.         Next
    
42. '        [i1].Resize(l, 2) = br '结果可以随时输出
    
43.     Next
    
44.     MsgBox Format(Timer - tms, "0.000s ") & cnt '最后统计耗时以及A组B组实际交互检查比对的次数
    
45. End Sub
    
46. Function f(j1, w1, j2, w2) '粗略经验公式
    
47.     f = (((j1 - j2) * 110965) ^ 2 + ((w1 - w2) * 95839) ^ 2) ^ 0.5
    
48. End Function
    
49. Function ff(j1, w1, j2, w2) '精确计算公式
    
50. '    t = Sin(w1) * Sin(w2) + Cos(w1) * Cos(w2) * Cos(j1 - j2)
    
51. '    If Round(t, 12) = 1 Then ff = 0 Else ff = Application.Acos(t) * 111195
    
52.     ff = Application.Acos(Sin(w1) * Sin(w2) + Cos(w1) * Cos(w2) * Cos(j1 - j2)) * 111195
    
53. End Function

复制代码

香川群子

如果A组、B组个数不同，那么应该以数量多的做B组，数量少的做A组。

这样计算效率更高。
因为B组超范围的数据会自动被舍掉不用计算，
那么A组的个数是必须检查的，所以A组数量应该设置的更少一些。

…………
当然，如果楼主的需求就是需要检查数量多的A组的结果，那么就没办法交换A组B组位置了。

香川群子

如果A组、B组个数不同，那么应该以数量多的做B组，数量少的做A组。

这样计算效率更高。
因为B组超范围的数据会自动被舍掉不用计算，
那么A组的个数是必须检查的，所以A组数量应该设置的更少一些。

…………
当然，如果楼主的需求就是需要检查数量多的A组的结果，那么就没办法交换A组B组位置了。