Modifying  @ruebenko's answer in the `StackOverflow` Q/A linked in Faysal's comment 
[(Mathematica fast 2D binning algorithm)][1] to get the means of `z`-values for each bin (using yet another undocumented setting for the option `"TreatRepeatedEntries"` ):

     zvalues = data[[All, 3]];
     epsilon = 1*^-10;
     indexes = 1 + Floor[(1 - epsilon) 256 Rescale[data[[All, {1, 2}]]]];
     System`SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> (Mean[{#}] &)}];
     binmeansZ = SparseArray[indexes -> zvalues];
     System`SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> First}];

A picture:

     MatrixPlot[binmeansZ]

![enter image description here][2]

Timing:

     Timing[System`SetSystemOptions["SparseArrayOptions" ->
        {"TreatRepeatedEntries" -> (Mean[{#}] &)}]; 
     binmeansZ = SparseArray[indexes -> zvalues];
     System`SetSystemOptions[
       "SparseArrayOptions" -> {"TreatRepeatedEntries" -> First}];]
     (* {0.031200, Null} *)


  [1]: http://stackoverflow.com/q/8178714/884752
  [2]: https://i.sstatic.net/jZAP0.png